#
## Copyright (c) 1995-2020 University Corporation for Atmospheric Research
## All rights reserved
#
my $pkgdoc = <<'EOD';
#/**----------------------------------------------------------------------
# @file       Cosmic2Parser.pm
#
# Routines for processing COSMIC-2 TRIG data, version 1 dated 7/22/2015.
# These routines are also used for the TRIG simulator.
#
# @author     Doug Hunt
# @since      10/9/2014
# @version    $URL$ $Id$
# cdaacTask   no
# -----------------------------------------------------------------------*/
EOD

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

require 5.010_000;  # Needs features from perl version 5.10

use strict;
use warnings;
use BJTools;
use TimeClass;
use PDL;
use PDL::NiceSlice;
use Carp;
use BJConstants;
use SatInfo;

# Read Bernese config file to get information on GLONASS frequency slot mapping.
# Do this once at load time.
# original --------------------------
#my $Ivs = 'I14'; # Hard-wired to IGS 14 coordinates
#my $satinfo = SatInfo->new("SATELLIT.$Ivs.bern")->parse;

# Jun modified ------Use new file provided by Yong Chen, include E svn--------------
my $Ivs = 'I20'; # Hard-wired to IGS 14 coordinates
my $satinfo = SatInfo->new("SATELLIT_$Ivs.SAT")->parse;
#
## Global data
#

#
# Descriptions and low and high ranges for each TRIG State Of Health data type.
# From JPL spreadsheet C2A_TriG_SOH.xlsx dated Nov. 18, 2016
# Red limit minimum temperatures modified based on document:  Payload_Limits_Table_Draft_TGRS-IVM-RFB_v4
# on 4/3/2018.
no warnings 'once';
%Cosmic2Parser::trig_SOH_desc = (
  'T005' => ['Temp 2 SciP','-10','70'],
  'T006' => ['Temp 1 SciP','-10','70'],
  'T014' => ['Temp 1 Sci DSP','-10','70'],
  'T021' => ['Temp DC-DC 3.3V PDC','-10','70'],
  'T022' => ['Temp DC-DC 5V PDC','-10','70'],
  'T023' => ['Temp DC-DC 15V PDC','-10','70'],
  'T024' => ['Temp OCXO CDC','-10','70'],
  'T025' => ['Temp Regulator CDC','-10','70'],
  'T026' => ['Temp CDCM7005 CDC','-10','70'],
  'T027' => ['Temp TCXO CDC','-10','70'],
  'T032' => ['Temp1_1 RFDC1','-10','70'],
  'T033' => ['Temp1_2 RFDC1','-10','70'],
  'T034' => ['Temp1_3 RFDC1','-10','70'],
  'T035' => ['Temp1_4 RFDC1','-10','70'],
  'T036' => ['Temp1_5 RFDC1','-10','70'],
  'T037' => ['Temp1_6 RFDC1','-10','70'],
  'T038' => ['Temp1_7 RFDC1','-10','70'],
  'T039' => ['Temp1_8 RFDC1','-10','70'],
  'T040' => ['Temp1_9 RFDC1','-10','70'],
  'T041' => ['Temp1_10 RFDC1','-10','70'],
  'T042' => ['Temp1_11 RFDC1','-10','70'],
  'T043' => ['Temp1_12 RFDC1','-10','70'],
  'T044' => ['Temp1_13 RFDC1','-10','70'],
  'T045' => ['Temp1_14 RFDC1','-10','70'],
  'T046' => ['Temp1_15 RFDC1','-10','70'],
  'T047' => ['Temp1_16 RFDC1','-10','70'],
  'T048' => ['Temp2_1 RFDC2','-10','70'],
  'T049' => ['Temp2_2 RFDC2','-10','70'],
  'T050' => ['Temp2_3 RFDC2','-10','70'],
  'T051' => ['Temp2_4 RFDC2','-10','70'],
  'T052' => ['Temp2_5 RFDC2','-10','70'],
  'T053' => ['Temp2_6 RFDC2','-10','70'],
  'T054' => ['Temp2_7 RFDC2','-10','70'],
  'T055' => ['Temp2_8 RFDC2','-10','70'],
  'T056' => ['Temp2_9 RFDC2','-10','70'],
  'T057' => ['Temp2_10 RFDC2','-10','70'],
  'T058' => ['Temp2_11 RFDC2','-10','70'],
  'T059' => ['Temp2_12 RFDC2','-10','70'],
  'T060' => ['Temp2_13 RFDC2','-10','70'],
  'T061' => ['Temp2_14 RFDC2','-10','70'],
  'T062' => ['Temp2_15 RFDC2','-10','70'],
  'T063' => ['Temp2_16 RFDC2','-10','70'],
  'T101' => ['Temp 2 NAVP','-10','70'],
  'T102' => ['Temp 1 NAVP','-10','70'],
  'U000' => ['Description Board','0','0'],
  'V000' => ['SOH_ADC_VA SciP','2.97','3.63'],
  'V001' => ['3.3V_SW SciP','2.97','3.63'],
  'V002' => ['1.35V SciP','1.215','1.485'],
  'V003' => ['1.5V SciP','1.35','1.65'],
  'V004' => ['PPC750FX_IO SciP','1.62','1.98'],
  'V007' => ['DGND SciP','-0.1','0.1'],
  'V008' => ['SOH_ADC_VA Sci DSP','2.97','3.63'],
  'V009' => ['2.5V Sci DSP','2.25','2.75'],
  'V010' => ['1.5V Sci DSP','1.35','1.65'],
  'V011' => ['1.2V_A Sci DSP','1.08','1.32'],
  'V015' => ['DGND Sci DSP','-0.1','0.1'],
  'V016' => ['10V PDC',9,11],
  'V017' => ['5V PDC','4.5','5.5'],
  'V018' => ['3.3V PDC','2.97','3.63'],
  'V019' => ['4V_RFDC0-1 PDC','3.6','4.4'],
  'V020' => ['4V_RFDC2-3 PDC','3.6','4.4'],
  'V096' => ['SOH_ADC_VA NAVP','2.97','3.63'],
  'V097' => ['3.3V_SW NAVP','2.97','3.63'],
  'V098' => ['1.2V NAVP','1.08','1.32'],
  'V099' => ['1.5V NAVP','1.35','1.65'],
  'V100' => ['2.5V NAVP','2.25','2.75'],
  'V103' => ['DGND NAVP','-0.1','0.1']
);
my $nan = PDL->new('nan')->sclr; # Not a Number
use warnings 'once';

# Hash of most recent times a PRN (including conid, eg 'G01') is tracked.
# Used for detecting the first measurement of a track, which can be corrupted.
# D. Hunt 9/10/2018
my %OBSDq12ab_last_time_tracked;
my %OBSDq12cd_last_time_tracked;

my $bf_factor = sqrt(3); # When beam forming, multiply noiseIQ by this factor

#/**----------------------------------------------------------------------
# @sub       parse_OBSDq12a
#
# Subroutine to parse the new C-2 observable record:  OBSDq12a
# Input is the packet, with 'OBSDq12a' 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.
#
# NOTE:  This subroutine also parses the OBSDq12b packet, which is just
# like OBSDq12a, but with extra delay error values at the end, which
# are not used in normal processing.
#
# @parameter  $self -- BJParser object
# @           $pkt -- Contents of OBSDqfit packet (without 'OBSDqfit')
# @                   as a perl string
# @return     $self->{OUTFILES}{lorate} updated
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub parse_OBSDq12a {

  my $self = shift;
  my $pkt  = shift;  # Should include entire packet except the first 8 chars:  'OBSDq12a'

  my $time = BJParser::parseStrip ("L", \$pkt); # 4 bytes, big endian unsigned integer

  my $tmpdir = $self->{TMPDIR};
  my $hrtype = defined($self->{HRTYPE}) ? $self->{HRTYPE} : 1; # 2 = closed loop, 1 = open loop or 3 = L2C open loop

  # C = unsigned byte, c = signed byte, S = unsigned short, d = double, f = float.  See 'perldoc pack'
  my ($conid, $prn, $ant, $fdma_chan, $obstype, $trackstat, $sample_int, $primary_chan,
      $signal_type1, $snr1, $phase1, $phaserate1, $phaseraterate1, $range1, $s4_1, $sigphi1,
      $signal_type2, $snr2, $phase2, $phaserate2, $phaseraterate2, $range2, $s4_2, $sigphi2)
      = BJParser::parseStrip ("CCCcCCCC CSddddff CSddddff", \$pkt);

  # Signal type1:  0x03 = GPS or GLONASS L1/CA
  # Signal type2:  0x30 = GPS L2P, 0x38 = GPS L2C (combined), 0x23 = GLONASS L2
  $conid = chr($conid); # convert to character
  $ant = 1 if ($ant > 3);  # !!! kludge for bad ant ids in JPL 10.22.210
  $primary_chan += 0; # canonicalize

  $prn -= 100 while ($prn > 100); # GLONASS PRNs can either be in the 100's or the 200's!

  my $fullprn = sprintf "%1s%02d", $conid, $prn;

  # print "SATTYPE$conid, T$time PRN$prn ANT$ant, FDMA$fdma_chan\n"; # DEBUG
  #
  ## Output scintillation parameters
  #
  $self->{OUTFILES}{'lorate_scin'}{$ant}{$time}{$fullprn} = [$s4_1, $sigphi1, $s4_2, $sigphi2];

  #
  ## Append GNSS records to the lorate_12a files
  #

  # For the current test data, all phases and ranges are in cycles.  For the real data, they will be
  # in time (secs, secs/sec, secs/sec**2, etc).
  die "Sample interval = $sample_int too high (> 100)" if ($sample_int > 100);

  # Note that Galileo and Beidou paths are for use only in the Day In The Life test.  These
  # constellations are not yet fully supported.
  my ($f1, $f2);
  if ($conid eq 'R') {
    #($f1, $f2) = SatInfo::glonass_freq($fdma_chan - 8); # sometimes 8 is not added to fdma_chan
    ($f1, $f2) = SatInfo::glonass_freq($fdma_chan);
  } elsif ($conid eq 'C') {
    ($f1, $f2) = ($B1, $B2);  # Beidou frequencies
  } elsif ($conid eq 'E') {
    #($f1, $f2) = ($E1, $E5b); # Galileo frequencies
    ($f1, $f2) = ($E1, $E5a); # Galileo frequencies, Jun modified
  } else {
    ($f1, $f2) = ($F1, $F2);  # GPS frequencies
  }

  my $phase1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f1; # seconds -> cycles
  my $phase2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f2;
  my $range1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
  my $range2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;

  #
  ## Sometimes the first epoch of a track is bogus.  Look for anomalous phase1resid at the
  ## first time of a packet.  If it is really odd, replace all 1st residuals with zero
  ## D. Hunt 7/10/2018
  #

  if (!defined($OBSDq12ab_last_time_tracked{$fullprn}) || ($time - $OBSDq12ab_last_time_tracked{$fullprn} > 10)) {

    # 10 sigma outlier test at beginning of track
    my ($avg, $stdv) = $phase1resid->slice("1:-1")->stats;
    if (abs($phase1resid->at(0) - $avg) > 10*$stdv) {
      print "$fullprn $time, Failed 1st residual check: $phase1resid\n";
      $phase1resid->slice('0') .= 0;
      $phase2resid->slice('0') .= 0;
      $range1resid->slice('0') .= 0;
      $range2resid->slice('0') .= 0;
    }
  }
  $OBSDq12ab_last_time_tracked{$fullprn} = $time;

  # If this is an OBSDq12b packet, these extra fields are present.  Not sure what to do with them...
  # The reading code is commented out for speed.  D. Hunt 8/4/2017
  #my ($delayError1, $delayError2);
  #if (length($pkt) > 0) {
  #  $delayError1 = pdl(BJParser::parseStrip ("s"x10, \$pkt)); # ResTau1
  #  $delayError2 = pdl(BJParser::parseStrip ("s"x10, \$pkt)); # ResTau2
  #}

  # Convert phases from seconds to cycles
  $phase1 *= $f1; # seconds * cycles/second = cycles
  $phase2 *= $f2; # seconds * cycles/second = cycles
  $phaserate1 *= $f1; # sec/sec * cycles/sec = cycles/sec
  $phaserate2 *= $f2; # sec/sec * cycles/sec = cycles/sec
  $phaseraterate1 *= $f1; # sec/sec**2 * cycles/sec = cycles/sec**2
  $phaseraterate2 *= $f2; # sec/sec**2 * cycles/sec = cycles/sec**2

  # Convert ranges from seconds to meters
  $range1 *= $CLIGHT; # seconds * meters/second = meters
  $range2 *= $CLIGHT; # seconds * meters/second = meters

  # dt = (-9, -8, ..., 0)
  my $dt   = -1 * sequence($sample_int)->slice("-1:0"); # optimize by pre-declaring???
  my $phs1 = (($phase1 + $phaserate1*$dt + $phaseraterate1*$dt**2/2) + $phase1resid); # L1 cycles
  my $phs2 = (($phase2 + $phaserate2*$dt + $phaseraterate2*$dt**2/2) + $phase2resid); # L2 cycles
  my $rng1 = ($range1 + $range1resid); # meters
  my $rng2 = ($range2 + $range2resid); # meters

  my $outrec_lr;
  my $maxrng = 50_000_000; # 50,000 km max range.
  for (my $i=0;$i<$sample_int;$i++) {
    my $t = $time + $dt->at($i);

    # Sometimes at the beginning of a track, the first residuals are garbage (unrealistically large)
    # Cull these out.
    if ($rng1->at($i) > $maxrng || $rng1->at($i) < 0 || $rng2->at($i) > $maxrng || $rng2->at($i) < 0) {
      print "Bad range found for PRN $fullprn at time = $t, skipping\n";
      next;
    }

    # Sometimes the phase values are NaNs.  Not sure why.  Skip these.  D. Hunt 2/19/2020
    if ($phs1->at($i) =~ /nan/i || $phs2->at($i) =~ /nan/i) {
      print "NaN value for phase found for PRN $fullprn at time = $t, skipping\n";
      next;
    }

    # Observable data format:
    # time (GPS seconds), prn (G|Rxx), antenna index (hex), GLONASS FDMA channel, or 0,
    # obstype (hex--does not have stable definition--I currently ignore this one),
    # tracking status (also currently undefined)
    # sample interval (seconds, hex), Primary Channel number
    # SIG1:, signal type 1 (hex), SNR, phase (cyc), pseudorange (meters),
    # SIG2:, signal type 2 (hex), SNR, phase (cyc), pseudorange (meters)
    # Note that these are 10-second SNRs.
    $self->{OUTFILES}{'lorate'}{$ant}{$t}{$fullprn}{$primary_chan} =
      {SATTYPE => $conid,
       FDMACHAN => $fdma_chan, OBSTYPE => $obstype, TRACKSTAT => $trackstat,
       SAMPLEINT => $sample_int, PRIMARYCHAN => $primary_chan,
       SIGTYPE1 => $signal_type1, SNR1 => $snr1, PHS1 => $phs1->at($i), RNG1 => $rng1->at($i),
       SIGTYPE2 => $signal_type2, SNR2 => $snr2, PHS2 => $phs2->at($i), RNG2 => $rng2->at($i) };
  }

  return $time;

}


#/**----------------------------------------------------------------------
# @sub       parse_OBSDq12c
#
# Subroutine to parse the new C-2 observable records:  OBSDq12c (GPS)
# and OBSDq12d (GLONASS).
# Input is the packet, with 'OBSDq12c/d' 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 OBSDq12c/d packet (without 'OBSDq12c/d')
# @                   as a perl string
# @return     $self->{OUTFILES}{lorate} updated
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub parse_OBSDq12c {

  my $self = shift;
  my $pkt  = shift;  # Should include entire packet except the first 8 chars:  'OBSDq12a'

  my $time = BJParser::parseStrip ("L", \$pkt); # 4 bytes, big endian unsigned integer

  my $tmpdir = $self->{TMPDIR};
  my $hrtype = defined($self->{HRTYPE}) ? $self->{HRTYPE} : 1; # 2 = closed loop, 1 = open loop or 3 = L2C open loop

  # C = unsigned byte, c = signed byte, S = unsigned short, d = double, f = float.  See 'perldoc pack'
  my ($conid, $prn, $ant, $fdma_chan, $obstype, $trackstat, $sample_int, $primary_chan,
      $signal_type1, $snr1avg, $phase1, $phaserate1, $phaseraterate1, $range1, $s4_1,
      $signal_type2, $snr2avg, $phase2, $phaserate2, $phaseraterate2, $range2, $s4_2, $azimuth, $elev, $lsa)
      = BJParser::parseStrip ("CSCcCCCC CSddddf CSddddfffs", \$pkt);

  # Signal type1:  0x03 = GPS or GLONASS L1/CA
  # Signal type2:  0x30 = GPS L2P, 0x38 = GPS L2C (combined), 0x23 = GLONASS L2
  $conid = chr($conid); # convert to character
  $ant = 1 if ($ant > 3);  # !!! kludge for bad ant ids in JPL 10.22.210
  $primary_chan += 0; # canonicalize

  $prn -= 100 while ($prn > 100); # GLONASS PRNs can either be in the 100's or the 200's!

  my $fullprn = sprintf "%1s%02d", $conid, $prn;

  # print "SATTYPE$conid, T$time PRN$prn ANT$ant, FDMA$fdma_chan\n"; # DEBUG
  #
  ## Output scintillation parameters
  #
  $self->{OUTFILES}{'lorate_scin'}{$ant}{$time}{$fullprn} = [$s4_1, $Dplib::err, $s4_2, $Dplib::err];

  #
  ## Append GNSS records to the lorate_12a files
  #

  # For the current test data, all phases and ranges are in cycles.  For the real data, they will be
  # in time (secs, secs/sec, secs/sec**2, etc).
  die "Sample interval = $sample_int too high (> 100)" if ($sample_int > 100);

  # Note that Galileo and Beidou paths are for use only in the Day In The Life test.  These
  # constellations are not yet fully supported.
  my ($f1, $f2);
  if ($conid eq 'R') {
    #($f1, $f2) = SatInfo::glonass_freq($fdma_chan - 8); # sometimes 8 is not added to fdma_chan
    ($f1, $f2) = SatInfo::glonass_freq($fdma_chan);
  } elsif ($conid eq 'C') {
    ($f1, $f2) = ($B1, $B2);  # Beidou frequencies
  } elsif ($conid eq 'E') {
    #($f1, $f2) = ($E1, $E5b); # Galileo frequencies
    ($f1, $f2) = ($E1, $E5a); # Galileo frequencies
  } else {
    ($f1, $f2) = ($F1, $F2);  # GPS frequencies
  }

  # Jun suspicious!
  my ($phase1resid, $phase2resid, $range1resid, $range2resid, $snr1, $snr2, $restau1, $restau2);
  if ($conid eq 'G') {
    $phase1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f1; # seconds -> cycles
    $phase2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f2;
    $range1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    $range2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    $snr1 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    $snr2 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
  } elsif ($conid eq 'R') {
    $phase1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f1; # seconds -> cycles
    $phase2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f2;
    $range1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    $range2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    $restau1 = pdl(BJParser::parseStrip ("S"x10, \$pkt)); # Not sure what to do with these!!!
    $restau2 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    $snr1 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    $snr2 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    #} elsif ($conid eq 'E') { # Jun added Galileo	  
    #$phase1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f1; # seconds -> cycles
    #$phase2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $f2;
    #$range1resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    #$range2resid = pdl(BJParser::parseStrip ("f"x10, \$pkt)) * $CLIGHT;
    ##$restau1 = pdl(BJParser::parseStrip ("S"x10, \$pkt)); # Not sure what to do with these!!!
    ##$restau2 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    #$snr1 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
    #$snr2 = pdl(BJParser::parseStrip ("S"x10, \$pkt));
  }

  #
  ## Sometimes the first epoch of a track is bogus.  Look for anomalous phase1resid at the
  ## first time of a packet.  If it is really odd, replace all 1st residuals with zero
  ## D. Hunt 7/10/2018
  #

  if (!defined($OBSDq12cd_last_time_tracked{$fullprn}) || ($time - $OBSDq12cd_last_time_tracked{$fullprn} > 10)) {

    print "$conid\n";
    # 10 sigma outlier test at beginning of track
    my ($avg, $stdv) = $phase1resid->slice("1:-1")->stats;
    if (abs($phase1resid->at(0) - $avg) > 10*$stdv) {
      print "$fullprn $time, Failed 1st residual check: $phase1resid\n";
      $phase1resid->slice('0') .= 0;
      $phase2resid->slice('0') .= 0;
      $range1resid->slice('0') .= 0;
      $range2resid->slice('0') .= 0;
    }
  }
  $OBSDq12cd_last_time_tracked{$fullprn} = $time;

  # Convert phases from seconds to cycles
  $phase1 *= $f1; # seconds * cycles/second = cycles
  $phase2 *= $f2; # seconds * cycles/second = cycles
  $phaserate1 *= $f1; # sec/sec * cycles/sec = cycles/sec
  $phaserate2 *= $f2; # sec/sec * cycles/sec = cycles/sec
  $phaseraterate1 *= $f1; # sec/sec**2 * cycles/sec = cycles/sec**2
  $phaseraterate2 *= $f2; # sec/sec**2 * cycles/sec = cycles/sec**2

  # Convert ranges from seconds to meters
  $range1 *= $CLIGHT; # seconds * meters/second = meters
  $range2 *= $CLIGHT; # seconds * meters/second = meters

  # dt = (-9, -8, ..., 0)
  my $dt   = -1 * sequence($sample_int)->slice("-1:0"); # optimize by pre-declaring???
  my $phs1 = (($phase1 + $phaserate1*$dt + $phaseraterate1*$dt**2/2) + $phase1resid); # L1 cycles
  my $phs2 = (($phase2 + $phaserate2*$dt + $phaseraterate2*$dt**2/2) + $phase2resid); # L2 cycles
  my $rng1 = ($range1 + $range1resid); # meters
  my $rng2 = ($range2 + $range2resid); # meters

  my $outrec_lr;
  my $maxrng = 50_000_000; # 50,000 km max range.
  for (my $i=0;$i<$sample_int;$i++) {
    my $t = $time + $dt->at($i);

    # Sometimes at the beginning of a track, the first residuals are garbage (unrealistically large)
    # Cull these out.
    if ($rng1->at($i) > $maxrng || $rng1->at($i) < 0 || $rng2->at($i) > $maxrng || $rng2->at($i) < 0) {
      print "Bad range found for PRN $fullprn at time = $t, skipping\n";
      next;
    }

    # Observable data format:
    # time (GPS seconds), prn (G|Rxx), antenna index (hex), GLONASS FDMA channel, or 0,
    # obstype (hex--does not have stable definition--I currently ignore this one),
    # tracking status (also currently undefined)
    # sample interval (seconds, hex), Primary Channel number
    # SIG1:, signal type 1 (hex), SNR, phase (cyc), pseudorange (meters),
    # SIG2:, signal type 2 (hex), SNR, phase (cyc), pseudorange (meters)
    $self->{OUTFILES}{'lorate'}{$ant}{$t}{$fullprn}{$primary_chan} =
      {SATTYPE => $conid,
       FDMACHAN => $fdma_chan, OBSTYPE => $obstype, TRACKSTAT => $trackstat,
       SAMPLEINT => $sample_int, PRIMARYCHAN => $primary_chan,
       AZIMUTH => $azimuth, ELEV => $elev, LSA => $lsa, # Note these are 10 second values repeated
       SIGTYPE1 => $signal_type1, SNR1 => $snr1->at($i), PHS1 => $phs1->at($i), RNG1 => $rng1->at($i),
       SIGTYPE2 => $signal_type2, SNR2 => $snr2->at($i), PHS2 => $phs2->at($i), RNG2 => $rng2->at($i)};
  }

  return $time;

}


# Allocate space for reading in I and Q values.  Do this once at load time.
my $I1 = zeroes(short,100);
my $I2 = zeroes(short,100);
my $Q1 = zeroes(short,100);
my $Q2 = zeroes(short,100);

# opnGns v2 Epoch Data Record low rate descriptor
my @lrdesc = ((pack 'S>', 0b000_000_00011_00000),  # RINEX3 = 'C1C', C/A pseudorange
              (pack 'S>', 0b000_000_00011_01101),  # RINEX3 = 'C1C', C/A pseudorange, with model flag ('M')
              (pack 'S>', 0b100_111_00001_00000),  # 4=reserved, freq = 8, track = A -- MSL height
              (pack 'S>', 0b100_111_00010_00000),  # 4=reserved, freq = 8, track = B -- azimuth offset of ray from vel. vector in deg.
              (pack 'S>', 0b100_111_00011_00000),  # 4=reserved, freq = 8, track = C -- latitude of ray at LSA
              (pack 'S>', 0b100_111_00100_00000),  # 4=reserved, freq = 8, track = D -- longitude of ray at LSA
              (pack 'S>', 0b100_111_00101_00000),  # 4=reserved, freq = 8, track = E -- beam form phase offset 1 (center - left), L1 milli-cyc
              (pack 'S>', 0b100_111_00110_00000),  # 4=reserved, freq = 8, track = F -- beam form phase offset 2 (center - right), L1 milli-cyc
            );

# suspicious Jun
# opnGns v2 Epoch Data Record high rate descriptor for L2P
my @hrdesc_l2p_gps = map { pack ('S>A', @$_) }
  ([0b001_000_00011_00000, 'd'], # RINEX3 = 'L1C', L1 CA phase
   [0b001_001_10000_00000, 'd'], # RINEX3 = 'L2P', L2 P phase
   [0b001_000_00011_01101, 'd'], # RINEX3 = 'L1C', L1 CA phase, 'model' flag set
   [0b011_000_00011_00000, 'S'], # RINEX3 = 'S1C', L1 C/A SNR
   [0b011_001_10000_00000, 'S'], # RINEX3 = 'S2P', L2 P SNR
 );
my @hrdesc_l2c_gps = map { pack ('S>A', @$_) }
  ([0b001_000_00011_00000, 'd'], # RINEX3 = 'L1C', L1 CA phase
   [0b001_001_01100_00000, 'd'], # RINEX3 = 'L2L' (L2C(L) bit-free pilot signal
   [0b001_000_00011_01101, 'd'], # RINEX3 = 'L1C', L1 CA phase, 'model' flag set
   [0b011_000_00011_00000, 'S'], # RINEX3 = 'S1C', L1 C/A SNR
   [0b011_001_01100_00000, 'S'], # RINEX3 = 'S2L' (L2C(L) bit-free pilot signal SNR
 );
my @hrdesc_l2c_glonass = map { pack ('S>A', @$_) }
  ([0b001_000_00011_00000, 'd'], # RINEX3 = 'L1C', L1 CA phase
   [0b001_001_00011_00000, 'd'], # RINEX3 = 'L2C' (L2C(L) bit-free pilot signal
   [0b001_000_00011_01101, 'd'], # RINEX3 = 'L1C', L1 CA phase, 'model' flag set
   [0b011_000_00011_00000, 'S'], # RINEX3 = 'S1C', L1 C/A SNR
   [0b011_001_00011_00000, 'S'], # RINEX3 = 'S2C' (L2C(L) bit-free pilot signal SNR
 );

# Jun added
my @hrdesc_l2c_galileo = map { pack ('S>A', @$_) }
  ([0b001_000_00011_00000, 'd'], # RINEX3 = 'L1C', L1 CA phase
   [0b001_100_10001_00000, 'd'], # RINEX3 = 'L5Q',
   [0b001_000_00011_01101, 'd'], # RINEX3 = 'L1C', L1 CA phase, 'model' flag set
   [0b011_000_00011_00000, 'S'], # RINEX3 = 'S1C', L1 C/A SNR
   [0b011_100_10001_00000, 'S'], # RINEX3 = 'S5Q' 
 );

#/**----------------------------------------------------------------------
# @sub       parse_OCCDq12a
#
# Subroutine to parse the new C-2 observable record:  OBSDq12a
# Input is the packet, with 'OBSDq12a' 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_OCCDq12a {

  my $self = shift;
  my $pkt  = shift; # Should include entire packet except the first 8 chars:  'OCCDq12a'

  print "Warning: OCCDq12a parsing no longer supported\n";
  # die "OCCDq12a parsing no longer supported";

}


#/**----------------------------------------------------------------------
# @sub       parse_OCCDq12b
#
# Subroutine to parse the new C-2 observable record:  OCCDq12b
# Input is the packet, with 'OCCDq12b' 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 OCCDq12b packet (without 'OCCDq12b')
# @                   as a perl string
# @return     $self->{OUTFILES}{'hirate'} updated
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub parse_OCCDq12b {

  my $self = shift;
  my $pkt  = shift; # Should include entire packet except the first 8 chars:  'OCCDq12b'

  my $fixc2l2 = (exists($self->{FIXC2L2}) && $self->{FIXC2L2});
  my $time = BJParser::parseStrip ("L", \$pkt); # GPS seconds

  # PPS == 'packets per second' (50, 100, etc)
  my ($conid, $prn, $ant, $fdma_chan, $obstype, $trackstat, $pps, $primary_chan,
      $height, $azimuth, $latitude, $longitude, $beamOffset1, $beamOffset2,
      $signal_type1, $snr1, $range1, $restau1,
      $signal_type2, $snr2, $range2, $restau2,
      $fractime, $fitphase1, $fitphaserate1, $fitphaseraterate1, $fitphase2, $fitphaserate2, $fitphaseraterate2
    ) = BJParser::parseStrip ("CCCcCCCC ssssss CSdd CSdd ddddddd", \$pkt);
  $conid = chr($conid); # convert to character (G = GPS, R = GLONASS, etc)
  #$ant           = chr($ant);            # This is an ASCII character!

  while ($prn > 100) { $prn -= 100 }  # Sometimes GLONASS PRNs are in the 100s or 200s

  # GPS L2C flag.  If no external map is available, assume L2C rather than L2P
  # suspicious Jun
  my $isgpsl2c  = $conid eq 'R' ? 0
                : ((keys %{$self->{L2CMAP}} > 0) && !$self->{L2CMAP}{$prn}) ? 0
                : 1;

  push (@{$self->{OCCTABLE}{"$conid$prn"}}, [$time, $ant, $height, $azimuth, $latitude, $longitude]);

  # suspicious Jun
  my ($f1, $f2);
  if ($conid eq 'R') {
    #($f1, $f2) = SatInfo::glonass_freq($fdma_chan - 8);
    ($f1, $f2) = SatInfo::glonass_freq($fdma_chan);
  } elsif ($conid eq 'C') {
    ($f1, $f2) = ($B1, $B2);  # Beidou frequencies
  } elsif ($conid eq 'E') {
    #($f1, $f2) = ($E1, $E5b); # Galileo frequencies
    ($f1, $f2) = ($E1, $E5a); # Jun modified
  } else {
    ($f1, $f2) = ($F1, $F2);  # GPS frequencies
  }
  my $lamb1 = $CLIGHT/$f1;
  my $lamb2 = $CLIGHT/$f2;

  # Read in high rate data efficiently (??? benchmark this)
  ${$I1->get_dataref} = substr ($pkt, 0, 200, ''); $I1->upd_data; $I1->bswap2;
  ${$Q1->get_dataref} = substr ($pkt, 0, 200, ''); $Q1->upd_data; $Q1->bswap2;
  ${$I2->get_dataref} = substr ($pkt, 0, 200, ''); $I2->upd_data; $I2->bswap2;
  ${$Q2->get_dataref} = substr ($pkt, 0, 200, ''); $Q2->upd_data; $Q2->bswap2;

  my $dfaz1 = atan2($Q1,$I1) * $RAD2CYC/$f1; # seconds (rad * cyc/rad * sec/cyc = sec)

  # Early COSMIC-2 data had L2 problems that showed up as phasor flips like nav bit flips.  Attempt
  # to fix these via 2 quadrant extraction
  my $dfaz2;
  if ($fixc2l2 && ($conid eq 'G') && !$isgpsl2c) {
    # Use atan (2 quadrant) if GPS L2P and requested to do so.
    $dfaz2 = atan(double($Q2)/double($I2)) * $RAD2CYC/$f2; # seconds
  } else {
    # Otherwise, use 4 quadrant
    $dfaz2 = atan2($Q2,$I2) * $RAD2CYC/$f2;
  }

  my $dt    = (1/$pps) * sequence($pps); # 0, 0.01, 0.02, ..., 0.99
  my $mdl1  = ($fitphase1 + $fitphaserate1*$dt + $fitphaseraterate1*$dt**2/2); # phase model (sec)
  my $mdl2  = ($fitphase2 + $fitphaserate2*$dt + $fitphaseraterate2*$dt**2/2);
  my $hrfaz1 = $CLIGHT * ($mdl1 - $dfaz1); # full phase, meters
  my $hrfaz2 = $CLIGHT * ($mdl2 - $dfaz2);
  $mdl1 *= $CLIGHT; # Convert phase model to meters as well.  D. Hunt 3/16/2018

  # From emails from Tom Meehan dated 5/19/2018 and 6/25/2018:
  # T = ~0.01 sec (small dead time for correlator R/W; closer to 0.0099)
  # sineFactor = 0.75
  # sampleRate = 20456000 samples/sec
  # noiseIQ = std deviation over integration  of I,Q correlator values = sqrt [ 2 * sampleRate * T *  sineFactor ]
  # (noiseIQ ~ 552 counts)
  # Amplitude is sqrt ( I^2 + Q^2 )
  # SNRv for each 10 msec point is  Amplitude  / noiseIQ.
  # However we reduce the number of bits to fit I & Q into 16 bit integers by dividing by an amplitude scale factor, AmpScale.
  # AmpScale is passed as a variable in the lower 6 bits of the ObsType field.
  # AmpScale is set to 12 for beamformed signals.
  #
  # 10msec SNRv = AmpScale * (sqrt(I**2 + Q**2) / 552 (L1, L2C, Glonass L1, L2)
  # 10msec SNRv = AmpScale * (sqrt(I**2 + Q**2) / 1299 (GPS L2P)

  # High rate CA SNR, P2 SNR
  $obstype = 0x04 if ($obstype == 0);     # Default to non beam-formed signal (helps for old DITL datasets)
  my $ampscale = ($obstype & 0b00111111); # should be 12 for beamformed signals, 4 for non beamformed signals
  my $bf_on    = ($obstype & 0b00001111) == 0x0C; # Beam forming on

  my $noiseIQ1 = 552;
  # suspicious Jun
  #my $noiseIQ2 = ($signal_type2 == 0x30) ? 1299 : 552;     # 0x30 is L2P

  # Modified Jun: update $noiseIQ2 for Galileo S5Q
  #my $noiseIQ2;
  #if ($signal_type2 == 0x30) {
  #	  $noiseIQ2 = 1299.0;
  #}
  #elsif ($signal_type2 == 0x91) {
  #	  $noiseIQ2 = 957.0;
  #}
  #else {
  #	  $noiseIQ2 = 552.0;
  #}

  # Modified Jun: update $noiseIQ2 for GNSS E/G/R
  my $noiseIQ2;
  if ($signal_type2 == 0x30) {  #GPS L2P
	  $noiseIQ2 = 2253.0;
  }
  elsif ($signal_type2 == 0x2c) { #GPS L2L
	  $noiseIQ2 = 957.0;
  }
  elsif ($signal_type2 == 0x03) { #Glonass L2C
	  $noiseIQ2 = 552.0;
  }
  elsif ($signal_type2 == 0x91) { # Galileo  S5Q
	  $noiseIQ2 = 957.0;
  }
  else {
	  $noiseIQ2 = 552.0;
  }



  # Email from Tom Meehan on 7/23/2019 says that noiseIQ for beam forming should be sqrt(3) larger.
  if ($bf_on) {
    $noiseIQ1 *= $bf_factor; # sqrt(3)
    $noiseIQ2 *= $bf_factor; # sqrt(3)
  }

  my $snrfact  = sqrt($I1->nelem); # 10.  Make 100 Hz SNRs have same scale as 1 second SNRs.  D. Hunt 7/5/2018
  my $hrsnr1   = $snrfact * $ampscale/$noiseIQ1 * sqrt($I1**2 + $Q1**2) * 10; # multiply v/v by 10 to use more of a ushort's range
  my $hrsnr2   = $snrfact * $ampscale/$noiseIQ2 * sqrt($I2**2 + $Q2**2) * 10;

  $dt += $fractime; # full delta time

  #
  ## Now format this as an opnGns data record
  #
  my $satid   = sprintf "%1s%02d", $conid, $prn; # eg 'G22'
  my ($svn, $freq_slot) = $satinfo->find_svn($satid, $time);
  my $status  = 0x01; # open loop (COSMIC-2 is always open loop)
  $status    += (0x02 * $bf_on); # Beam forming flag, 0b0000 00X0:  0 = no beam forming 1 = beam forming
  my $rec     = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $pps);

  # Format Epoch Data Record (v2)  Use @lrdesc and @hrdesc defined at a package level for speed
  # Add low rate record count, high rate record count, then low rate descriptors and high rate descriptors
  # suspicious Jun
  # Jun added E
  if ($conid eq 'G' && (($signal_type2 == 0x30) || !$isgpsl2c)) { # L2P (0x30) instead of L2CL (0x2C)
    $rec .= pack ('CC', scalar(@lrdesc), scalar(@hrdesc_l2p_gps)) . join ('', @lrdesc) . join ('', @hrdesc_l2p_gps);
  } elsif ($conid eq 'R') {
    $rec .= pack ('CC', scalar(@lrdesc), scalar(@hrdesc_l2c_glonass)) . join ('', @lrdesc) . join ('', @hrdesc_l2c_glonass);
  } elsif ($conid eq 'E') {
    $rec .= pack ('CC', scalar(@lrdesc), scalar(@hrdesc_l2c_galileo)) . join ('', @lrdesc) . join ('', @hrdesc_l2c_galileo);
  } else { # default to GPS L2C
    $rec .= pack ('CC', scalar(@lrdesc), scalar(@hrdesc_l2c_gps)) . join ('', @lrdesc) . join ('', @hrdesc_l2c_gps);
  }

  my $range_meters = ($range1+$restau1) * $CLIGHT;
  my $range_model_meters = $range1 * $CLIGHT;
  my $lrdata  = pack ("(dddddddd)>", $range_meters, $range_model_meters, $height, $azimuth, $latitude, $longitude, $beamOffset1, $beamOffset2);
  my $tf      = long((($dt+0.5)/1e-7)->rint); # this must be positive!!! if not, fiddle the trailer time offset value to be > 0.5
  my $mostsig = byte(($tf & 0x00ff0000) >> 16);
  my $midsig  = byte(($tf & 0x0000ff00) >>  8);
  my $lowsig  = byte ($tf & 0x000000ff);
  my $hrvals  = cat (double($mostsig), double($midsig), double($lowsig), $hrfaz1, $hrfaz2, $mdl1, $hrsnr1, $hrsnr2);
  $hrvals->inplace->setbadtonan;
  my $hrrec   = pack ("(CCCdddSS)>" x $pps, $hrvals->xchg(0,1)->flat->list); # 31 byte records

  # Output hash
  # This is the hash format compatible with Level1aTools::sort_opnGns_v2
  $self->{OUTFILES}{'hirate'}{$ant}{$conid}{$prn}{$time} = ($rec . $lrdata . $hrrec);

  return $time;

}


#/**----------------------------------------------------------------------
# @sub       parse_SCNThrat
#
# Subroutine to parse the new C-2 scintillation data record:  SCNThrat
# Input is the packet, with 'SCNThrat' removed and the CRC removed
# (this routine is only called if the CRC is good).
#
# Output is a structure in $self containing the decoded SNR or phase values.
# These values are written to an opnScn file by trigZero2one.pl.
#
# @parameter  $self -- BJParser object
# @           $pkt -- Contents of SCNThrat packet (without 'OBSDqfit')
# @                   as a perl string
# @return     Scintillation parameters appended to a file or structure:
# @           Structure: $self->{OUTFILES}{$file}
# @           High rate data added to a structure:
# @           $self->{OUTFILES}{'SCNThrat'}{$ant}{$fullprn}{$time}{$obstype} = $values;
# @           where $values are the L1 or L2 S4 or sigma phi high rate data.
# @exception  Will toss an exception in case of real trouble.
# ----------------------------------------------------------------------*/
sub parse_SCNThrat {

  my $self = shift;
  my $pkt  = shift;  # Should include entire packet except the first 8 chars:  'OBSDq12a'

  my $time = BJParser::parseStrip ("L", \$pkt); # 4 bytes, big endian unsigned integer

  my $tmpdir = $self->{TMPDIR};

  # C = unsigned byte, c = signed byte, S = unsigned short, d = double, f = float.  See 'perldoc pack'
  my ($conid, $prn, $ant, $fdma_chan, $obstype, $num_points, $primary_chan, $signal_type,
      $snr, $x, $v, $a, $S4, $sig_phi, $scale_resids)
      = BJParser::parseStrip ("CCCcCSCC Sdddfff", \$pkt);
  $conid = chr($conid); # convert to character
  $obstype = chr($obstype);
  while ($prn > 100) { $prn -= 100 }  # Sometimes GLONASS PRNs are in the 100s or 200s
  my $fullprn = sprintf "%1s%02d", $conid, $prn;

  #print "SCNThrat $time $fullprn $obstype\n";

  #
  ## There is a format difference for GLONASS under firmware version 4.3.3.  Check for this and handle differently
  #

  # Email from Tom Meehan, 11/12/2019:

  # We have been using the TGScintHighRate11A packet type
  # for high-rate scintillation data with both GPS and GLONASS.
  # V4.3.3 modifies only the GLO data that fills out this packet. The
  # packet is unchanged.

  # New format:
  #   Packets are output *every 1-second instead of every 10-secs*. L1 packet, then L2 packet.
  #   100 Hz I,Q data (200 raw I,Q, values scaled to fit into 8-bit int).
  #
  #   Instead of 1,2,3,4; ObsType denotes whether L1 (Obstype=5) or L2 (ObsType=6).
  #   This ObsType notation should assure backwards
  #   compatibility with earlier scint data.
  #
  #   SigPhi field now contains fractional secs at start of interval.
  #   Increments by exactly 0.01 sec as with Atm Occ packets (TGAtmoOcc12B).
  #   The x,v,a fields are used to create a phase model at 100 Hz as we do with
  #   TGAtmoOcc12B data packet.
  #
  #   High rate phase is created by adding the model to the fractional phase
  #   from atan(Q/I) while keeping track of cycles. This format should allow
  #   adding nav-bit half-cycles as with atmo-RO.
  #
  #   High rate SNR is just the RSS of the scaled I&Q values. Where
  #   scaleResids * Resids is used, divide RSS value by noise
  #   to get 1-sec SNR.
  #
  #   10 msec noise (N10ms) is sqrt(2 * 20456000 * 0.75 * 0.01).
  #   1-sec SNRv = sqrt (100 ) * RSS / N10ms.
  #
  #   For C2, 10 msec RSS of 35000 I get a 1-sec SNRv of ~632 V/V.
  #
  #   Doug, please note the scale factor of 4 or 12 (BF= no or yes) present in atmo-RO
  #   data is not present here. The scaleResids field encompasses the entire
  #   adjustment from 8-bits in packet to 32-bits in open-loop DSP correlator.

  # If new format (v4.3.3) data
  if ($obstype > 4 && $conid eq 'R') {

    my ($f1, $f2) = SatInfo::glonass_freq($fdma_chan);
    my $freq   = ($obstype == 5) ? $f1 : $f2;
    my $freqid = ($obstype == 5) ? 'L1' : 'L2';

    # Read in high rate data
    my $IQ = pdl(BJParser::parseStrip ("c"x200, \$pkt));
    my $I  = $IQ->slice("0:-1:2") * $scale_resids;
    my $Q  = $IQ->slice("1:-1:2") * $scale_resids;

    # Is this the right handling of $scale_resids???

    my $dfaz  = atan2($Q,$I) * $RAD2CYC/$freq; # seconds (rad * cyc/rad * sec/cyc = sec)

    my $pps   = 100; # Data frequency
    my $dt    = (1/$pps) * sequence($pps); # 0, 0.01, 0.02, ..., 0.99
    my $mdl   = ($x + $v*$dt + $a*$dt**2/2); # phase model (sec)
    my $hrfaz = $CLIGHT * ($mdl - $dfaz); # full phase, meters
    $mdl     *= $CLIGHT; # Convert phase model to meters as well.  D. Hunt 3/16/2018

    my $N10ms = sqrt(2 * 20456000 * 0.75 * 0.01); # make this a constant at package level?
    my $rss   = sqrt($I**2 + $Q**2);
    my $hrsnr = $rss/$N10ms * 10; # multiply v/v by 10 to use more of a ushort's range

    # print "SCNThratII: $fullprn $time $freqid $sig_phi\n";
    $dt += $sig_phi; # For new packet, sig_phi is interpreted as fractime.  Full delta time

    # Add scintillation data for this observation type, PRN and epoch to a data structure
    # for writing to an output file later.
    $self->{OUTFILES}{'SCNThrat433'}{$ant}{$fullprn}{$time}{$freqid} = [$dt, $hrfaz, $mdl, $hrsnr, $signal_type];

    return $time;
  }


  # If new format (v4.3.4) GPS data
  # The obstype 7 and 8 data are L1 or L2 one-byte 50 Hz phase residuals in one-second packets.
  # They must be combined with the 10-second obstype 1 and 2 packets to get the SNRs :-(
  # See emails from Tom Meehan dated 2/26/2020.
  if ($obstype > 6 && $conid eq 'G') {

    my $freq   = ($obstype == 7) ? $F1 : $F2;
    my $freqid = ($obstype == 7) ? 'L1' : 'L2';

    # Read in high rate data
    my $resid = pdl(BJParser::parseStrip ("c"x50, \$pkt)) * $scale_resids; # cycles

    my $pps   = 50; # Data frequency
    my $dt    = (1/$pps) * sequence($pps) - 0.5; # -0.5, -0.48, ..., 0, ..., +0.48
    my $mdl   = ($x + $v*$dt + $a*$dt**2/2); # phase model (sec)
    my $hrfaz = ($mdl + $resid/$freq) * $CLIGHT; # full phase, meters
    $mdl     *= $CLIGHT; # Convert phase model to meters as well.  D. Hunt 3/16/2018

    # Add scintillation data for this observation type, PRN and epoch to a data structure
    # for writing to an output file later.
    $self->{OUTFILES}{'SCNThrat434'}{$ant}{$fullprn}{$time}{$freqid} = [$dt, $hrfaz, $mdl, $signal_type];

    return $time;
  }


  # Some GLONASS simulator scintillation data has bogus high residual scales.
  # Set these values to NaN to flag them as bad for the scintillation phase routines.  D. Hunt 6/17/2019
  # This represents an unrealisticly large jump in phase (about 30 meters/sec) per increment
  # of residual.
  my $max_residual = 1e-9;
  if ($obstype >= 3 && $obstype < 5 && abs($scale_resids) > $max_residual){
    print "$conid$prn $time, scale_resids too large: $scale_resids.  Setting to NaN\n";
    $scale_resids = $nan;
  }

  # print "$conid, $prn, $ant, $fdma_chan, $obstype, $num_points, $primary_chan, $signal_type, $snr, $x, $v, $a, $S4, $sig_phi, $scale_resids\n";
  $prn -= 100 while ($prn > 100); # GLONASS PRNs can either be in the 100's or the 200's!

  # Note that Galileo and Beidou paths are for use only in the Day In The Life test.  These
  # constellations are not yet fully supported.
  my ($f1, $f2, $scale_snr);
  $scale_snr = 1;
  if ($conid eq 'R') {
    $scale_snr = 55/8; # will change with next software release (may change to 1) !!!
  }

  my $resids = pdl(BJParser::parseStrip ("c"x$num_points, \$pkt));

  #
  ## Append scintillation parameters to the 'scinhr[1234].txt' file
  ## Note that $sig_phi is always zero.  'S4' means:
  ## L1_S4       for obstype = 1
  ## L2_S4       for obstype = 2
  ## L1_Sigmaphi for obstype = 3
  ## L2_Sigmaphi for obstype = 4
  #

  my $outrec_scin = sprintf("%ld SCIN %1s%02d %02x %f\n", $time, $conid, $prn, $ant, $S4);
  my $file = "scin$obstype.txt";
  $self->{OUTFILES}{$file} .= $outrec_scin;  # NOTE:  This is currently not used!!!

  #
  ## Decode compressed scintillation data.  See Angie Dorsey's email of 11/21/2016
  #

  my $tf     = 10 * sequence($num_points)/$num_points; # Time fraction
  my $model  = $x + $v*$tf + 0.5*$a*$tf**2;            # Quadratic model
  my $values = $model + $scale_resids * $resids;
  $values   /= $scale_snr if ($obstype <= 2);          # Only scale SNRs, not phases

  # Not sure what units SNRs (obstype 1,2) are in.  Convert phases (obstype 3,4) to meters
  # from seconds
  $values *= $CLIGHT if ($obstype > 2);

  if ($obstype <= 2 && any($values < 0)) {
    $values->where($values < 0) .= 0;
    print "SCNThrat version 1 L$obstype SNR < 0:  T$time PRN$fullprn.  Setting to 0\n";
  }

  # If these are GPS phases, they still have the intermedidate frequency in them.
  # Take this out (based on emails from Tom Meehan, 3/13/2019)
  if ($obstype == 3 && $conid eq 'G') {
    my $l1IF = 308_000; # Hz
    my $start_phase = $values->at(0);
    my $dop = append(0, ($values - $values->rotate(1))->slice("1:-1"));
    $dop = -1 * ($dop - ($l1IF * $LAMB1)/($num_points/10)); # Not sure why -1.  Works when compared with OBSDq12a data
    $values = dcumusumover($dop) + $start_phase;
  } elsif ($obstype == 4 && $conid eq 'G') {
    my $l2IF = 240_000; # Hz
    my $start_phase = $values->at(0);
    my $dop = append(0, ($values - $values->rotate(1))->slice("1:-1"));
    $dop = -1 * ($dop - ($l2IF * $LAMB2)/($num_points/10)); # Not sure why -1.  Works when compared with OBSDq12a data
    $values = dcumusumover($dop) + $start_phase;
  }

  # Output for debug
  #open my $fh, ">>scint_OT$obstype";
  #wcols "%10.5f  %4d", $values, $resids, $fh;
  #close $fh;

  # Add scintillation data for this observation type, PRN and epoch to a data structure
  # for writing to an output file later.
  $self->{OUTFILES}{'SCNThrat432'}{$ant}{$fullprn}{$time}{$obstype} = [$values, $signal_type];

  return $time;

}


#/**----------------------------------------------------------------------
# @sub       parse_RCVMadct
#
# Subroutine to parse the C-2 TriG State of Health record
# Input is the packet, with 'RCVMadct' removed and the CRC removed
# (this routine is only called if the CRC is good).
#
# Output is a structure in $self containing the decoded SOH values.
#
# @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_RCVMadct {

  my $self = shift;
  my $pkt  = shift;  # Should include entire packet except the first 8 chars:  'RCVMadct'

  my $time = BJParser::parseStrip ("L", \$pkt); # 4 bytes, big endian unsigned integer

  # Loop, extracting values, until packet is empty
  while (1) {
    last if (length($pkt) <= 0);
    my ($val, $type) = unpack ("f>a", substr($pkt, 0, 5, ''));
    my ($name) = ($pkt =~ m/(\w+)\x0/);
    $pkt = substr($pkt, length($name)+1);
    # print "T$time, $type$name = $val\n";
    my $key = $type . sprintf ("%03d", $name);
    $self->{OUTFILES}{'RCVMadct'}{$time}{$key} = $val;
    last if (!defined($name));
  }

  return $time;

}

#/**----------------------------------------------------------------------
# @sub       parse_RCVMcmdr
#
# Subroutine to parse the C-2 TriG command response packet.
# Input is the packet, with 'RCVMcmdr' removed and the CRC removed
# (this routine is only called if the CRC is good).
#
# Output is a structure in $self containing the decoded SOH values.
#
# @parameter  $self -- BJParser object
# @           $pkt -- Contents of RCVMcmdr packet (without 'RCVMcmdr')
# @                   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_RCVMcmdr {

  my $self = shift;
  my $pkt  = shift;  # Should include entire packet except the first 8 chars:  'RCVMcmdr'

  my $time = $self->{LASTTIME};

  my $date    = defined($self->{BYDATE})  ? '.'.TimeClass->new->set_gps($time)->get_yrdoy_gps : '';
  my $lv0name = defined($self->{ADDNAME}) ? '.'.Dplib::simplefn($self->{INFILE}) : '';
  my $tmpdir = $self->{TMPDIR};

  my ($accepted, $library_id, $cmd_code) = unpack "Ca4a4", substr ($pkt, 0, 9, '');
  my $status_code = BJParser::parseStrip ("L", \$pkt);

  $pkt =~ s/\x00$//; # get rid of trailing null
  my $message = $pkt;

  # Append " log> " to all lines in the log message
  my $outrec = " cmdr> T$time ACC:$accepted $library_id$cmd_code $status_code $message\n";

  my $file = "cmdr.txt$date$lv0name";

  $self->{OUTFILES}{$file} .= $outrec;

  return;

}


#/**----------------------------------------------------------------------
# @sub       dump_occs
#
# Generate a list of occultation times and locations based on OCCDq12a/b tracking info.
#
# @parameter  $self      Parser object
# @           $outfile   Name of file to write to
# @exception  If there is a problem...
# ----------------------------------------------------------------------*/
sub dump_occs {

  my $self    = shift; # parser object
  my $outfile = shift;

  open my $fh, '>', $outfile or die "Cannot open $outfile for writing";

  foreach my $prn (sort keys %{$self->{OCCTABLE}}) {
    my @obs = @{$self->{OCCTABLE}{$prn}};
    foreach my $ob (@obs) {
      printf {$fh} "%2s %15.1f %1d %8.3f %8.3f %8.2f %8.2f\n", $prn, @$ob;
    }
  }
  close $fh;
  print "FILE CREATION: ", Dplib::simplefn($outfile), "\n";
}


#/**----------------------------------------------------------------------
# @sub       plot_occs
#
# Generate a plot of occultation locations based on OCCDq12a/b tracking info.
#
# @parameter  $self      Parser object
# @           $outfile   Name of file to write to
# @exception  If there is a problem...
# ----------------------------------------------------------------------*/
sub plot_occs {

  my $self    = shift; # parser object
  my $outfile = shift;

  #
  ## Create map plot
  #

  my @pagesize = (1200, 600);  # pixels
  my $planet = PDL::Planet->new->create(zeroes(PDL::byte,3,@pagesize)+255);

  # Create a background map to draw upon
  my $pl = PDL::Graphics::PLplot::Map->new (DEV  => 'memcairo',
                                            MEM  => $planet->rgb,
                                            PAGESIZE => [@pagesize]);

  # Draw the map
  my $viewport = [0.05, 0.9, 0.15, 0.85];
  $pl->worldmap (PROJECTION => 'LINEAR',
                 MAPBOX => [-180, 180, -90, 90],
                 VIEWPORT => $viewport,
                 JUST => 0, # Important! Otherwise the VIEWPORT is not used.
                 RESOLUTION => 'low');

  $pl->close;

  $pl = PDL::Graphics::PLplot::Map->new (DEV  => 'memcairo',
                                         MEM  => $planet->rgb,
                                         PAGESIZE => [@pagesize]);

  my %occ_count;
  my $all_h = PDL::null;
  foreach my $prn (sort keys %{$self->{OCCTABLE}}) {
    $all_h = $all_h->append(pdl($self->{OCCTABLE}{$prn})->slice("(2),:"));
  }
  my ($minheight, $maxheight) = $all_h->minmax;
  foreach my $prn (sort keys %{$self->{OCCTABLE}}) {

    my $recs = pdl($self->{OCCTABLE}{$prn});
    my $t    = $recs->slice("(0),:");
    my $gaps = which ($t - $t->rotate(1) > 1);
    my @idx  = (0, $gaps->list, $t->nelem - 1);
    my $symbol = ($prn =~ /G/) ? 851 : 852;
    for (my $i=0;$i<@idx-1;$i++) {
      my $si  = $idx[$i];
      my $ei  = $idx[$i+1];
      my $h   = $recs->slice("(2),$si:$ei");
      my $lat = $recs->slice("(4),$si:$ei");
      my $lon = $recs->slice("(5),$si:$ei");
      $lon->where($lon > 180) -= 360; # convert to -180 to 180 longitude
      $pl->xyplot ($lon, $lat,
                   BOX => [-180, 180, -90, 90],
                   XBOX => '', YBOX => '', SYMBOL => $symbol, SYMBOLSIZE => 0.5,
                   VIEWPORT => $viewport,
                   PALETTE => 'REDGREEN', PLOTTYPE => 'POINTS', COLORMAP => $h, ZRANGE => [$minheight, $maxheight]);
      $occ_count{substr($prn,0,1)}++;
      # printf "%3s lat %3d lon %3d\n", $prn, $lat->avg, $lon->avg;
    }
  }
  (my $title = $outfile) =~ s/\.png//;
  $title = Dplib::simplefn($title);

  $pl->text('Occultations from', TEXTPOSITION => ['T', 3.5, 0.5, 0.5], CHARSIZE => 1, COLOR => 'BLACK');
  $pl->text($title, TEXTPOSITION => ['T', 2.0, 0.5, 0.5], CHARSIZE => 1, COLOR => 'BLACK');
  $pl->text('Square = GPS, Triangle = GLONASS', CHARSIZE => 0.8, TEXTPOSITION => ['B', 1.0, 0.0, 0.0]);

  $pl->text("Copyright (C) 1999-2016 University Corporation for Atmospheric Research",
            TEXTPOSITION => ['L', 6.0, 0, 0,], COLOR => 'BLUE', CHARSIZE => 0.4);
  $pl->text("All rights reserved",
            TEXTPOSITION => ['L', 4.0, 0, 0,], CHARSIZE => 0.4);

  $pl->colorkey (pdl($minheight, $maxheight), 'v', YBOX => 'BNSTI', YLAB => 'Straight Line Height, km',
                 CHARSIZE => 0.6, VIEWPORT => [0.95, 0.98, 0.15, 0.85]);

  $pl->close;

  open my $fh, '>', "$outfile\_map.png" or die "Cannot open $outfile for writing";
  print {$fh} $planet->write_png;
  close $fh;
  print "FILE CREATION: ", Dplib::simplefn("$outfile\_map.png"), "\n";

  #
  ## Occultation height vs. time plot
  #

  undef $pl;
  $pl = PDL::Graphics::PLplot->new(DEV => 'pngcairo', FILE => "$outfile\_height.png", PAGESIZE => [@pagesize], SUBPAGES => [1,2]);

  my $all_t = PDL::null;
  foreach my $prn (sort keys %{$self->{OCCTABLE}}) {
    $all_t = $all_t->append(pdl($self->{OCCTABLE}{$prn})->slice("(0),:"));
  }
  my ($mint, $maxt) = $all_t->minmax;

  my @colors   = (qw(BLACK GREEN WHEAT BLUE RED AQUAMARINE GREY BLUEVIOLET YELLOW PINK BROWN CYAN TURQUOISE));
  my $coloridx = 0;
  my @first = (1,1,1);
  my @vp = (0,
            [0.10, 0.90, 0.05, 0.90],   # top
            [0.10, 0.90, 0.20, 0.95]);  # bottom
  foreach my $prn (sort keys %{$self->{OCCTABLE}}) {

    my $panel = ($prn =~ /G/) ? 1 : 2;
    my $recs  = pdl($self->{OCCTABLE}{$prn});
    my $t     = $recs->slice("(0),:") - $mint;
    my $h     = $recs->slice("(2),:");
    my $xbox  = my $ybox = '';
    $xbox     = 'BCST'  if ($first[$panel] && $panel == 1);
    $xbox     = 'BCNST' if ($first[$panel] && $panel == 2);
    $ybox     = 'BCNST' if ($first[$panel]);
    $first[$panel] = 0;
    $pl->xyplot ($t, $h,
                 COLOR    => $colors[$coloridx],
                 BOX      => [0, $maxt-$mint, $minheight, $maxheight],
                 SYMBOL   => 850, SYMBOLSIZE => 0.5,
                 XBOX     => $xbox, YBOX => $ybox,
                 VIEWPORT => $vp[$panel],
                 SUBPAGE  => $panel,
                 PLOTTYPE => 'POINTS');
    $coloridx = (($coloridx + 1) % scalar(@colors));
  }

  ($title = $outfile) =~ s/\.png//;
  $title = Dplib::simplefn($title);

  $pl->text("Occultations from $title", TEXTPOSITION => ['T', 0.3, 0.5, 0.5], SUBPAGE  => 1,
            CHARSIZE => 1, COLOR => 'BLACK');

  my $startstamp = TimeClass->new->set_gps($mint)->get_stamp_gps;
  $pl->text("Seconds since $startstamp", TEXTPOSITION => ['B', 3.0, 0.5, 0.5], SUBPAGE  => 2,
            CHARSIZE => 1, COLOR => 'BLACK');

  $pl->text(sprintf ("Total occultations: GPS = %3d, GLONASS = %3d",  @occ_count{'G', 'R'}),
            CHARSIZE => 1, COLOR => 'BLUE', TEXTPOSITION => ['B', 3.0, 0.0, 0.3]);

  $pl->text("HSL (km) GPS",     TEXTPOSITION => ['L', 3.0, 0.4, 0.5],
            SUBPAGE  => 1, CHARSIZE => 1, COLOR => 'BLACK');
  $pl->text("HSL (km) GLONASS", TEXTPOSITION => ['L', 3.0, 0.5, 0.5],
            SUBPAGE  => 2, CHARSIZE => 1, COLOR => 'BLACK');

  $pl->close;

  print "FILE CREATION: ", Dplib::simplefn("$outfile\_height.png"), "\n";

}

1;