#!/usr/bin/perl -w

use strict;
use Getopt::Long;
use Carp;
use POSIX;

sub turn2str ($;$);
sub readheader($);
sub mjd2cal($);
sub cal2mjd($$$;$);

my ($invalid, $power, $cal, $complex, $seconds, $frame, $version,
    $nbits, $refepoch, $representation, $antid, $sampleblocklength,
    $nchan, $threadid, $groupid, $secondaryid, $period, $numsamples,
    $sync, $framelength, $protocol, $metaid);

my $once = 0;
my $check = 0;
my $dosync = 0;
my $help = 0;
my $skip = 0;
my $headeronly = 0;
my $summary = 0;

GetOptions('once'=>\$once, 'check'=>\$check, 'help'=>\$help, 'skip=i'=>\$skip, 'sync'=>\$dosync,
	   'headeronly'=>\$headeronly, 'summary'=>\$summary);

my ($lastframe, $lastsec);
my %lastframe = ();

if ($help || @ARGV==0) {
  print<<EOF;
Usage: codifheader.pl [options] <codiffile>

Options:
   -once          Print only the first header 
   -check         Do check frames increase monotonically with no gaps (single thread only?)
   -sync          Check SYNC word is correct and report byte offset if not
   -skip <bytes>  Skip <bytes> bytes at the start of each file
   -summary       Just print summary of headers
   -headeronly    Read a file which has just CODIF headers, no data segment
EOF
}

foreach (@ARGV) {
  open(CODIF, $_) || die "Could not open $_: $!\n";

  print "Reading $_\n\n";

  my $first = 1;
  %lastframe = ();
  while (1) {

    if ($skip>0) {
      my $status = sysseek CODIF, $skip, SEEK_SET;
      if (! defined $status) {
	warn "Error trying to skip $skip bytes - skipping file\n";
	last;
      } elsif ($status != $skip) {
	warn "Failed to skip $skip bytes - is this a CODIF file?\n";
	last;
      }
      $skip=0;
    }

    ($frame, $seconds, $refepoch, $nbits, $power, $invalid, $complex,
     $cal, $representation, $version, $protocol, $period, $threadid,
     $groupid, $secondaryid, $antid, $nchan, $sampleblocklength, $framelength,
     $numsamples, $sync, $metaid) = readheader(*CODIF);

    if (!defined $invalid) {
      print "   empty file\n" if ($first);
      close(CODIF);
      last;
    }

    my $hexsync = sprintf("0x%X", $sync);
    my $timestr = turn2str(fmod($seconds/60/60/24, 1.0));
    # Derived parameters
    my $samplesperframe;
    if ($nbits==0 || $nchan==0) {
      $samplesperframe = '--';
    } else {
      $samplesperframe = $framelength/($nbits*$nchan*($complex+1))*8;
    }
    my $frameperperiod;
    my $frameuSec;
    my $secoffset;
    my $timestr2;
    if ($samplesperframe eq '--') {
       $frameperperiod = '--';
       $frameuSec = '--';
       $secoffset = '--';
       $timestr2 = '--';
    } else {
       $frameperperiod = $numsamples/$samplesperframe;
       $frameuSec = $period/$frameperperiod*1e6;
       $secoffset = $frame*($frameuSec/1e6);
       $timestr2 = turn2str(fmod(($seconds+$secoffset)/60/60/24, 1.0),9);
    }

    if (!$check && !$dosync) {
      print "-------------------\n" if (!($first | $summary));
      my $mjd =  cal2mjd(1, ($refepoch%2)*6+1, 2000 + int($refepoch/2));
      $mjd += int($seconds/(60*60*24));
      my ($day, $month, $year, $ut) = mjd2cal($mjd);
      my $date = sprintf("%02d/%02d/%04d", $day, $month, $year);

    
      if ($summary) {
	  printf("%8d/%6d $timestr2 %2d %3d\n", $seconds, $frame, $threadid, $groupid);
      } else {
	  printf("FRAME#:  %8d  $timestr2\n", $frame);
	  printf("SECONDS: %8d  $timestr   $date\n", $seconds);
	  print<<EOF;

EPOCH:       $refepoch
NBITS:       $nbits
POWER:       $power
INVALID:     $invalid
COMPLEX:     $complex
CALON:       $cal
REPR:        $representation
VERSION:     $version
PROTOCOL:    $protocol
PERIOD:      $period

THREADID:    $threadid
GROUPID:     $groupid
SECONDARYID: $secondaryid
STATIONID:   $antid

NCHAN:       $nchan
SAMPLEBLOCK: $sampleblocklength
DATALENGTH:  $framelength

#SAMPLES:    $numsamples

SYNC:        $hexsync
METADATA ID: $metaid

Frame uSec = $frameuSec
Samples/Frame = $samplesperframe

EOF
      }
    } elsif ($check) {
      if ($sync!=0xFEEDCAFE) {
	print("Skippig bad sync ($hexsync)\n");
	next;
      }
      my $thisID = "$secondaryid-$groupid-$threadid";
      if (!exists $lastframe{$thisID}) {
	$lastframe{$thisID} = [$seconds, $frame];
      } else {
	my $lastsec = $lastframe{$thisID}->[0];
	my $lastframeno = $lastframe{$thisID}->[1];
	if ($frame-$lastframeno!=1) {
	  if ($seconds==$lastsec) {
	    printf("Skipped %d frames at $seconds $timestr/$lastframeno--$frame $timestr2:  $thisID\n", $frame-$lastframeno-1);
	  } elsif ($seconds-$lastsec==$period) {
	    if ($frame!=0) {
	      printf("Skipped > %d frames at $seconds $timestr/$lastframeno--$frame $timestr2:  $thisID\n", $frame);
	    }
	  } else {
	      printf("Skipped  %d seconds at $timestr2 (sync=$hexsync):  $thisID\n", $seconds-$lastsec);
	  }
	}
	$lastframe{$thisID} = [$seconds, $frame];
      }
    } elsif ($dosync) {
      if ($sync!=0xFEEDCAFE) {
	my $pos = sysseek(CODIF, 0, 1) - 64; # Account for header
	die "Bad sync ($hexsync) at offset $pos\n";
      }
    }

    if (!$headeronly) {
      my $status = sysseek(CODIF, $framelength, 1);
      if (!defined $status) {
	close(CODIF);
	last;
      }
    }
    last if ($once);
    $first = 0;
  }

  print "All good!\n" if ($dosync || $check);
}

sub readheader($) {
  my $codif = shift;
  my $buf;
  my $nread = sysread($codif, $buf, 64);
  if (! defined($nread)) {
    die("Error reading $_: $!");
    return undef;
  } elsif ($nread==0) { # EOF
    return undef;
  } elsif ($nread!=64) {
    die "Error: Only read $nread bytes from header\n";
  }

  my ($word0, $word1, $word2, $word3, $numbersamples, $word5, $meta1, $meta2) = unpack 'QQQQQQQQ', $buf;

  # Word 0
  my $seconds = ($word0>>32)&0xFFFFFFFF;
  $frame = $word0 & 0xFFFFFFFF;
  
  # Word 1
  # Byte 0
  my $refepoch = $word1 & 0xFF;
  # Byte 1
  my $nbits = ($word1>>8)&0xFF;
  # Byte 2
  my $power = ($word1>>16)&0x1;
  my $invalid = ($word1>>17)&0x1;
  my $complex = ($word1>>18)&0x1;
  my $cal = ($word1>>19)&0x1;
  my $representation = ($word1>>20)&0xF;
  # Byte 3
  my $version = ($word1>>24)&0x1F;
  my $protocol = ($word1>>29)&0x7;
  # Byte 4-5 unused
  my $period = ($word1>>48)&0xFFFF;
  
  # Word 2
  my $threadid = $word2 & 0xFFFF;
  my $groupid = ($word2>>16) & 0xFFFF;
  my $secondaryid = ($word2>>32) & 0xFFFF;
  my $antid = unpack 'A2', pack('n', ($word2>>48)&0xFFFF);
  
  # Word 3
  my $nchan = $word3 & 0xFFFF;
  my $sampleblocklength = ($word3>>16)&0xFF;
  my $framelength = ($word3>>32)&0xFFFFFFFF;
  $framelength *= 8;

  # Word 4 - done directly

  # Word 5
  my $sync = $word5 & 0xFFFFFFFF;
  my $metaid = ($word5>>32) & 0xFFFF;
  
  return ($frame, $seconds, $refepoch, $nbits, $power, $invalid,
	  $complex, $cal, $representation, $version, $protocol,
	  $period, $threadid, $groupid, $secondaryid, $antid, $nchan,
	  $sampleblocklength, $framelength, $numbersamples, $sync, $metaid);
}

sub turn2str ($;$) {
  my($turn) = @_;
  my $mode = 'H';
  if (($mode ne 'H') && ($mode ne 'D')) {
    carp 'turn2str: $mode must equal \'H\' or \'D\'';
    return undef;
  }
  my $strsep = ':';

  my ($angle, $str, $sign, $wholesec, $secfract, $min);

  if ($mode eq 'H') {
    $angle = $turn * 24;
  } else {
    $angle = $turn * 360;
  }

  if ($angle < 0.0) {
    $sign = -1;
    $angle = -$angle;
  } else {
    $sign = 1;
  }

  my $wholeangle = (int $angle);

  $angle -= $wholeangle;
  $angle *= 3600;

  # Get second fraction
  $wholesec = int $angle;
  $secfract = $angle - $wholesec;

  my $sig = 0;
  $sig = $_[1] if (defined $_[1]);
  $wholesec %= 60;
  $min = ($angle-$wholesec - $secfract)/60.0;
  $secfract = int ($secfract * 10**$sig + 0.5); # Add 0.5 to ensure rounding

  # Check we have not rounded too far
  if ($secfract >= 10**$sig) {
    $secfract -= 10**$sig;
    $wholesec++;
    if ($wholesec >= 60.0) {
      $wholesec -= 60;
      $min++;
      if ($min >= 60.0) {
	$min -= 60;
	$wholeangle++;
      }
    }
  }

  my $angleform = '%02';

  my ($sep1, $sep2, $sep3);
  if ($strsep eq 'HMS') {
    if ($mode eq 'H') {
      $sep1 = 'H';
    } else {
      $sep1 = 'D';
    }
    $sep2 = 'M';
    $sep3 = 'S';
  } elsif ($strsep eq 'hms') {
    if ($mode eq 'H') {
      $sep1 = 'h';
    } else {
      $sep1 = 'd';
    }
    $sep2 = 'm';
    $sep3 = 's';
  } elsif ($strsep eq 'deg') { # What if $mode eq 'H'??
    $sep1 = 'd';
    $sep2 = "'";
    $sep3 = '"';
  } else {
    $sep1 = $sep2 = $strsep;
    $sep3 = '';
  }

  if ($sig > 0) {
    $str = sprintf("${angleform}d$sep1%02d".
		   "$sep2%02d.%0${sig}d$sep3", 
		   $wholeangle, $min, $wholesec, $secfract);
  } else {
    $str = sprintf("${angleform}d$sep1%02d".
		   "$sep2%02d$sep3", 
		   $wholeangle, $min, $wholesec);
  }

  if ($sign == -1) {
    $str = '-'.$str;
  }
  return $str;
}


# Pinched from Astro::Time

sub cal2mjd($$$;$) {
  my($day, $month, $year, $ut) = @_;
  $ut = 0.0 if (!defined $ut);

  my ($m, $y);
  if ($month <= 2) {
    $m = int($month+9);
    $y = int($year-1);
  } else {
    $m = int($month-3);
    $y = int($year);
  }
  my $c = int($y/100);
  $y = $y-$c*100;
  my $x1 = int(146097.0*$c/4.0);
  my $x2 = int(1461.0*$y/4.0);
  my $x3 = int((153.0*$m+2.0)/5.0);
  return($x1+$x2+$x3+$day-678882+$ut);
}

# Return the nearest integer (ie round)
sub nint ($) {
  my ($x) = @_;
  ($x<0.0) ? return(ceil($x-0.5)) : return(floor($x+0.5))
}

sub mjd2cal($) {

  my $mjd = shift;

  my $ut = fmod($mjd,1.0);
  if ($ut<0.0) {
    $ut += 1.0;
    $mjd -= 1;
  }

  use integer;  # Calculations require integer division and modulation
  # Get the integral Julian Day number
  my $jd = nint($mjd + 2400001);

  # Do some rather cryptic calculations

  my $temp1 = 4*($jd+((6*(((4*$jd-17918)/146097)))/4+1)/2-37);
  my $temp2 = 10*((($temp1-237)%1461)/4)+5;

  my $year = $temp1/1461-4712;
  my $month =(($temp2/306+2)%12)+1;
  my $day = ($temp2%306)/10+1;

  return($day, $month, $year, $ut);

}