#!/usr/bin/perl
#
my @tau;
my @adev;
my @mdev;
my @tdev;

# Allan($tau,\@data) fills global arrays tau,adev,mdev,tdev
sub Allan {
	my $modified = 0;
	my $x;
	my $d;
	my @data;
	my $dd;
	my $avg;
	my $s;
	my @y;
	my $ix;
	my $lix;
	my $n;
	my $k;
	my $l;
	my $m;

	my ($tau0,@y) = @_;
	my $ndata = $#y;

	$lix = 0;
	for ($x = 1.0; $x < $ndata/3; $x *= 1.1) {
		$ix = $x + .5;
#		if ($ix == $lix) 
#			continue;
		$lix = $ix;
		$m = $ix;
		$s = 0;
		$n = 0;
		$dd = 0;
		for ($k=0;$k<$ndata-2*$ix;$k++) {
			$y[$k]=($data[$k]+$data[$k+2*$ix]-2*$data[$k+$ix])/$tau;
			}
		if ($modified) {
			$d = 0;
			for ($l = 0; $l < $ix; $l++) { 
				$d += $y[l];
				}
			$s += $d * $d;
			for ($k=0;$k<$ndata-3*$ix;$k++) {
				$d -= $y[$k];
				$d += $y[$k+$ix];
				$s += $d * $d;
				$n++;
				}
			$avg = $s/(2.0*($ndata-3*$m+1)*$m*$m*$m*$m);
			}
		else {
			for ($k=0;$k<$ndata-2*$ix;$k+=$ix) {
				$s += $y[$k] * $y[$k];
				$n++;
				}
			$avg = $s/(2*$n*$m*$m);
			}
	if ($n>1) {
		printf("%g %e\n",$ix * $tau, sqrt($avg));
	}
}

# Add data points from stdin
while ($reading=<STDIN>) {
	if (substr($reading,0,1) =~ "[0-9]") {
		($x,$y) = split(/\s/,$reading);
		push(@xarray,$x);
		push(@yarray,$y);
		}
	}
avar(300,@yarray);
