#!/bin/env perl
#
##  Copyright (c) 1995-2014 University Corporation for Atmospheric Research
## All rights reserved
#
my $pkgdoc = <<'EOD';
#/**----------------------------------------------------------------------    
# @file       dump_opnGns_external.pl
# ASCII dump of one opnGns file (high rate GNSS data file)
# 
# @author     Doug Hunt
# @since      12/03/2014
# @version    $URL$ $Id$
# @cdaacTask
# @see        GNSS_High_Rate_Binary_Format_v2.2.docx
# @usage      dump_opnGns_external.pl opnGns_file [--lowrate]
# @           Dump a version 1 or version 2 opnGns file to ASCII.
# @           If --lowrate is specified, only dump low rate portion of data
# @           If --trailer is specified, then dump the trailer (no other data from the body)
# -----------------------------------------------------------------------*/
EOD

$|++;

use lib qw (. /ops/tools/lib);
use strict;
use warnings;

use PDL;
use Getopt::Long;
use Data::Dumper;
use vars qw($lowrate); # global

if (@ARGV < 1) {
  print $pkgdoc;
  exit -1;
}

my $infile  = shift;

# Handle command-line options
$lowrate  = 0;
my $trailer = 0;
GetOptions (
            "lowrate!"  => \$lowrate,  # show only low rate data
            "trailer!"  => \$trailer,
	    ) || die "Cannot parse command line options!";

# Time stamp on input file
my $inpath = $infile;

open my $ifh, $inpath or die "Cannot read $inpath";

# Read end block with version and GNSS list
seek ($ifh, -9, 2);
read ($ifh, my $buf = '', 9);

my ($version, @conids) = unpack ('CAAAAAAAA', $buf);

if ($trailer) {
  dump_trailer ($ifh, $version, \@conids);
  close ($ifh);
  exit;
}

if ($version == 1) {
  read_v1 ($ifh, \@conids);
} elsif ($version >= 2) {
  read_v2 ($ifh, \@conids, $version);
} else {
  die "Can only read version 1 or >= 2 opnGns files, found version = $version";
}

close ($ifh);

exit;

#/**----------------------------------------------------------------------
# @sub       read_v1
#
# Read in and print out a version 1 opnGns file.
#
# @parameter  $ifh -- Input file handle
# @           $conids -- Reference to list of constellation IDs
# @return     Nothing.  But prints contents of file
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub read_v1 {

  my $ifh = shift;
  my @conids = @{(shift)};

  @conids = grep /[A-Z]/, @conids;  # get rid of null entries

  my $desc_size  = 340;  # 336 bytes of descriptive trailer for each conid
  my $lr_hdr_len = 12;   # 12 bytes of fixed-length information for each low rate record

  # Seek to the beginning of the descriptive trailer section
  seek ($ifh, (-1 * $desc_size * scalar(@conids)) - 9, 2);

  #
  ## Read in trailer information on what data are stored by each constellation.
  ## Store this information in the %desc and %hdr hashes.
  #
  my %desc;  # data descriptions per conid
  my %hdr;   # header per conid
  foreach my $conid (@conids) {
    read ($ifh, my $rawdesc = '', $desc_size);

    # Extract low- and high-rate data descriptors.  Get rid of blank records
    my @lrdesc = grep { $_ } unpack ('A2'x16, substr ($rawdesc, 0,  32, ''));
    my @hrdesc = grep { $_ } unpack ('A3'x16, substr ($rawdesc, 0,  48, ''));

    $desc{$conid}{'LOWRATE'}  = [@lrdesc];
    $desc{$conid}{'HIGHRATE'} = [@hrdesc];
    $desc{$conid}{HROFFSET}   = (unpack 'l>', substr ($rawdesc, 0, 4, '')) * 1e-9; # seconds
    $desc{$conid}{'OFFSETS'}  = [unpack ('q>' x32, substr ($rawdesc, 0, 256, ''))];

    # Extract the high rate data packing format
    @{$desc{$conid}{HRPACK}} = map { substr ($_,2,1) } @hrdesc; # the perl 'pack' descriptors ('d', 'S', etc)

    # Compute the length of the high rate data section.  This is done by packing a zero
    # into each type found in HRPACK.
    $desc{$conid}{HRLEN} = pdl(map { length(pack($_, 0)) } @{$desc{$conid}{HRPACK}})->sum;

    # Compute header strings which describe the low-rate and high-rate data using RINEX3 data types
    @{$hdr{$conid}{LRHEAD}} = map { decode_desc ($_) } @lrdesc;
    @{$hdr{$conid}{HRHEAD}} = map { decode_desc ($_) } @hrdesc;

  }

  my $first = 1;

  seek ($ifh, 0, 0); # seek to beginning of file

  # Read through entire file
  while (1) {

    # Read in fixed-length portion of the epoch data record.
    my $bytes_read = read ($ifh, my $hdr = '', $lr_hdr_len);
    last if ($bytes_read != $lr_hdr_len);
    my ($tpkt, $status, $antid, $conid, $svn, $prn, $rate) = unpack ('(LCCASCS)>', $hdr);
    last if (!exists($desc{$conid})); # end of file if no constellation ID found

    # Read one double per data descriptor for this constellation ID.
    # Decode and print these values, along with the time stamp and flags.
    my $nlow = scalar (@{$hdr{$conid}{LRHEAD}});
    read ($ifh, my $lrvals = '', scalar (@{$desc{$conid}{'LOWRATE'}}) * 8);
    my @lrvals = unpack ("d>" x $nlow, $lrvals);

    # Print out low rate data with header.  Only print out header the first time
    # if --lowrate is specified
    if ($lowrate && $first) {
      printf ("GPS_sec stat ant svn prn rate"."%20s"x$nlow."\n", @{$hdr{$conid}{LRHEAD}});
      $first = 0;
    } elsif (!$lowrate) {
      printf ("GPS_sec stat ant svn prn rate"."%20s"x$nlow."\n", @{$hdr{$conid}{LRHEAD}});
    }
    printf "%10d %3d %2d %3d %1s%02d %3d"."%20.10f "x$nlow."\n",
      $tpkt, $status, $antid, $svn, $conid, $prn, $rate, @lrvals;

    # Read in high rate data
    my $hrlen = 3 + $desc{$conid}{HRLEN}; # 3 byte time offset added
    $bytes_read = read ($ifh, my $hr = '', $rate * $hrlen);
    last if ($bytes_read < $rate * $hrlen); # end of file if less high rate data than expected

    next if ($lowrate); # skip high rate printout if requested.

    # Print out high rate data
    my $nhigh = scalar(@{$desc{$conid}{HRPACK}});
    printf "  Time frac "."%20s"x$nhigh."\n", @{$hdr{$conid}{HRHEAD}};
    for (my $i=0;$i<$rate;$i++) {
      my $packspec = join '', @{$desc{$conid}{HRPACK}};
      my ($tmost, $tmid, $tmin, @hrvals) = unpack ("(CCC$packspec)>",
                                                 substr ($hr, 0, 3+$desc{$conid}{HRLEN}, ''));
      my $tf = (unpack ("N", pack ("CCCC", 0, $tmost, $tmid, $tmin)) * 1e-7) + $desc{$conid}{HROFFSET};
      printf "  %12.8f "."%20.10f "x$nhigh."\n", $tf, @hrvals;
    }

  }

}

# Place to memoize the length of the HR records, so you don't have to recompute it
# using 'pack' each time.
my %hrlen_memo;


#/**----------------------------------------------------------------------
# @sub       read_v2
#
# Read in and print out a version 2 opnGns file.
#
# @parameter  $ifh -- Input file handle
# @           $conids -- Reference to list of constellation IDs
# @return     Nothing.  But prints contents of file
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub read_v2 {

  my $ifh = shift;
  my @conids = @{(shift)};
  my $version = shift;

  @conids = grep /[A-Z]/, @conids;  # get rid of null entries

  my $lr_hdr_len = 14;   # 12 bytes of fixed-length information for each low rate record
  my $desc_size  = $version eq 2   ? 260 :
                   $version eq 255 ? 516 :
                   $version eq 3   ? 516 :  # 260 bytes of descriptive trailer for each conid for v2
                   -1;                      # 516 for v3

  die "Unknown version $version" if ($desc_size < 0);

  # Seek to the beginning of the descriptive trailer section
  seek ($ifh, (-1 * $desc_size * scalar(@conids)) - 9, 2);

  #
  ## Read in trailer information on what data are stored by each constellation.
  ## Store this information in the %desc hashe.
  #
  my %desc;  # data descriptions per conid
  foreach my $conid (@conids) {
    read ($ifh, my $rawdesc = '', $desc_size);
    $desc{$conid}{HROFFSET}   = (unpack 'l>', substr ($rawdesc, 0, 4)) * 1e-9; # seconds
  }

  my $first = 1;

  seek ($ifh, 0, 0); # seek to beginning of file

  # GPS, GLONASS, Galileo, BeiDou, SBAS considered valid
  my %valid_conid = (G => 1, R => 1, E => 1, C => 1, S => 1, J => 1);

  # Read through entire file
  while (1) {

    # Read in fixed-length portion of the epoch data record.
    my $bytes_read = read ($ifh, my $hdr = '', $lr_hdr_len);
    last if ($bytes_read != $lr_hdr_len);
    my ($tpkt, $status, $antid, $conid, $svn, $prn, $rate, $nlr, $nhr) = unpack ('(LCCASCSCC)>', $hdr);
    last if (!$valid_conid{$conid}); # end of file if no constellation ID found

    read ($ifh, my $lrhdr = '', $nlr*2);         # Variable length LR headers
    read ($ifh, my $hrhdr = '', $nhr*3) or last; # Variable length HR headers

    # Read one double per data descriptor for this constellation ID.
    # Decode and print these values, along with the time stamp and flags.
    read ($ifh, my $lrvals = '', $nlr * 8);
    my @lrvals = unpack ("d>" x $nlr, $lrvals);

    my @lrdesc = map { decode_desc ($_) } unpack ('(a2)*', $lrhdr);

    # Print out low rate data with header.  Only print out header the first time
    # if --lowrate is specified
    if ($lowrate && $first) {
      printf ("GPS_sec stat ant svn prn rate"."%20s"x$nlr."\n", @lrdesc);
      $first = 0;
    } elsif (!$lowrate) {
      printf ("GPS_sec stat ant svn prn rate"."%20s"x$nlr."\n", @lrdesc);
    }
    printf "%10d %3d %2d %3d %1s%02d %3d"."%20.10f "x$nlr."\n",
      $tpkt, $status, $antid, $svn, $conid, $prn, $rate, @lrvals;

    #
    ## Read in high rate data
    #

    # Compute length of high rate data section from high rate header
    my @pack_desc = map { substr ($_, 2, 1) } unpack("(A3)*", $hrhdr); # the perl 'pack' descriptors ('d', 'S', etc)
    my $packspec  = join '', @pack_desc;

    # Use a memo hash to look up high rate data length instead of computing it each time.
    my $hrlen = $hrlen_memo{$packspec} //= 3 + pdl(map { length(pack($_, 0)) } @pack_desc)->sum;  # 3 byte time offset + sum of high rate data lengths

    $bytes_read = read ($ifh, my $hr = '', $rate * $hrlen);
    last if ($bytes_read < $rate * $hrlen); # end of file if less high rate data than expected

    next if ($lowrate); # skip high rate printout if requested.

    # Print out high rate data
    my @hrdesc = map { decode_desc ($_) } unpack ('(A3)*', $hrhdr);


    printf "  Time frac "."%20s"x$nhr."\n", @hrdesc;
    for (my $i=0;$i<$rate;$i++) {
      my ($tmost, $tmid, $tmin, @hrvals) = unpack ("(CCC$packspec)>", substr ($hr, 0, $hrlen, ''));
      my $tf = (unpack ("N", pack ("CCCC", 0, $tmost, $tmid, $tmin)) * 1e-7) + $desc{$conid}{HROFFSET};
      printf "  %12.8f "."%20.10f "x$nhr."\n", $tf, @hrvals;
    }

  }

}


#/**----------------------------------------------------------------------
# @sub       decode_desc
#
# Given a 16 bit data descriptor, return the RINEX 3 code
# See Section 2 of the GNSS_High_Rate_Binary_Format document
#
# @parameter  $desc -- Two byte data descriptor.  Can have a pack byte at the end.
# @return     $r3   -- RINEX 3 code 'L1P', etc.
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub decode_desc {
  my $desc = shift;

  my @ot = (qw(C L D S X X X X)); # table for decoding RINEX 3 observation type
  my ($obstype, $freq_band, $trk_attr, $meta_data) =
    map { bin2dec($_) } unpack ("a3 a3 a5 a5", unpack ("B16", $desc));
  my $r3 = sprintf "%s%s%s", $ot[$obstype], $freq_band+1, chr(ord('A') + $trk_attr-1);

  # Add the letter code for the meta-data flag, if present
  my $flag = $meta_data ? chr(ord('A') + $meta_data-1) : '';
  $r3 .= " ($flag)" if ($flag);

  return $r3;
}

#/**----------------------------------------------------------------------
# @sub       dump_trailer
#
# Read in and print out the trailer any version file
#
# @parameter  $ifh -- Input file handle
# @           $version -- opnGns version number
# @           $conids -- ref of list of constellation IDs
# @return     Nothing.  But prints contents of trailer
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub dump_trailer {

  my $ifh = shift;
  my $version = shift;
  my @conids = @{(shift)};

  @conids = grep /[A-Z]/, @conids;  # get rid of null entries

  my $desc_size  = $version == 1   ? 340 :
                   $version == 2   ? 260 : # 260 bytes of descriptive trailer for each conid for v2
                   $version == 3   ? 516 : # 516 for v3
                   $version eq 255 ? 516 :
                                      -1;

  die "Unknown version $version" if ($desc_size < 0);

  # Seek to the beginning of the descriptive trailer section
  seek ($ifh, (-1 * $desc_size * scalar(@conids)) - 9, 2);

  #
  ## Read in trailer information on what data are stored by each constellation.
  ## Store this information in the %desc and %hdr hashes.
  #
  my %desc;  # data descriptions per conid
  my %hdr;   # header per conid
  $desc{VERSION} = $version;
  foreach my $conid (@conids) {
    read ($ifh, my $rawdesc = '', $desc_size);

    # Extract low- and high-rate data descriptors.  Get rid of blank records
    if ($version == 1) {
      my @lrdesc = grep { $_ } unpack ('A2'x16, substr ($rawdesc, 0,  32, ''));
      my @hrdesc = grep { $_ } unpack ('A3'x16, substr ($rawdesc, 0,  48, ''));

      $desc{$conid}{'LOWRATE'}  = [@lrdesc];
      $desc{$conid}{'HIGHRATE'} = [@hrdesc];

      # Extract the high rate data packing format
      @{$desc{$conid}{HRPACK}} = map { substr ($_,2,1) } @hrdesc; # the perl 'pack' descriptors ('d', 'S', etc)

      # Compute the length of the high rate data section.  This is done by packing a zero
      # into each type found in HRPACK.
      $desc{$conid}{HRLEN} = pdl(map { length(pack($_, 0)) } @{$desc{$conid}{HRPACK}})->sum;

      # Compute header strings which describe the low-rate and high-rate data using RINEX3 data types
      @{$hdr{$conid}{LRHEAD}} = map { decode_desc ($_) } @lrdesc;
      @{$hdr{$conid}{HRHEAD}} = map { decode_desc ($_) } @hrdesc;
    }

    my $maxprn  = $version eq 2 ? 32  :  64;
    my $consize = $version eq 2 ? 256 : 512;

    $desc{$conid}{HROFFSET}   = (unpack 'l>', substr ($rawdesc, 0, 4, '')) * 1e-9; # seconds
    $desc{$conid}{'OFFSETS'}  = [unpack ('q>' x$maxprn, substr ($rawdesc, 0, $consize, ''))];

  }

  $Data::Dumper::Indent = 1;
  print Data::Dumper->Dump([\%desc], [qw(*desc)]);

}

# Convert a binary string in ascii: '01100110' to a perl scalar value
# From perlmonks:  http://www.perlmonks.org/?node_id=2664
sub bin2dec { return unpack("N", pack("B32", substr("0" x 32 . shift, -32))); }