#!/usr/bin/perl -w
##
## stable-stat.pl
## version 0.1 -- generate statistics from file containing MJD and phase data
## 31 December 2004
##
## Copyright 2004 by John R. Ackermann  N8UR (jra@febo.com)
## Licensed under the GPL version 2 or later; see the file COPYING
## included with this distribution.  I request, but do not require, that
## any modifications that correct bugs or errors, or increase the program's
## functionality, be sent via email to the author at the address above.

use strict;
use Statistics::OLS;

# This needs to be made a command line argument
my $tau = 300;

my $i;
my $x;
my $y;
my $reading;
my @xarray;
my @yarray;
my $days;
my $seconds;
my $total_offset;
my $linear_offset;
my $sample_size;
my $ls;
my $intercept;
my $slope;
my $R_squared;
my $xmin;
my $xmax;
my $ymin;
my $ymax;

# ADEV calculation
sub adev {
	my ($tau,$sample_size,$avg,@yarray) = @_;
	my $sum = 0.0;
	my $i;
	my $y;
	my $last_y;
	my $dy;
	my $samples = 1;
	my $adev;

	if ($sample_size > 2) {
		for ($i=$avg; $i<$sample_size; $i = $i + $avg) {
			$y = $yarray[$i] - $yarray[$i-$avg];
			if ($i > $avg) {
				$dy = $y - $last_y;
				$sum += $dy * $dy;
				}
			$last_y = $y;
			$samples++;
			}
	}
	return $adev = (sqrt($sum / (2.0 * ($samples - 1))))/($tau*$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);
		}
	}

# set up linear regression	
$ls = Statistics::OLS->new(); 
$ls->setData (\@xarray, \@yarray);
$ls->regress();

# get useful information
($intercept, $slope) = $ls->coefficients();
$R_squared = $ls->rsq();
$sample_size = $ls->size();    
($xmin, $xmax, $ymin, $ymax) = $ls->minMax();

# do some calculations
$total_offset = $yarray[$sample_size-1]-$yarray[0];
$days = $xmax-$xmin;
$seconds = $days*86400;

# For some reason, the returned slope seems to be per day, not per sample.
# This could potentially break with differently formatted data.
$linear_offset = $slope/86400;

# output routine
print "\n";
printf "%i samples at tau %s seconds over %.2f days (%i missing samples)\n",
	$sample_size,$tau,$days,($seconds/$tau)-$sample_size;
print "\n";
printf "End-to-End drift:\t%4.3E\n",$total_offset;
printf "Drift per day:\t\t%4.3E\n",$total_offset/$days;
printf "Peak-to-Peak drift:\t%4.3E\n",abs($ymin-$ymax);
print "\n";
printf "Frequency Offset (linear):\t%4.3E\n",$linear_offset;
printf "Frequency Offset (endpoints):\t%4.3E\n",$total_offset/$seconds;
print "\n";

# Do ADEV calculations and output
print "Tau\t\tADEV\n";
for ($i=1; $i<($sample_size/3); $i = $i*2) {
	printf "%4.3E\t%4.3E\n",
		$tau*$i,&adev($tau,$sample_size,$i,@yarray);
	}