#
## Copyright (c) 1995-2020 University Corporation for Atmospheric Research
## All rights reserved
#
my $pkgdoc = <<'EOD';
#/**----------------------------------------------------------------------
# @file       CosmicParser.pm
#
# Routines for processing COSMIC GOX data
# 
# @author     Doug Hunt
# @since      09/30/2009
# @version    $URL: svn://ursa.cosmic.ucar.edu/trunk/src/BJutils/CosmicParser.pm $ $Id: CosmicParser.pm 13834 2011-05-10 15:36:40Z vanhove $
# cdaacTask   no
# -----------------------------------------------------------------------*/
EOD

package CosmicParser;
use base qw(BJParser); # inherit from BJParser.pm

require 5.010_000;  # Needs features from perl version 5.10

use strict;
use warnings;
use BJTools;
use TimeClass;
use vars qw($tmpdir);
use Config;
use PDL;
use PDL::NiceSlice;
use Carp;
use BJConstants;


#/**----------------------------------------------------------------------
# @sub       parse_OBSDqfit
#
# Subroutine to parse the main observable record:  OBSDqfit
# (for 'observation, quadratic fit')
# Input is the packet, with 'OBSDqfit' removed and the CRC removed
# (this routine is only called if the CRC is good).
#
# Output is a set of ASCII files in $$tmpdir (global) with low rate
# data records per antenna.  Also output is high rate data, per antenna,
# in binary files ready for conversion into opnGps format.
#
# @parameter  $self -- BJParser object
# @           $pkt -- Contents of OBSDqfit packet (without 'OBSDqfit')
# @                   as a perl string
# @return     $$tmpdir/lorate.txtXX, where XX is the antenna ID
# @           $$tmpdir/hirate.binXX, where XX is the antenna ID
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub parse_OBSDqfit {

  my $self = shift;
  my $pkt  = shift;

  my $bjdlength = length($pkt) + 8;  # length should include library ID and packet ID ('OBSDqfit')

  my $time = BJParser::parseStrip ("L", \$pkt); # 4

  my $date = defined($self->{BYDATE}) ? '.'.TimeClass->new->set_gps($time)->get_yrdoy_gps : '';
  my $tmpdir = $self->{TMPDIR};
  my $hrtype = defined($self->{HRTYPE}) ? $self->{HRTYPE} : 1; # 2 = closed loop, 1 = open loop or 3 = L2C open loop

  my ($prn, $antin, $obstype, $sampint,
      $CA_chan, $CA_snr, $CA_faz, $CA_tau,
      $P2_chan, $P2_snr, $P2_faz, $P2_tau,
      $P1_chan, $P1_snr, $P1_faz, $P1_tau) = BJParser::parseStrip ("CCCC CSdd CSdf Csff", \$pkt);

  # Compute status values from obstype and antin bytes
  my $flywheel  = ($obstype & 0x20) ? 1 : 0;
  my $l2c       = ($obstype & 0x08) ? 2 : 0;
  my $trkstatus = ($antin & 0b1111_1100) >> 2;
  my $flystatus = $flywheel | $l2c; # High rate status byte
     $antin     = ($antin & 0b0000_0011); # antin: 0b 0000 00XX, X contains antenna bits

  # Convert the pseudoranges to meters
  my $Pcodechip = $CLIGHT/10230000;  # P1 wavelength ~ 30m
  $CA_tau *= (10*$Pcodechip);
  $P1_tau *= $Pcodechip;
  $P2_tau *= $Pcodechip;

  # Observable data format (all on one line):
  # time (sec past J2000), OBSD, prn (prnxx), antenna index (hex),
  # obstype (hex--does not have stable definition--I currently ignore this one),
  # sample interval (seconds, hex), BlackJack packet size (bytes, decimal),
  # CA:, CA channel (hex), CA SNR, CA phase (cyc), CA pseudorange (meters),
  # P2:, (similar fields), P1:, (similar fields).
  my $outrec_lr = sprintf ("%ld OBSD prn%02d %02x %02x %02x %d CA: %02x %d %lf %lf  P2: %02x %d %f %f  P1: %02x %d %f %f %2d %1d\n",
                           $time, $prn,
                           $antin, $obstype, $sampint, $bjdlength,
                           $CA_chan, $CA_snr, $CA_faz , $CA_tau,
                           $P2_chan, $P2_snr, ($CA_faz - $P2_faz)*$L1_L2, $CA_tau - $P2_tau,
                           $P1_chan, $P1_snr, $CA_faz - $P1_faz, $CA_tau - $P1_tau, $trkstatus, $flywheel);

  if ($tmpdir eq 'memory') {
    $self->{OUTFILES}{"lorate.txt"} .= $outrec_lr;
  } else {
    open my $fh, '>>', "$$tmpdir/lorate.txt" or croak "Cannot open lorate.txt for appending";
    print {$fh} $outrec_lr;
    close $fh;
  }

  # As of Build 4.3, a scintillation parameter is added as a negative P1_snr
  if ($P1_snr < 0) {
    my $outrec_scin = sprintf("%ld SCIN %02d %02x %d %d\n", $time, $prn, $antin, $CA_snr, abs($P1_snr));
    my $file = "scin.txt";
    if ($tmpdir eq 'memory') {
      $self->{OUTFILES}{$file} .= $outrec_scin;
    } else {
      open my $fh, '>>', "$$tmpdir/$file" or croak "Cannot open $file for appending";
      print {$fh} $outrec_scin;
      close $fh;
    }
  }

  # Return if no high rate data present
  return if (length($pkt) == 0);

  #
  ## Now handle high rate data!
  #
  my ($CA_dop, $CA_ddop) = BJParser::parseStrip ("ff", \$pkt); # 8

  my %res = ();

  # Rtyp:  For phase: CA 0xfc, P1 0xf1, P2 0xf2.  For amp: CA 0xac, P1 0xa1, P2 0xa2.
  my %dtype = (0xfc, 'CA_faz', 0xf1, 'P1_faz', 0xf2, 'P2_faz', 0xac, 'CA_amp', 0xa1, 'P1_amp', 0xa2, 'P2_amp');

  while (length($pkt) > 0) {

    my ($Rtyp, $Rscale, $Rrate) = BJParser::parseStrip ("Ccc", \$pkt); # 3
    my $dtyp = $dtype{$Rtyp} || do { printf "Unknown residual type: $Rtyp at $time PRN$prn ANT$antin TRK$trkstatus\n"; return $time };
    if ($Rrate < 47 || $Rrate > 55) {
      printf "Rrate is $Rrate \(should be near 50\).  Mal-formed high-rate data of type $dtyp for PRN$prn at $time\n";
      return $time
    }

    # Read residuals into a hash of signed short PDLs which is indexed by 'Rtyp'
    $res{$dtyp} = zeroes(short,$Rrate);
    $res{$dtyp}->make_physical;
    ${$res{$dtyp}->get_dataref} = substr ($pkt, 0, $Rrate*2, '');
    $res{$dtyp}->upd_data;
    $res{$dtyp}->bswap2 if (!$BIGENDIAN);
  }
  return if (keys %res == 0); # no high rate resid

  # printf "%10d TRK%02d PRN%02d %s\n", $time, $trkstatus, $prn, join (' ', sort keys %res);

  #
  ## Compute high rate data vectors
  ##
  ## The number of high rate values this epoch (Rrate) is tricky to determine because
  ## under various circumstances extra values are stored at the end of various high rate
  ## vectors.  If it is defined, then CA_amp always contains the correct number of values.
  ## If it is not defined (in the case of clock link data) then CA_faz should contain
  ## the correct number
  #
  my $Rrate = defined($res{'CA_amp'}) ? $res{'CA_amp'}->nelem
            :                           $res{'CA_faz'}->nelem;
  # printf "Rate = %2d PRN%02d $time\n", $Rrate, $prn;

  my ($dt, $ttag_frac);

  #
  ## If two extra residuals are stored at the end of CA_faz, these are
  ## 'time tag fractions' used in computing high rate data times for open
  ## loop data.  These values are not used for clock link data (only open loop data)
  ## so the above definition of Rrate holds.
  #
  if (!exists($res{'CA_faz'})) {
    print "High rate CA_faz data does not exist for trkstatus = $trkstatus!  Skipping...\n";
    return $time;
  }
  if ($res{'CA_faz'}->nelem == $Rrate + 2) {
    my $ttag_frac1 = ushort($res{'CA_faz'}->slice("-2"))->sclr;
    my $ttag_frac2 = ushort($res{'CA_faz'}->slice("-1"))->sclr;
    $res{'CA_faz'}->inplace->reshape($Rrate); # chop off last two elements
    $ttag_frac  = ( ( 65536.0 * $ttag_frac1 ) + $ttag_frac2 ) / 19328000.0;

    #
    ## Compute the fractional high rate times for this seconds worth of data.
    ## For open loop data, this process involves the 'ttag_frac' encoded at the end
    ## of the CA_faz residuals.  For closed loop data, this is just even 50Hz time steps.
    #
    my $ttag_frac_last = $self->{TTAG_FRAC_LAST}[$prn] // $ttag_frac;
    $ttag_frac_last    = $ttag_frac if (!defined($self->{TIME_LAST_FRAC}[$prn]) ||
                                        $self->{TIME_LAST_FRAC}[$prn] != $time - 1); # reset if gap
    $self->{TTAG_FRAC_LAST}[$prn] = $ttag_frac;
    $self->{TIME_LAST_FRAC}[$prn] = $time;
    $dt = (-1 + $ttag_frac_last + (sequence($Rrate)+1) * (1 + $ttag_frac - $ttag_frac_last) / $Rrate);

  } else {
    $dt = (sequence($Rrate)/$Rrate) - 0.5; # ((double)i)/((double)Rrate) - 0.5;
  }


  #
  ## New for build 4.4:
  ## Determine if there is an extra L2 phase scaling factor on the end of the L2 phase
  ## residuals array.  This should be the case for build 4.4 and later data.
  ## If so, extract it and use it to uncompress the L2 phase values.  If not,
  ## use the old approach of uncompressing the values.
  #
  my $p2slope = 0;
  my $have_p2slope = ($res{'P2_faz'}->nelem == $Rrate + 1);
  if ($have_p2slope) {
    $p2slope = ( $res{'P2_faz'}->slice("-1")->sclr / 1000.0 ) * $L2_L1; # convert to units of L1 wavelengths (it will be converted back below!)
    $res{'P2_faz'}->inplace->reshape($Rrate); # chop off last element
  }

  # Attempt to correct 'wrapping'--bad compression when there is not enough dynamic range for residuals
  my ($ca_resid, $p2_resid);
  eval {
    $ca_resid = PDL::fixwrap($dt, $res{'CA_faz'});
    $p2_resid = PDL::fixwrap($dt, $res{'P2_faz'});
  };
  if ($@) {
    print "Bad return from PDL::fixwrap: $@\n";
    return $time;
  }

  my $ones = ones($Rrate); # do this once, save time

  # Variables to put into opnGps records.  Computed differently
  # depending upon track status word ($trkstatus).
  my ($ca, $p2, $casnr, $p2snr, $dfaz, $dfaz2, $camdl, $p2mdl);

  if ($trkstatus >= 21 && $trkstatus <= 26) { # clock/reference link data

    # High rate values for CA faz and P2 faz
    $ca = (( $CA_faz           +             $CA_dop *$dt + $CA_ddop*$dt**2/2)           + $ca_resid/8192) * $LAMB1;
    $p2 = ((($CA_faz - $P2_faz + ($p2slope + $CA_dop)*$dt + $CA_ddop*$dt**2/2) * $L1_L2) + $p2_resid/8192) * $LAMB2;

    # High rate CA SNR

    # New for build 4.5.3--High rate CA SNR on reference link data can be commanded on
    if (exists($res{'CA_amp'})) {
      $casnr = $SNRCONV50 * double($res{'CA_amp'}) * 10;
    } else {
      $casnr = $CA_snr * 10 * $ones; # No high rate CA SNRs for ref link data
    }
    $p2snr = $P2_snr * 10 * $ones;

    # Error values for open loop fields
    $dfaz  = $camdl = $ones * -999;
    $dfaz2 = $p2mdl = $ones * -999 if ($hrtype == 3);

  } elsif ($trkstatus == 33 || $trkstatus == 36) { # closed loop occultation data 33 = setting, 36 = rising

    if (!exists($res{'CA_amp'}) || !exists($res{'P1_amp'})) {
      print "Missing required high rate data for trkstatus = $trkstatus, skipping!\n";
      return $time;
    }

    # Not sure what I and Q mean for closed loop data...
    my $I_L1 = double($res{'CA_amp'});
    my $Q_L1 = double($res{'P1_amp'});

    # For L2C, there exists high rate SNRs for L2 even in closed loop mode
    # D. Hunt 15 Apr, 2011
    if (exists($res{'P2_amp'}) && exists($res{'P1_faz'})) {
      my $I_L2 = double($res{'P2_amp'});
      my $Q_L2 = double($res{'P1_faz'});
      $p2snr = $SNRCONV50 * sqrt($I_L2**2 + $Q_L2**2) * 10;
    } else {
      $p2snr = $P2_snr * 10 * $ones;
    }

    # High rate values for CA faz and P2 faz
    $ca = (( $CA_faz           +             $CA_dop *$dt + $CA_ddop*$dt**2/2)           + $ca_resid/8192) * $LAMB1;
    $p2 = ((($CA_faz - $P2_faz + ($p2slope + $CA_dop)*$dt + $CA_ddop*$dt**2/2) * $L1_L2) + $p2_resid/8192) * $LAMB2;

    # High rate CA SNR
    $casnr = $SNRCONV50 * sqrt($I_L1**2 + $Q_L1**2) * 10;

    # Error values for open loop fields
    $dfaz  = $camdl = $ones * -999;
    $dfaz2 = $p2mdl = $ones * -999 if ($hrtype == 3);

  } elsif ($l2c && ($trkstatus == 34 || $trkstatus == 35)) { # L2C open loop: 34 = setting, 35 = rising

    if (!exists($res{'CA_amp'}) || !exists($res{'P1_amp'}) || !exists($res{'P2_amp'}) || !exists($res{'P1_faz'})) {
      print "Missing required high rate data for trkstatus = $trkstatus, skipping!\n";
      return $time;
    }

    my $I_L1 = double($res{'CA_amp'});
    my $Q_L1 = double($res{'P1_amp'});
    my $I_L2 = double($res{'P2_amp'});
    my $Q_L2 = double($res{'P1_faz'});

    # For L2C, CA_faz = CA_model, CA_amp = I_L1
    #          P1_faz = Q_L2,     P1_amp = Q_L1
    #          P2_faz = P2_model, P2_amp = I_L2

    $dfaz  = atan2(-1*$Q_L1,$I_L1) * $RAD2CYC * $LAMB1;
    $dfaz2 = atan2(-1*$Q_L2,$I_L2) * $RAD2CYC * $LAMB2;
    $camdl =  (($CA_faz           +             $CA_dop *$dt + $CA_ddop*$dt**2/2)           + $ca_resid/8192) * $LAMB1;
    $p2mdl = ((($CA_faz - $P2_faz + ($p2slope + $CA_dop)*$dt + $CA_ddop*$dt**2/2) * $L1_L2) + $p2_resid/8192) * $LAMB2;

    # closed loop variables
    $ca    = $dfaz  + $camdl;
    $p2    = $dfaz2 + $p2mdl;

    # High rate CA SNR, P2 SNR
    $casnr = $SNRCONV50 * sqrt($I_L1**2 + $Q_L1**2) * 10;
    $p2snr = $SNRCONV50 * sqrt($I_L2**2 + $Q_L2**2) * 10;

  } elsif ($trkstatus == 34 || $trkstatus == 35) { # non-L2C open loop: 34 = setting, 35 = rising

    if (!exists($res{'CA_amp'}) || !exists($res{'P1_amp'})) {
      print "Missing required high rate data for trkstatus = $trkstatus, skipping!\n";
      return $time;
    }
    my $I_L1 = double($res{'CA_amp'});
    my $Q_L1 = double($res{'P1_amp'});

    $dfaz  = atan2(-1*$Q_L1,$I_L1) * $RAD2CYC * $LAMB1;
    $camdl =  (($CA_faz           +             $CA_dop *$dt + $CA_ddop*$dt**2/2)           + $ca_resid/8192) * $LAMB1;
    $p2    = ((($CA_faz - $P2_faz + ($p2slope + $CA_dop)*$dt + $CA_ddop*$dt**2/2) * $L1_L2) + $p2_resid/8192) * $LAMB2;
    $ca    = $dfaz + $camdl;

    # High rate CA SNR, one second P2 SNR
    $casnr = $SNRCONV50 * sqrt($I_L1**2 + $Q_L1**2) * 10;
    $p2snr = $P2_snr * 10 * $ones;

    $dfaz2 = $p2mdl = $ones * -999 if ($hrtype == 3);

  } else {
    return $time; # don't know how to handle this packet
  }

  #
  ## Reformat these high rate values into an opnGps record, either the normal open loop type (HRTYPE == 1)
  ## or the L2C open loop type (HRTYPE == 3)
  #

  my ($hrhdr, $hrrec);
  my $reserved = 0;

  if ($hrtype == 1) { # normal open loop data
    $hrhdr = pack ("SSIdCCCC", # 20 byte header
                   $Rrate, $prn, $time, $CA_tau, $reserved, $antin, $trkstatus, $flystatus);
    my $hrvals = cat ($dt, $ca, $p2, $casnr, $p2snr, $camdl, $dfaz);
    $hrrec  = pack ("fddSSdf" x $Rrate, $hrvals->xchg(0,1)->flat->list); # 36 byte records
  } elsif ($hrtype == 3) { # L2C open loop data
    $hrhdr = pack ("SSIdCCCCd", # 28 byte header
                   $Rrate, $prn, $time, $CA_tau, $reserved, $antin, $trkstatus, $flystatus, $CA_tau - $P2_tau);
    my $hrvals = cat ($dt, $ca, $p2, $casnr, $p2snr, $camdl, $dfaz, $p2mdl, $dfaz2);
    $hrrec  = pack ("fddSSdfdf" x $Rrate, $hrvals->xchg(0,1)->flat->list); # 48 byte records
  } else {
    die "HRTYPE of $hrtype is not supported!";
  }

  # Append record onto data record
  my $file = sprintf ("hirate.bin%02d", $antin);
  if ($tmpdir eq 'memory') {
    $self->{OUTFILES}{$file} .= $hrhdr;
    $self->{OUTFILES}{$file} .= $hrrec;
  } else {
    open my $fh, '>>', "$$tmpdir/$file" or croak "Cannot open $file for appending";
    print  {$fh} $hrhdr;
    print  {$fh} $hrrec;
    close   $fh;
  }

  return $time;

}

1;