#!/usr/local/bin/perl

use strict;
use warnings;
use Time::HiRes qw(time usleep);
use POSIX qw(strftime);
use POSIX qw(floor);
use Term::ANSIColor qw(:constants);
use Getopt::Long qw(GetOptions);
use constant {N => "\e[0K\n", T => "\e[36G"};

my %prefix = (
	'-8' => 'y',    '8' => 'Y',
	'-7' => 'z',    '7' => 'Z',
	'-6' => 'a',    '6' => 'E',
	'-5' => 'f',    '5' => 'P',
	'-4' => 'p',    '4' => 'T',
	'-3' => 'n',    '3' => 'G',
	'-2' => 'µ',    '2' => 'M',
	'-1' => 'm',    '1' => 'k',
	'0' => '' ,
);

sub average {
	my($data) = @_;
	if (not @$data) {
		die("Empty array");
	}
	my $total = 0;
	foreach (@$data) {
		$total += $_;
	}
	my $average = $total / @$data;
	return $average;
}

sub stdev {
	my($data) = @_;
	if(@$data == 1){
		return 0;
	}
	my $average = &average($data);
	my $sqtotal = 0;
	foreach(@$data) {
		$sqtotal += ($average-$_) ** 2;
	}
	my $std = ($sqtotal / (@$data-1)) ** 0.5;
	return $std;
}

sub absmax {
	my($data) = @_;
	if(@$data == 1){
		return 0;
	}
	my $absmax=0;
	foreach(@$data) {
		if(abs($_)>$absmax) {
			$absmax=abs($_);
		}
	}
	return $absmax;
}

sub absmin {
	my($data) = @_;
	if(@$data == 1){
		return 0;
	}
	my $absmin=1e15;
	foreach(@$data) {
		if(abs($_)<$absmin) {
			$absmin=abs($_);
		}
	}
	return $absmin;
}

sub rms0 {
	my($data) = @_;
	if(@$data == 1) {
		return 0;
	}
	my $sqtotal = 0;
	foreach(@$data) {
		$sqtotal += $_ ** 2;
	}
	my $std = ($sqtotal / (@$data-1)) ** 0.5;
	return $std;
}

sub print_sec {
	my $num = shift;
	if($num == 0) {
		return '0 s';
	}
	my $sign = ($num < 0) ? '-' : '';
	$num = abs $num;
	my $e = floor(log($num) / log(1000));
	my $mult = 1000**$e;
	$num = 0 + sprintf '%.5g', ($num / $mult);
	return $sign . $num . ' ' . $prefix{$e} . 's';
}

sub printms {
	my $num = shift;
	print_sec($num/1000);
}

my $n=0;
my $delay = 0;
my $gga = 0;

GetOptions (
	"delay=i" => \$delay,
	"n=i" => \$n,
	"gga" => \$gga,
) or die ("Error in command line arguments\n");

if($delay) {
	printf "\e[?25l\e[2J"; # hide cursor & clear screen
}

do {
	if($delay) {
		printf "\e[H";
	}
	my $date = strftime "%FT%T", localtime;
	print REVERSE "GPSDO NTP time server status at $date (local)\n", RESET;

	$_ = `ntpq -np`;
	s/\n\*/\n\e[1;37;42m\*/g;
	s/\no/\n\e[1;37;42mo/g;
	s/\n\+/\n\e[0;30;43m\+/g;
	s/\n-/\n\e[0;37;41m\-/g;
	s/\n/\e[0m\e[0K\n/g;
	print $_;

	my $r1 = qr/maximum error:\s+(\S+)\n.*estimated error:\s+(\S+)\n/s;
	my $r2 = qr/stratum=(\S+),.*refid=(\S+),.*offset=(\S+),.*sys_jitter=(\S+),/s;
	my $r3 = qr/up (.*),\s+\S+\s+user.*load averages:\s+(.*)\n/s;

	my (@p) = (`ntpq -c kerninfo` =~ /$r1/);

	print N."Maximum error:".T.printms($p[0]).N;
	print "Estimated error:".T.$p[1].N;

	my (@q) = (`ntpq -crv` =~ /$r2/);
	print "Stratum:".T.$q[0].N;
	print "Reference ID:".T.$q[1].N;

	print BOLD UNDERLINE "Combined time offset:".T.printms($q[2]), RESET, N;
	print UNDERLINE "Combined system jitter:".T.printms($q[3]), RESET, N;

	$_ = `ntpq -ccl`;
	/.*timecode="(.*?)".*fudgetime2=(.*?),/gs;
	print "Fudgetime2:".T.$2.N;
	print "Timecode:".N.$1.N;

# $GPGGA,122800.00,5421.11103,N,01843.31182,E,1,09,1.00,3.1,M,32.6,M,,*5F

	if($gga) {
		my @nmea = split ",", $1;

		my $lat = $nmea[2];
		my $latdegmin = floor($lat);
		my $latdeg = floor($latdegmin/100);
		my $latmin = floor($latdegmin%100);
		my $latsec = 0 + sprintf("%.5g", 60*($lat-$latdegmin));
		my $lon = $nmea[4];
		my $londegmin = floor($lon);
		my $londeg = floor($londegmin/100);
		my $lonmin = floor($londegmin%100);
		my $lonsec = 0 + sprintf("%.5g", 60*($lon-$londegmin));
		
		print "Latitude:".T."$latdeg° $latmin' $latsec\" $nmea[3]".N;
		print "Longitude:".T."$londeg° $lonmin' $lonsec\" $nmea[5]".N;
		print "Altitude:".T.$nmea[9]." ".$nmea[10].N;
		print "Fix:".T.$nmea[6].N;
		print "Number of satellites:".T.$nmea[7].N;
		print "HDOP:".T.$nmea[8].N;
	}

	$_ = `ntpq -csysstats`;
	/.*uptime:.*?([0-9]+)/g;

	print "NTPd uptime:".T;
	my $ntpu = $1;
	if(floor($ntpu/86400)) {
		print floor($ntpu/86400)." days ";
	}
	$ntpu %= 86400;
	if(floor($ntpu/3600)) {
		print floor($ntpu/3600)." hours ";
	}
	$ntpu %= 3600;
	if(floor($ntpu/60)) {
		print floor($ntpu/60)." minutes ";
	}
	$ntpu %= 60;
	print $ntpu . " seconds".N;

	my (@u) = (`uptime` =~ /$r3/);
	print "System uptime:".T.$u[0].N;
	print "Load average (1, 5, 15 min):".T.$u[1].N;

	$_ = `sysctl -a | grep temperature`;
	s/hw.*: (.*)C\nhw.*: (.*)C/$1°, $2°C/gm;

	print "CPU temperatures:".T.$_;

	my $peerstats = "";
	if($n) {
		$peerstats = `tail -$n /var/log/ntpstats/peerstats`;
	} else {
		$peerstats = `cat /var/log/ntpstats/peerstats`;
	}
	my @gpsoffs;
	foreach my $line (split /[\r\n]+/, $peerstats) {
		if ($line =~ /127\.127\.20\.0/) {
			my @fields = split " ", $line;
			push(@gpsoffs, $fields[4]);
		}
	}
	print "GPS (NMEA) offset average (n=$n):".T.print_sec(average(\@gpsoffs)).N;
	print "GPS (NMEA) offset stdev (n=$n):".T.print_sec(stdev(\@gpsoffs)).N;

	my $loopstats = "";
	if($n) {
		$loopstats = `tail -$n /var/log/ntpstats/loopstats`;
	} else {
		$loopstats = `cat /var/log/ntpstats/loopstats`;
	}

	my @clockoffs;
	my @clockjitt;
	foreach my $line (split /[\r\n]+/, $loopstats) {
		my @fields = split " ", $line;
		push(@clockoffs, $fields[2]);
		push(@clockjitt, $fields[4]);
	}
	print "Clock offset average (n=$n):".T.print_sec(average(\@clockoffs)).N;
	print "Clock offset stdev (n=$n):".T.print_sec(stdev(\@clockoffs)).N;
	print UNDERLINE "Clock offset RMS (n=$n):".T.print_sec(rms0(\@clockoffs)), RESET, N;
	print "Clock offset absmax (n=$n):".T.print_sec(absmax(\@clockoffs)).N;
	print "Clock offset absmin (n=$n):".T.print_sec(absmin(\@clockoffs)).N;

	print UNDERLINE "Clock jitter average (n=$n):".T.print_sec(average(\@clockjitt)).RESET.N;
	print "Clock jitter stdev (n=$n):".T.print_sec(stdev(\@clockjitt)).N;
	print "Clock jitter absmax (n=$n):".T.print_sec(absmax(\@clockjitt)).N;
	print "Clock jitter absmin (n=$n):".T.print_sec(absmin(\@clockjitt)).N;

	if($delay) {
		printf "\e[0m\e[0K";
		sleep $delay;
	}
} while($delay);