#!/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;
use GD;
use GD::Text;
use GD::Text::Align;
use GD::Graph::mixed;
use GD::Graph::colour;


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

# Left and center footers; there's also a timestamp in the lower right corner
my $left_footer = "jra\@febo.com";
my $center_footer = "";

open (PIC1, ">test-phase.png") ||
        die "Can't open image file!\n";
open (PIC2, ">test-sigma.png") ||
        die "Can't open image file!\n";


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;
my $w;
my $h;
my @xsigma;
my @ysigma;
my $avar;
my $l10 = log(10);

# 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) {
	$avar = &adev($tau,$sample_size,$i,@yarray);
	printf "%4.3E\t%4.3E\n",$tau*$i,$avar;
	push(@xsigma,$tau*$i);
	# convert to log format for plotting
	push(@ysigma,log($avar)/$l10);
	}	

# normalize y values
for ($i=0; $i<$sample_size; $i++) {
	$yarray[$i] = ($yarray[$i] -= $ymin)*10e5;
	}

# plotting routines

# set some dimensions
my $l_margin = 30;
my $r_margin = 25;
my $t_margin = 50;
my $b_margin = 35;
$x = 650;
$y = 500;

my $plot = new GD::Graph::mixed($x,$y);

$plot->set (
	l_margin => $l_margin,
	r_margin => $r_margin,
	t_margin => $t_margin,
	b_margin => $b_margin,
	transparent => 0,
	y_tick_number => 10,
	y_long_ticks => 1,
	x_tick_number => 10,
	x_label_skip => 2,
	x_long_ticks => 1
	);

my $foo;
$foo = $plot->plot([\@xarray,\@yarray]);

my $black = $foo->colorAllocate(0,0,0);
my $red = $foo->colorAllocate(255,0,0);

#-----------
# add text and number formatting features
my $gd_text = GD::Text->new() or die GD::Text::error();

my $title = "HP 5065A vs. GPS";
$gd_text->set_font(gdGiantFont);
$gd_text->set_text($title);
($w, $h) = $gd_text->get('width', 'height'); $foo->string(gdGiantFont,($x/2)-($w/2),10,$title,$black);

my $annotation1 = "via UT+ and HP 5370B";
$gd_text->set_font(gdLargeFont);
$gd_text->set_text($annotation1);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdLargeFont,($x/2)-($w/2),27,$annotation1,$black);

my $xlabel = "MJD";
$gd_text->set_font(gdSmallFont);
$gd_text->set_text($xlabel);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdSmallFont,($x/2)-($w/2),$y-30,$xlabel,$black);

my $ylabel = "Microseconds";
my $align = GD::Text::Align->new($foo);
$align->set_font(gdSmallFont);
$align->set_text($ylabel);
$align->set (color=>$black);
$align->draw(20, 286, 3.1416/2);

$gd_text->set_font(gdTinyFont);
$gd_text->set_text($left_footer);
$foo->string(gdTinyFont,20,$y-10,$left_footer,$black);

$gd_text->set_font(gdTinyFont);
$gd_text->set_text($center_footer);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdTinyFont,($x/2)-($w/2)+($l_margin/2),$y-10,
        $center_footer,$black);

my $timestamp;
$timestamp = gmtime() . " (UTC)";
$gd_text->set_font(gdTinyFont);
$gd_text->set_text($timestamp);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdTinyFont,$x-$w-12,$y-10,$timestamp,$black);

binmode PIC1;
print PIC1 $foo->png;
close PIC1;
################################################################################
# plotting routines for ADEV
################################################################################

sub yform {
	my $y = shift;
	print $y,"\n";
	return int(10**$y+.5)*1e6;
	}

sub Max {
# takes an array ref - returns the max
# $max = Max(\@TempRand) or $max = Max([2, 199, 238]);
    
    my $list = shift;
    my $max = $list->[0];
    foreach (@$list) {
          $max = $_ if ($_ > $max);
          # (to do min) $min = $_ if ($_ < $min);
         }
    
    return($max);
   }

my $ym = int(log(abs(Max(\@ysigma)))/$l10+.99999);

# set some dimensions
$l_margin = 30;
$r_margin = 25;
$t_margin = 50;
$b_margin = 35;
$x = 650;
$y = 500;

$plot = new GD::Graph::mixed($x,$y);

$plot->set (
	l_margin => $l_margin,
	r_margin => $r_margin,
	t_margin => $t_margin,
	b_margin => $b_margin,
	transparent => 0,
	y_max_value => -10,
	y_min_value => -13,
	y_tick_number => 10, 
	y_number_format => \&yform,
	y_long_ticks => 1,
	x_tick_number => 10,
	x_label_skip => 2,
	x_long_ticks => 1
	);

$foo = $plot->plot([\@xsigma,\@ysigma]);

$black = $foo->colorAllocate(0,0,0);
$red = $foo->colorAllocate(255,0,0);

#-----------
# add text and number formatting features
$gd_text = GD::Text->new() or die GD::Text::error();

$title = "HP 5065A vs. GPS";
$gd_text->set_font(gdGiantFont);
$gd_text->set_text($title);
($w, $h) = $gd_text->get('width', 'height'); $foo->string(gdGiantFont,($x/2)-($w/2),10,$title,$black);

$annotation1 = "via UT+ and HP 5370B";
$gd_text->set_font(gdLargeFont);
$gd_text->set_text($annotation1);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdLargeFont,($x/2)-($w/2),27,$annotation1,$black);

$xlabel = "Tau";
$gd_text->set_font(gdSmallFont);
$gd_text->set_text($xlabel);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdSmallFont,($x/2)-($w/2),$y-30,$xlabel,$black);

$ylabel = "Allen Variance";
$align = GD::Text::Align->new($foo);
$align->set_font(gdSmallFont);
$align->set_text($ylabel);
$align->set (color=>$black);
$align->draw(20, 286, 3.1416/2);

$gd_text->set_font(gdTinyFont);
$gd_text->set_text($left_footer);
$foo->string(gdTinyFont,20,$y-10,$left_footer,$black);

$gd_text->set_font(gdTinyFont);
$gd_text->set_text($center_footer);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdTinyFont,($x/2)-($w/2)+($l_margin/2),$y-10,
        $center_footer,$black);

$timestamp = gmtime() . " (UTC)";
$gd_text->set_font(gdTinyFont);
$gd_text->set_text($timestamp);
($w, $h) = $gd_text->get('width', 'height');
$foo->string(gdTinyFont,$x-$w-12,$y-10,$timestamp,$black);

binmode PIC2;
print PIC2 $foo->png;
close PIC2;

exit 0;