#!/usr/bin/perl -w

# aging.pl
# 4 November 2006
#
# Copyright 2005 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 IO::File;
use Getopt::Std;
use n8ur qw(round);

my $opt;
my $opt_string;
my $data_file;
my $reading;
my $tau;
my $interval;
my $verbose;
my $df;
my $x;
my $y;
my @xarray;
my @yarray;
my $offset;
my $last_offset;
my $counter;
my $pointer;
my $last_pointer;
my $i;
my $aging;
my $start;
my $finish;

#----------
# usage and option initialization
$opt_string = 'h:f:t:r:vi:';
sub usage() {
print STDERR << "EOF";

usage: $0 [-hv] -t tau -i interval -r readings for offset -f filename

-h	: this (help) message

-t	: tau in teconds; if not specified, will be calculated from data

-i	: number of readings between calculations; default is 1 day

-f	: data file name; if "-" read from STDIN

-v	: verbose output

EOF
}

#----------------------
getopts( "$opt_string", \my %opt ) or usage() and exit;
usage() and exit if $opt{h};
usage() and exit if !$opt{f};

$tau = 0;
if ($opt{t}) {
	$tau = $opt{a};
	}

$interval = 0;
if ($opt{i}) {
	$interval = $opt{i};
	}

$verbose = 0;
if ($opt{v}) {
	$verbose = 1;
	}

if ($opt{f}) {
	$data_file = $opt{f};
	}
else {
	die "No data file specified!\n";
	}

###############################################################################
# main program
# #############################################################################

# Add data points
$df = new IO::File $data_file, "r";
if (defined $df) {
	while ($reading = <$df>) {
		if (substr($reading,0,1) =~ "[0-9]") {
			($x,$y) = split(/\s/,$reading);
			push(@xarray,$x);
			push(@yarray,$y);
		}
	}
	$df->close;
} else {
	exit;
	}

# calculate tau if necessary -- average from first ten MJD values
# trying to average over full file instead
if ($tau == 0) {
	$tau = (($xarray[$#xarray]-$xarray[0])/($#xarray-1))* 86400;
	if ($tau >= 500) {
		# long taus tend to round up, so truncate instead
		$tau = int($tau);
		}
	if (($tau < 500) && ($tau >= 10)) {
		$tau = round(0,$tau);
		}
	if ($tau < 10) {
		$tau = round(1,$tau);
		}
}

# set interval to 1 day if not set on command line
if ($interval == 0) {
	$interval = 86400/$tau;
	}

$pointer = 0;
$last_pointer = 0;
$offset = 0;

# loop so long as we have enough data left to do the offset measurement
while (($pointer + $interval) < ($#xarray)) {

	# are we at a measurement point?
	if ( ($pointer == 0) || ($pointer == ($last_pointer + $interval) ) ) {
		$last_pointer = $pointer;

		# calculate the end-point offset over $offset_readings values
		# starting at current array position
		$start = $yarray[$pointer];
		$finish = $yarray[$pointer + $interval];

		# save $last_offset
		$last_offset = $offset;

		# calculate the offset -- delta t over t
		$offset = ($finish - $start)/($interval * $tau);
		
		# if we have a valid previous offset, subtract,
		# normalize to 24 hour value, and print
		if ($last_offset != 0) {
			if ($verbose) {
			printf "Offset: %.4e    Aging/day: %.4e\n", $offset,
			 ($offset-$last_offset)*(86400/($interval*$tau));
			}
			else {
				printf "%.4e\n",
			 	($offset-$last_offset)*(86400/($interval*$tau));
			}
		}
		else {
			if ($verbose) {
				printf "Offset: %.4e\n", $offset,
			}
		}
		
	}		
	$pointer++;
}
exit 0;