#!/usr/bin/perl

use lib $ENV{HOME}.'/lib';
use lib $ENV{HOME}.'/lib/perl';

use Param;

my $p = Param::Parse
(
	{'info','simple time series calculation. I create a signal out of the given function (depending on time $n and/or recursively the last value) and add noise (gaussian or equal distributed) to it if you wish'},
	'length',100,'n','length of time series (in steps)',
	'function','$t','f','straight signal ($t is the time, $x last value)',
	'noise',0,'N','factor / standard deviation for added noise',
	'gaussvar',1,'g','use gaussian (normal distributed) noise (0: just use equal distributed noise)',
	'startx',0,'x','start value for last x for recursion (meaningless otherwise)',
	'withstart',0,'w','include the start x',
	'step',1,'s','step to increment time (sampling)',
	'starttime',0,'t','starting time',
	'withtime',0,'W','include a column showing the time'
);

die "step > 0!" unless $p->{step} > 0;
my $noise = $p->{gaussvar} > 0 ? 'Gauss()' : '(rand()-0.5)';
my $sum = 0;
my @list = ();
my $fs = 'sub { my ($t,$x) = @_; return '.$p->{function}.($p->{noise} ? ' + '.$p->{noise}.' * '.$noise : '').'; }';
print STDERR "time series using function: $fs\n";
my $f = eval $fs;
my $x = $p->{startx};
print STDOUT ($p->{withtime} ? $p->{starttime}."\t" : '').$x."\n" if $p->{withstart};

for(my $n = 1; $n <= $p->{length}; ++$n)
{
	$x = &$f($p->{starttime}+$n*$p->{step},$x);
	print $p->{starttime}+$n*$p->{step}."\t" if $p->{withtime};
	print $x."\n";
}

my $cache = undef;
sub Gauss
{
	#Polar-Methode
	my $x;
	if(defined $cache)
	{
		$x = $cache;
		undef $cache;
	}
	else
	{
		my $u1,$u2;
		my $v;
		do
		{
			$u1 = rand();
			$u2 = rand();
			$v = (2*$u1-1)**2+(2*$u2-1)**2;
		}
		while($v >= 1);
		$cache = (2*$u2-1)*sqrt(-2*log($v)/$v);
		$x = (2*$u1-1)*sqrt(-2*log($v)/$v);
	}
	return $x;
}