# ## Copyright (c) 1995-2020 University Corporation for Atmospheric Research ## All rights reserved # my $pkgdoc = <<'EOD'; #/**---------------------------------------------------------------------- # @file Level1aTools.pm # # Tools for converting raw data into level 1a data formats. # # @author Doug Hunt # @since 10/07/2003 # @cdaacTask no # @version $URL: svn://svn1.cosmic.ucar.edu/trunk/src/BJutils/Level1aTools.pm $ $Id: Level1aTools.pm 22447 2016-12-01 03:58:34Z dhunt $ # -----------------------------------------------------------------------*/ EOD package Level1aTools; use strict; use warnings; use FileHandle; use PDL; use BJConstants; use File::Basename; use File::Temp; #/**---------------------------------------------------------------------- # @sub blackJack2rinex # # Convert an ASCII file in the ad-hoc BJfmtl format to RINEX 2.20 format # similar to a podObs file. # # @parameter $infile Input ASCII file. # @ $outfile Output Hatenaka compressed, gzipped RINEX file. # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub blackJack2rinex { my $infile = shift; my $outfile = shift; my $antid = shift; my $leoname = shift; my $mission = shift; my $interval = shift || 0; # if no interval specified, do not restrict based on interval my $options = shift || {compress => 1}; my ($outmem, $rinex_body, $rinex_header); my ($fileyr, $filedoy) = ($outfile =~ /(\d\d\d\d)\.(\d\d\d)/); if ($options->{rinex_header}) { $rinex_body = $options->{rinex_body}; $rinex_header = $options->{rinex_header}; } else { $$rinex_body = ''; $$rinex_header = ''; } my $filegps = TimeClass->new->set_yrdoyhms_gps($fileyr, $filedoy, 0, 0, 0)->get_gps; open my $infh, '<', $infile or die "Cannot open input file $infile"; my $duplicate_count = 0; # count of duplicate records for a given time stamp my $endtime = 0; my %rec = (); # hash by PRN of data for one time step. Eash entry is another hash # with all the data types shown below in %hash. my %lastRec = (); # hash with all data for the previous record. Useful in determining # which duplicate record to use if there is more than one for a given # PRN. The rule is that tracks should, when possible, continue to # use the same obstype. my @times = (); RECORD: while (<$infh>) { # use nifty hash slice for record my %hash = (); @hash{'time', 'obsd', 'prn', 'antidx', 'obstype', 'interval', 'pcktsize', 'ca', 'cachan', 'casnr', 'caphase', 'capseudo', 'p2', 'p2chan', 'p2snr', 'p2phase', 'p2pseudo', 'p1', 'p1chan', 'p1snr', 'p1phase', 'p1pseudo'} = split; next unless (defined ($hash{obsd}) && $hash{obsd} =~ /OBSD/); # skip nav data next unless ($hash{time} =~ /\d+/); # make sure time is numeric next unless (abs($filegps - $hash{time}) < 604800); # skip records with times way off expected time. D. Hunt 4/11/2007 # If requested, allow 1 second times that are a multiple of 10 seconds if ($interval && $options->{allinterval}) { next if (($hash{time} % $interval) != 0); # reject records that are not a multiple of the requested interval. D. Hunt 3/20/2019 } elsif ($interval) { next if ((hex($hash{interval}) != $interval)); # reject records with data rate != to that requested D. Hunt 7/27/2004 } $hash{prn} = substr($hash{prn},3,2); # get rid of 'prn' in eg 'prn08' my $prn = 'G' . $hash{prn}; # Add G to PRN number $hash{cachan} = hex($hash{cachan}); # check for out of order records if ($hash{time} < $endtime) { warn "Out of order record: time = $hash{time}, prn = $hash{prn}"; next RECORD; } # write rinex record when all prns for this epoch gathered if ($endtime != $hash{time} && $endtime != 0) { push (@times, $endtime); $$rinex_body .= writeRinex (\%rec, $endtime); %lastRec = %rec; undef %rec; } # If this is a duplicate entry for this PRN, we decide if we should use this record # or the currently existing record. We want to keep the obstype the same as the last record # if possible. if ($rec{$prn}) { # Comment these print statements out because they make too # much noise in the output files. Now there is the possibility # of processing 'double dump' cosLv0 files (dump_new, dump_prev) # so there will be many of these messages. D. Hunt 9/25/2006 $duplicate_count++; # Don't add %hash to %rec if this would cause a break in obstypes next RECORD if (defined ($lastRec{$prn}) && ($lastRec{$prn}{obstype} ne $hash{obstype})); } # Hash of prns per time step. Used to save obstypes from the last time step # in order to help decide which observations to use if there are two observations # for a given time step and PRN. $rec{$prn} = \%hash; $endtime = $hash{time}; } # Write last record $$rinex_body .= writeRinex (\%rec, $endtime); push (@times, $endtime); close ($infh); my ($tstart, $tend) = @times[0,-1]; return if ($tend == 0); # no records written, bail out if (exists ($options->{rinex_header})) { $$rinex_header = writeRinexHeader ($antid, $leoname, $tstart, $tend, $mission); return; # $$rinex_body is already filled up } my $outfh; if($options->{compress}){ open $outfh, "| rnx2crx - -s -e 100 | gzip - > $outfile"; }else{ $outfile =~ s/_crx\.gz$//i; open $outfh, "> $outfile"; } print {$outfh} writeRinexHeader ($antid, $leoname, $tstart, $tend, $mission); print {$outfh} $$rinex_body; close ($outfh); print "Found $duplicate_count records for duplicate times\n" if ($duplicate_count); return; } #/**---------------------------------------------------------------------- # @sub 012A_to_rinex # # Convert an ASCII '012A' file from the TRIG in a perl hash structure # to the RINEX 2.11 format. # # @parameter $infile Input ASCII file. # @ $outfile Output Hatenaka compressed, gzipped RINEX file. # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub twoA_to_rinex { my $data = shift; my $outfile = shift; my $antid = shift; my $leoname = shift; my $mission = shift; my $options = shift || {COMPRESS => 1}; my $rnxver = $options->{RNXVER} // 2.20; # Default to old-style, less-correct RINEX files which work my $firmware_version = $options->{FIRMWARE_VERSION} // 'v434+'; # default to new firmware my ($outmem, $rinex_body, $rinex_header); my ($fileyr, $filedoy) = ($outfile =~ /podCrx_(\d\d\d\d)\.(\d\d\d)/); if ($options->{rinex_header}) { $rinex_body = $options->{rinex_body}; $rinex_header = $options->{rinex_header}; } else { $$rinex_body = ''; $$rinex_header = ''; } my $filegps = TimeClass->new->set_yrdoyhms_gps($fileyr, $filedoy, 0, 0, 0)->get_gps; # ## Early version of TRIG firmware tracked multiple channels for the same time/antenna/PRN ## This causes problems with reference link data jumping from channel to channel, and ## other issues. If the firmware is v433 or earlier, get rid of duplicate channels. # my $culled_data; if ($firmware_version eq 'v434+') { $culled_data = choose_all_tracks($data); # Do not cull duplicate channels } else { # older firmware # de-conflict multiple tracking of the same PRN in a different way for post-processing if (defined($ENV{versionMode}) && $ENV{versionMode} eq 'post') { $culled_data = choose_tracks($data); # Better with PP } else { $culled_data = choose_tracks2($data); # Simple, good with RT data } } my @times = sort keys %$culled_data; foreach my $time (@times) { if ($rnxver == 2.20) { $$rinex_body .= writeRinex220 ($culled_data->{$time}, $time, $options); } elsif ($rnxver == 2.11) { $$rinex_body .= writeRinex211 ($culled_data->{$time}, $time, $options); } } my ($tstart, $tend) = @times[0,-1]; if (exists ($options->{rinex_header})) { $$rinex_header = writeRinexHeader ($antid, $leoname, $tstart, $tend, $mission, {RNXVER => $rnxver}); return; # $$rinex_body is already filled up } my $outfh; if($options->{COMPRESS}){ $SIG{PIPE} = 'IGNORE'; # If there is a problem with this open, SIGPIPE can be thrown. open $outfh, "| rnx2crx - -s -e 100 | gzip - > $outfile" or die "can't fork: $!"; }else{ $outfile =~ s/_crx\.gz$/_rnx/i; open $outfh, "> $outfile"; } print {$outfh} writeRinexHeader ($antid, $leoname, $tstart, $tend, $mission, {RNXVER => $rnxver}) or die "can't write: $!"; print {$outfh} $$rinex_body or die "can't write: $!"; close ($outfh) or die "can't close: status=$?"; return; } #/**---------------------------------------------------------------------- # @sub choose_tracks # # Get rid of duplicate data for the same PRN and antenna. The TRIG # sometimes tracks the same PRN on two or more channels for the same # time. Choose the best track for a given time and PRN # based on continuity of channel or the length of track for that channel # Pass back a culled dataset. # # @parameter # @ $data Data hashref: $data->{gpstime}{fullprn}{primary_chan} # @return $culleddata Data hashref: $culleddata->{gpstime}{fullprn} # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub choose_tracks { my $data = shift; my %count; # $count{fullprn}{primary_chan} = N foreach my $time (keys %$data) { foreach my $prn (keys %{$$data{$time}}) { foreach my $chan (keys %{$$data{$time}{$prn}}) { $count{$prn}{$chan}++; } } } # Dataset to pass out my %culled_data; my $duplicate_count = 0; my $epoch_count = 0; my $prn_epoch_count = 0; # Store which channel used for each PRN and time. Try to maintain channel continuity my %chan_used; # $chan_used{fullprn}{gpstime} = channel_id foreach my $time (sort keys %$data) { $epoch_count++; foreach my $prn (sort keys %{$$data{$time}}) { $prn_epoch_count++; # Sort channels by overall count for each PRN my @chans = sort { $count{$prn}{$a} <=> $count{$prn}{$b} } keys %{$$data{$time}{$prn}}; $duplicate_count += (scalar(@chans)-1); # if @chans == 1, no duplicates for this time/prn # Use the same channel as last epoch, or if not known, use the channel most used for this PRN this dump my $last_chan = $chan_used{$prn}{$time-1}; my $chan_to_use; # If there is a defined channel from the last time step, and if that channel is available # this time step, then use it. if (defined($last_chan) && defined($$data{$time}{$prn}{$last_chan})) { $chan_to_use = $last_chan; # keep continuity # If there is a defined channel from the last time step, but that channel is *not* available # this time step, then insert a gap. Do not allow this track to jump from one channel to another } elsif (defined($last_chan)) { $chan_to_use = undef; # If there is no defined channel from the last time step, then use the channel with the most # data } else { $chan_to_use = $chans[-1]; } $chan_used{$prn}{$time} = $chan_to_use; $culled_data{$time}{$prn} = $$data{$time}{$prn}{$chan_to_use} // undef; # use highest count (last in sort) } } print "Epochs: $epoch_count epochs/prn: $prn_epoch_count duplicates: $duplicate_count\n"; return \%culled_data; } #/**---------------------------------------------------------------------- # @sub choose_all_tracks # # Starting with TRIG firmware version 4.3.4, there should be no # duplicate tracks. # # Just copy all data to the output structure. # # @parameter # @ $data Data hashref: $data->{gpstime}{fullprn}{primary_chan} # @return $culleddata Data hashref: $culleddata->{gpstime}{fullprn} # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub choose_all_tracks { my $data = shift; # Dataset to pass out my %culled_data; my $duplicate_count = 0; print "No data culling: v434+\n"; foreach my $time (sort keys %$data) { foreach my $prn (sort keys %{$$data{$time}}) { my @chans_tracked = keys %{$$data{$time}{$prn}}; $duplicate_count += (scalar(@chans_tracked)-1); # if @chans == 1, no duplicates for this time/prn my $chan_to_use = $chans_tracked[0]; $culled_data{$time}{$prn} = $$data{$time}{$prn}{$chan_to_use}; } } print "Warning: Found duplicate channels in v434 data\n" if ($duplicate_count > 0); return \%culled_data; } #/**---------------------------------------------------------------------- # @sub choose_tracks2 # # Get rid of duplicate data for the same PRN and antenna. The TRIG # sometimes tracks the same PRN on two or more channels for the same # time. Choose the track for a given time and PRN # based on sole use of the most prevalent channel. If this epoch does not # include the most prevalent channel, do not record it. Should get rid of # all phase jumps from switching channels, but at the cost of some data loss. # # @parameter # @ $data Data hashref: $data->{gpstime}{fullprn}{primary_chan} # @return $culleddata Data hashref: $culleddata->{gpstime}{fullprn} # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub choose_tracks2 { my $data = shift; my %count; # $count{fullprn}{primary_chan} = N foreach my $time (keys %$data) { foreach my $prn (keys %{$$data{$time}}) { foreach my $chan (keys %{$$data{$time}{$prn}}) { $count{$prn}{$chan}++; } } } # Dataset to pass out my %culled_data; my $duplicate_count = 0; my $epoch_count = 0; my $prn_epoch_count = 0; my $total_data_count = 0; foreach my $time (sort keys %$data) { $epoch_count++; foreach my $prn (sort keys %{$$data{$time}}) { $prn_epoch_count++; my @chans_tracked = keys %{$$data{$time}{$prn}}; $duplicate_count += (scalar(@chans_tracked)-1); # if @chans == 1, no duplicates for this time/prn # Only use most prevalent channel my @chans = sort { $count{$prn}{$a} <=> $count{$prn}{$b} } keys %{$count{$prn}}; # old_dedup my $chan_to_use = $chans[-1]; next if (!defined($chan_to_use)); # Skip if no preferred channel next if (!defined($$data{$time}{$prn}{$chan_to_use})); # skip if preferred channel not available # print "Adding $prn $time $chan_to_use\n"; $culled_data{$time}{$prn} = $$data{$time}{$prn}{$chan_to_use}; $total_data_count++; } } print "Epochs: $epoch_count epochs/prn: $prn_epoch_count duplicates: $duplicate_count total_data: $total_data_count\n"; return \%culled_data; } #/**---------------------------------------------------------------------- # @sub writeRinex # # Write one time step of RINEX data # # @parameter # @ $rec Record hashref by PRN with input data for one time step # @return $data The RINEX text of the record # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub writeRinex { my %rec = %{(shift)}; my $time = shift; my $data = ''; my @obs = qw(L1 L2 C1 P1 P2 LA SA S2); # From [IGSMAIL-3281]: RINEX Modifications for LEO Satellites (Vers. 2.20) # Werner Gurtner # Lou Estey # April 12, 2001 # L1, L2: L1, L2 phase, derived from P code # P1, P2: P1, P2 pseudorange from P code # C1: Pseudorange using C/A-Code on L1 # LA: Phase measurements on L1 derived from C/A code tracking # SA: Raw signal strength or SNR for LA # S2: Raw signal strengths or SNR values as given by the receiver for the L2 phase observation my $dataformat = "%14.3f"."%16.3f" x 4 ."\n"."%14.3f"."%16.3f" x (@obs - 6) ."\n"; my @prnlist = sort keys %rec; # Inlined from TimeClass for speed my @mkta = gmtime($time + $TimeClass::GPSSEC); my @epoch = (substr($mkta[5] + 1900, 2, 2), $mkta[4] + 1, $mkta[3], $mkta[2], $mkta[1], $mkta[0]); my $epochformat = " %0.2d"."%3d" x 4 ."%11.7f"." 0%3d"; if (@prnlist <= 12) { $epochformat .= "%3s" x @prnlist ."\n"; } else { # Allow for continuation line in case of more than 12 PRNs. $epochformat .= "%3s" x 12 ."\n"." " x 32 ."%3s" x (@prnlist - 12) ."\n"; } $data .= sprintf $epochformat, @epoch, scalar(@prnlist), @prnlist; foreach my $prn (@prnlist) { $data .= sprintf $dataformat, @{$rec{$prn}}{qw(p1phase p2phase capseudo p1pseudo p2pseudo caphase casnr p2snr)}; } return $data; } #/**---------------------------------------------------------------------- # @sub writeRinex211 # # Write one time step of RINEX data (RINEX 2.11 format) # # @parameter # @ $rec Record hashref by PRN with input data for one time step # @ $time Time of record in GPS seconds # @ $opt Options hashref: CONID => 'G|R' -- restrict to specified constellation # @return $data The RINEX text of the record # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub writeRinex211 { my %rec = %{(shift)}; my $time = shift; my $opt = shift // {}; my $conid = $$opt{CONID} // ''; my $data = ''; my @obs = qw(L1 L2 C1 C2 P2 S1 S2); # New for RINEX 2.11 my $dataformat1 = "%14.3f%16.3f%16.3f%16s%16.3f\n"."%14.3f"."%16.3f" x (@obs - 6) ."\n"; my $dataformat2 = "%14.3f%16.3f%16.3f%16.3f%16s\n"."%14.3f"."%16.3f" x (@obs - 6) ."\n"; my @prnlist = sort keys %rec; @prnlist = grep /$conid/, @prnlist if ($conid); # Inlined from TimeClass for speed my @mkta = gmtime($time + $TimeClass::GPSSEC); my @epoch = (substr($mkta[5] + 1900, 2, 2), $mkta[4] + 1, $mkta[3], $mkta[2], $mkta[1], $mkta[0]); my $epochformat = " %0.2d"."%3d" x 4 ."%11.7f"." 0%3d"; my $nleft = scalar(@prnlist); my $vals_this_line = ($nleft < 12) ? $nleft : 12; $epochformat .= "%3s" x $vals_this_line ."\n"; $nleft -= $vals_this_line; while (1) { last if ($nleft <= 0); $vals_this_line = ($nleft < 12) ? $nleft : 12; $epochformat .= " " x 32 ."%3s" x ($vals_this_line) ."\n"; $nleft -= $vals_this_line; } $data .= sprintf $epochformat, @epoch, scalar(@prnlist), @prnlist; foreach my $prn (@prnlist) { # Now write RINEX 2.11, using the SIGTYPE1, and SIGTYPE2 fields read in from the TRIG: # SIGTYPE1 = 0x03 = GPS L1/CA, # SIGTYPE2 = 0x30 = GPS L2P, 0x38 = GPS L2C (combined), 0x23 = GLONASS L2 # Jan: P1 should be blank # Add C2 field: For sigtype 0x38 (L2C), set C2 = RNG2, P2 = blank # For sigtype 0x30 (L2P), set C2 = blank, P2 = RNG2 # For sigtype 0x23 (GLONASS L2), C2 = RNG2, P2 = blank # Get rid of LA # { # Note: The DUMM key is not present in the hash, so results in a zero value (legal in RINEX for 'missing') no warnings 'uninitialized'; if ($rec{$prn}{SIGTYPE2} == 0x30) { # L2P $data .= sprintf $dataformat1, @{$rec{$prn}}{qw(PHS1 PHS2 RNG1 DUMM RNG2 SNR1 SNR2)}; } else { # L2C or GLONASS $data .= sprintf $dataformat2, @{$rec{$prn}}{qw(PHS1 PHS2 RNG1 RNG2 DUMM SNR1 SNR2)}; } } } return $data; } #/**---------------------------------------------------------------------- # @sub writeRinex220 # # Write one time step of RINEX data (RINEX 2.20 format) # # @parameter # @ $rec Record hashref by PRN with input data for one time step # @ $time Time of record in GPS seconds # @ $opt Options hashref: CONID => 'G|R' -- restrict to specified constellation # @return $data The RINEX text of the record # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub writeRinex220 { my %rec = %{(shift)}; my $time = shift; my $opt = shift // {}; my $conid = $$opt{CONID} // ''; my $data = ''; my @obs = qw(L1 L2 C1 P1 P2 LA SA S2); # From [IGSMAIL-3281]: RINEX Modifications for LEO Satellites (Vers. 2.20) # Werner Gurtner # Lou Estey # April 12, 2001 # L1, L2: L1, L2 phase, derived from P code # P1, P2: P1, P2 pseudorange from P code # C1: Pseudorange using C/A-Code on L1 # LA: Phase measurements on L1 derived from C/A code tracking # SA: Raw signal strength or SNR for LA # S2: Raw signal strengths or SNR values as given by the receiver for the L2 phase observation my $dataformat = "%14.3f"."%16.3f" x 4 ."\n"."%14.3f"."%16.3f" x (@obs - 6) ."\n"; my @prnlist = sort keys %rec; @prnlist = grep /$conid/, @prnlist if ($conid); # Inlined from TimeClass for speed my @mkta = gmtime($time + $TimeClass::GPSSEC); my @epoch = (substr($mkta[5] + 1900, 2, 2), $mkta[4] + 1, $mkta[3], $mkta[2], $mkta[1], $mkta[0]); my $epochformat = " %0.2d"."%3d" x 4 ."%11.7f"." 0%3d"; my $nleft = scalar(@prnlist); my $vals_this_line = ($nleft < 12) ? $nleft : 12; $epochformat .= "%3s" x $vals_this_line ."\n"; $nleft -= $vals_this_line; while (1) { last if ($nleft <= 0); $vals_this_line = ($nleft < 12) ? $nleft : 12; $epochformat .= " " x 32 ."%3s" x ($vals_this_line) ."\n"; $nleft -= $vals_this_line; } $data .= sprintf $epochformat, @epoch, scalar(@prnlist), @prnlist; foreach my $prn (@prnlist) { $data .= sprintf $dataformat, @{$rec{$prn}}{qw(PHS1 PHS2 RNG1 RNG1 RNG2 PHS1 SNR1 SNR2)}; } return $data; } #/**---------------------------------------------------------------------- # @sub writeRinexHeader # # Write the header for a RINEX file # # @parameter # @ $antid Antenna ID (integer) # @ $tstart Start time of data in GPS seconds # @ $tend End time of data in GPS seconds # @ $mission Mission name # @ $opt Optional options # @return $data The RINEX header text # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub writeRinexHeader { my $antid = shift; my $leoname = shift; my $tstart = shift; my $tend = shift; my $mission = $ENV{cdaac_mission} = shift; my $opt = shift // {}; my $data = ''; # Get program name and version. my $prgname = $0; my $longheader = `strings $prgname | grep -m 1 '\$Id:'` // ''; ($prgname = basename($prgname)) =~ s/\.pl//; $prgname = $mission if (length($prgname) >= 13); my ($f, $version, $date, $time, $user) = ($longheader =~ m{\$Id:\s+([\w\.]+)\s+(\d+)\s+([\d\-]+)\s+([\d\:Z]+)\s+(\w+)}); if (defined($version)) { $prgname .= " v$version"; } die "PGM name ($prgname) and version\# too long" if (length($prgname) > 20); # Defaults for RINEX header if no other info. my $runby = "COSMIC CDAAC"; my $marktype = "SPACEBORNE"; my $rectype = ""; my $anttype = ""; my $recnum = ""; my $antnum = ""; my $recver = ""; my $rnxver = $opt->{RNXVER} // 2.2; my @obs = $rnxver == 2.2 ? qw(L1 L2 C1 P1 P2 LA SA S2) : qw(L1 L2 C1 C2 P2 S1 S2); my @wavefac = (1,1); my @c_mass = (0,0,0); my @antdel = (0,0,0); my @antbor = (0,0,0); # Read Antenna parameters from the Config.pm file. my $config = do "$mission/Config.pm"; die "Cannot find Config.pm file for $mission" if (!defined($config)); my %antennas = %{$config->{'satellite'}{'antennas'}}; my ($antname) = grep { $antid == $antennas{$_}{'id'} } keys %antennas; my $nleo = ($leoname =~ /^CHAM/) ? 0 : ($leoname =~ /^CNFS/) ? 0 : ($leoname =~ /^C2E/) ? (substr $leoname, 3)-1 # COSMIC-2: 1-6 convert to index: 0-5 : ($leoname =~ /^C/) ? (substr $leoname, 1)-1 # COSMIC: 1-6, convert to index: 0-5 : ($leoname =~ /^GRC/) ? (substr $leoname, 3)-1 # GRACE 1-2 (GRACE-A or GRACE-B), convert to index, 0-1 : 0; @antdel = ref($antennas{$antname}{'offset'}->[0]) eq 'ARRAY' ? @{$antennas{$antname}{'offset'}->[$nleo]} : @{$antennas{$antname}{'offset'}}; @antbor = ref($antennas{$antname}{'boresight'}->[0]) eq 'ARRAY' ? @{$antennas{$antname}{'boresight'}->[$nleo]} : @{$antennas{$antname}{'boresight'}}; # Determine MARKER NAME. Special logic for cosmic_GOX_reader # C006 L26 (POD2) MARKER NAME my $bernid = $config->{'satellite'}{'id_num'}[$nleo]; my $markname = sprintf "$leoname L%02d ($antname)", $bernid-900; $rectype = $config->{'satellite'}{'rectype'}; $anttype = $config->{'satellite'}{'antennas'}{$antname}{'fullname'}; # OBSERVER / AGENCY my %obsagency = (gpsmet => ["GPS/MET", "UCAR"], champ => ["CHAMP", "GFZ"], grace => ["GRACE", "GFZ"], tsx => ["TSRX", "GFZ"], tdx => ["TDMX", "GFZ"], metopa => ["MTPA", "EUMETSAT"], metopb => ["MTPB", "EUMETSAT"], metopc => ["MTPC", "EUMETSAT"], sacc => ["SAC-C", "CONAE/NASA"], cosmic => ["FORMOSAT3/COSMIC", "NSPO/UCAR"], coseq => ["FORMOSAT7/COSMIC-2E", "NSPO/UCAR"], cos2e => ["FORMOSAT7/COSMIC-2E", "NSPO/UCAR"], kompsat5 => ["KOMPSAT5", "KASI/KARI"], cnofs => ["C/NOFS", "US AIR FORCE"], ); my $proto = $mission; die "Proto mission $proto not recognized" unless exists $obsagency{$proto}; my $observer = $obsagency{$proto}->[0]; my $agency = $obsagency{$proto}->[1]; # Get times of first and last observation. my @firstepoch = TimeClass->new->set_gps($tstart)->get_ymdhms_gps; my @lastepoch = TimeClass->new->set_gps($tend)->get_ymdhms_gps; # Get date and time of this processing. my $timestamp = TimeClass->new->now->get_bernestamp_gps; # Write RINEX header. $data .= sprintf "%9.2f"." " x 11 ."%-20s" x 2 ."%1s\n", $rnxver, "OBSERVATION DATA", "M (MIXED)", "RINEX VERSION / TYPE"; $data .= sprintf "%-20s" x 3 ."%1s\n", $prgname, $runby, $timestamp, "PGM / RUN BY / DATE"; $data .= sprintf "%-60s%1s\n", $markname, "MARKER NAME"; $data .= sprintf "%-60s%1s\n", $marktype, "MARKER TYPE"; $data .= sprintf "%-20s%-40s%1s\n", $observer, $agency, "OBSERVER / AGENCY"; $data .= sprintf "%-20s" x 3 ."%1s\n", $recnum, $rectype, $recver, "REC \# / TYPE / VERS"; $data .= sprintf "%-20s%-40s%1s\n", $antnum, $anttype, "ANT \# / TYPE"; $data .= sprintf "%14.4f" x 3 ." " x 18 ."%1s\n", @antdel, "ANTENNA: DELTA X/Y/Z"; $data .= sprintf "%14.4f" x 3 ." " x 18 ."%1s\n", @antbor, "ANTENNA: B.SIGHT XYZ"; $data .= sprintf "%14.4f" x 3 ." " x 18 ."%1s\n", @c_mass, "CENTER OF MASS: XYZ"; $data .= sprintf "%6d" x 2 ." " x 48 ."%1s\n", @wavefac, "WAVELENGTH FACT L1/2"; $data .= sprintf "%6d" . "%6s" x @obs . " " x (54 - 6*@obs) ."%1s\n", scalar(@obs), @obs, "\# / TYPES OF OBSERV"; $data .= sprintf "%6d" x 5 ."%13.7f" ." " x 5 ."%3s" ." " x 9 . "%1s\n", @firstepoch, "GPS", "TIME OF FIRST OBS"; $data .= sprintf "%6d" x 5 ."%13.7f" ." " x 5 ."%3s" ." " x 9 . "%1s\n", @lastepoch, "GPS", "TIME OF LAST OBS"; $data .= sprintf " " x 60 ."%1s\n", "END OF HEADER"; return $data; } #/**---------------------------------------------------------------------- # @sub sortOpnGps # # Sort high rate data into the correct order after BJfmtl_cosmic # is called. Add headers and trailers. # # # @parameter $infile Input file in opnGps format, but out of order and # @ with no headers/trailers. # @ $outfile Output sorted file. # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub sortOpnGps { my $infile = shift; my $outfile = shift; my $opt = shift || {}; my $version = defined($opt->{VERSION}) ? $opt->{VERSION} : 1; # default to open loop format my $accept_odd_rates = $opt->{ODD_RATES}; open (IN, $infile) || die "Cannot open $infile"; my %data = (); # ## Read in data and stash in a hash by PRN and time # my $hrlen = ($version == 1) ? 36 # open loop : ($version == 2) ? 24 # closed loop : ($version == 3) ? 48 # l2c open loop : 0; die "bad version number: $version" if ($hrlen == 0); my $hdrlen = ($version == 3) ? 28 # l2c open loop : 20; while (1) { read (IN, my $hdr = '', $hdrlen) or last; my ($rate, $prn, $t) = unpack ('SSL', $hdr); my $nok = $accept_odd_rates ? ($prn < 1 || $prn > 32 || $t < 1) : ($rate < 46 || $rate > 52 || $prn < 1 || $prn > 32 || $t < 1); die "Bad rate, prn or time: $rate, $prn, $t" if ($nok); read (IN, my $hr = '', $rate * $hrlen) or last; $data{$prn}{$t} = $hdr . $hr; } close (IN); # ## Print data to output file sorted by PRN and time. ## Keep track of byte offsets of beginning of each PRN section # my @offsets = (-1) x 33; my $byteCount = 0; open (OUT, ">$outfile") || die "Cannot open $outfile for writing"; foreach my $prn (1..32) { $offsets[$prn] = (exists $data{$prn}) ? $byteCount : 0; foreach my $time (sort { $a <=> $b } keys %{$data{$prn}}) { print OUT $data{$prn}{$time}; $byteCount += length($data{$prn}{$time}); } } # Append PRN index block to end of opnGps file print OUT pack ("l" x 32, @offsets[1..32]); # Append Trailer block describing version and column layout # Version 1, high rate data in "fddSSdf" format, one second data in "SSLdCCCC" format # (see perldoc -f pack) if ($version == 1) { print OUT pack ("Ca31a32", 1, "fddSSdf", "SSLdCCCC"); # 64 bytes of trailer } elsif ($version == 2) { # closed loop data (no CAmodel or dfaz information) print OUT pack ("Ca31a32", 2, "fddSS", "SSLdCCCC"); # 64 bytes of trailer } elsif ($version == 3) { # L2C open loop data (P2 range, P2 phase model and P2 delta phase included) print OUT pack ("Ca31a32", 3, "fddSSdfdf", "SSLdCCCCd"); # 64 bytes of trailer } close OUT; return; } #/**---------------------------------------------------------------------- # @sub sort_opnGns # # Sort high rate data into the correct order after BJfmtl_cosmic # is called. Add headers and trailers. # # # @parameter $infile Input file in opnGns format, but out of order and # @ with no headers/trailers. Can also be a perl scalar containing data. # @ $outfile Output sorted file. # @ $lrdesc A hash by GNSS type, referencing a list of low rate data descriptors, up to # @ 16 unsigned shorts. GNSS types are 'G' for GPS, 'R' for GLONASS, etc. # @ $hrdesc A hash by GNSS type, referencing a list of high rate data descriptors, up to # @ 16 unsigned shorts with a perl 'pack' data type for each one. # @ Data descriptors look like: # @ 3 bits obs type: 0=C, 1=L, 2=D, 3=S, 4-7 reserved # @ 3 bits frequency band: 0-7 add +1 to interpret # @ 5 bits tracking attribute: 0, 1=A to 26=Z # @ 5 bits metadata: 0, 1=A to 26=Z. This is only defined for CDAAC. 'M' = phase model # @ $opt Options hashref. { HROFFSET => offset_in_seconds } default = -0.5 # @see RINEX 3.02 specification, section 5: http://igscb.jpl.nasa.gov/igscb/data/format/rinex302.pdf # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub sort_opnGns { my $infile = shift; my $outfile = shift; my $lrdesc = shift; my $hrdesc = shift; my $opt = shift || {}; my $lr_hdr_len = 12; my @conids = sort keys %$lrdesc; # ## Compute the length of the high rate data section for each constellation ID # my %hrlen; foreach my $conid (@conids) { my @pack_desc = map { substr ($_, 2, 1) } @{$hrdesc->{$conid}}; # the perl 'pack' descriptors ('d', 'S', etc) $hrlen{$conid} = pdl(map { length(pack($_, 0)) } @pack_desc)->sum; } my $infh; if (length($infile) > 512) { # in-memory file open $infh, '<', \$infile or die "Cannot open in-memory file $infile"; } else { # Regular file open $infh, '<', $infile or die "Cannot open input file $infile"; } my %data; while (1) { read ($infh, my $hdr = '', $lr_hdr_len) or last; my ($tpkt, $status, $antid, $conid, $svn, $prn, $rate) = unpack ('(LCCASCS)>', $hdr); # Read one double per data descriptor for this constellation ID read ($infh, my $lrvals = '', scalar (@{$lrdesc->{$conid}}) * 8); my $hrlen = 3 + $hrlen{$conid}; # 3 byte time offset added read ($infh, my $hr = '', $rate * $hrlen) or last; $data{$conid}{$prn}{$tpkt} = $hdr . $lrvals . $hr; } close ($infh); # ## Print data to output file sorted by Constellation, PRN and time. ## Keep track of byte offsets of beginning of each PRN section # @conids = sort keys %data; my %offsets; open my $ofh, '>', $outfile or die "Cannot open $outfile for writing"; my $byteCount = 0; foreach my $conid (@conids) { my @offsets = (-1) x 33; foreach my $prn (1..32) { $offsets{$conid}[$prn] = (exists $data{$conid}{$prn}) ? $byteCount : -1; foreach my $time (sort { $a <=> $b } keys %{$data{$conid}{$prn}}) { print {$ofh} $data{$conid}{$prn}{$time}; $byteCount += length($data{$conid}{$prn}{$time}); } } } my $hroffset = $opt->{HROFFSET} // -0.5; foreach my $conid (@conids) { # Append low- and high-rate data descriptors to file my @lrdesc = @{$lrdesc->{$conid}}; my @hrdesc = @{$hrdesc->{$conid}}; print {$ofh} join ('', (@lrdesc, map { pack 'S>', 0 } (0) x (16 - scalar(@lrdesc)))); print {$ofh} join ('', (@hrdesc, map { pack 'S>C', (0,0) } (0) x (16 - scalar(@hrdesc)))); print {$ofh} pack 'l>', $hroffset * 1e9; # Set high rate data offset to -0.5e9 nanoseconds for JPL RO receiver # Append PRN index block to end of opnGns file print {$ofh} pack ("q>" x 32, @{$offsets{$conid}}[1..32]); } # Append Trailer block with version ID and constellation IDs my $version = 1; print {$ofh} pack ('C', $version) . pack ("A" x 8, (@conids, ((' ') x (8 - scalar(@conids))))); close $ofh; return; } #/**---------------------------------------------------------------------- # @sub sort_opnGns_v2 # # Sort high rate data into the correct order after Cosmic2Parser is called. # is called. Add headers and trailers. Write a version 2 opnGns file. # # # @parameter $data Input data consisting of a hashref by time and satid with # @ records in opnGns format # @ Can *also* be a hashref sorted like %data below. In this case # @ skip the first half of this routine # @ $outfile Output sorted file. # @ $opt Options hashref. { HROFFSET => offset_in_seconds } default = -0.5 # @see RINEX 3.02 specification, section 5: http://igscb.jpl.nasa.gov/igscb/data/format/rinex302.pdf # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub sort_opnGns_v2 { my $infile = shift; my $outfile = shift; my $opt = shift || {}; my $lr_hdr_len = 14; my $infh; my %data; if (ref($infile) eq 'HASH') { %data = %$infile; # treat 'infile' as the %data hash } else { if (length($infile) > 512) { # in-memory file open $infh, '<', \$infile or die "Cannot open in-memory file $infile"; } else { # Regular file open $infh, '<', $infile or die "Cannot open input file $infile"; } while (1) { read ($infh, my $hdr = '', $lr_hdr_len) or last; my ($tpkt, $status, $antid, $conid, $svn, $prn, $rate, $nlr, $nhr) = unpack ('(LCCASCSCC)>', $hdr); read ($infh, my $lrhdr = '', $nlr*2); # Variable length LR headers read ($infh, my $hrhdr = '', $nhr*3) or last; # Variable length HR headers # Read one double per data descriptor for this constellation ID read ($infh, my $lrvals = '', $nlr * 8); # 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 $hrlen = 3 + pdl(map { length(pack($_, 0)) } @pack_desc)->sum; # 3 byte time offset + sum of high rate data lengths read ($infh, my $hr = '', $rate * $hrlen) or last; $data{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr; } close ($infh); } # ## Print data to output file sorted by Constellation, PRN and time. ## Keep track of byte offsets of beginning of each PRN section # my @conids = sort keys %data; my %offsets; open my $ofh, '>', $outfile or die "Cannot open $outfile for writing"; my $byteCount = 0; foreach my $conid (@conids) { #my @offsets = (-1) x 33; # original my @offsets = (-1) x 37; # Jun modified, to include E36 #foreach my $prn (1..32) { # original foreach my $prn (1..36) { # Jun modified, to include E36 $offsets{$conid}[$prn] = (exists $data{$conid}{$prn}) ? $byteCount : -1; foreach my $time (sort { $a <=> $b } keys %{$data{$conid}{$prn}}) { print {$ofh} $data{$conid}{$prn}{$time}; $byteCount += length($data{$conid}{$prn}{$time}); } } } my $hroffset = $opt->{HROFFSET} // -0.5; #my $hroffset = $opt->{HROFFSET} // -1e-9; #Jun modified, to be consistent with UCAR #my $hroffset = $opt->{HROFFSET} // 0.0; #Jun modified, to be consistent with UCAR, for Galileo foreach my $conid (@conids) { print {$ofh} pack 'l>', $hroffset * 1e9; # Set high rate data offset to -0.5e9 nanoseconds for JPL RO receiver # Append PRN index block to end of opnGns file #print {$ofh} pack ("q>" x 32, @{$offsets{$conid}}[1..32]); #original print {$ofh} pack ("q>" x 36, @{$offsets{$conid}}[1..36]); # Jun modified, to include E36 } # Append Trailer block with version ID and constellation IDs my $version = 2; print {$ofh} pack ('C', $version) . pack ("A" x 8, (@conids, ((' ') x (8 - scalar(@conids))))); close $ofh; return; } #/**---------------------------------------------------------------------- # @sub sortLrObs # # Sort Low rate data ASCII files, as output by the GRACE parsing software # and the BJfmtl* family of C codes from JPL. # # Reads/writes from disk files or perl strings. # # @parameter $infile Input file in opnGps format, but out of order and # @ with no headers/trailers. # @ $outfile Output sorted file. # @return none # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub sortLrObs { my $infile = shift; my $outfile = shift; my %data = (); # ## Read in data and stash in a hash by PRN and time # my $infh; if (ref($infile) eq 'SCALAR') { open $infh, '<', $infile or die "Cannot open in-memory input scalar"; } else { open $infh, '<', $infile or die "Cannot $infile for reading"; } while (my $line = <$infh>) { my ($time, $prn) = ($line =~ /^(\d+) OBSD prn(\d+)/); $data{$time}{$prn} = $line; } close $infh; # ## Print data to output file sorted by time. # my $outfh; if (ref($outfile) eq 'SCALAR') { open $outfh, '>', $outfile or die "Cannot open in-memory output scalar"; } else { open $outfh, '>', $outfile or die "Cannot open $outfile for writing"; } foreach my $time (sort { $a <=> $b } keys %data) { foreach my $prn (sort { $a <=> $b } keys %{$data{$time}}) { print {$outfh} $data{$time}{$prn}; } } close $outfh; return; } #/**---------------------------------------------------------------------- # @sub pod2opn # # Convert an ASCII JPL-format file into a one-second data rate opnGps file. # # @parameter $infile -- JPL ASCII format file to convert # @output $outfile -- opnGps format output file # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub pod2opn { my $infile = shift; my $outfile = shift; open my $infh, '<', $infile or die "Cannot open $infile for reading"; open my $outfh, '>', $outfile or die "Cannot open $outfile for writing"; my $dt = 0; my $rate = 1; my $antin = 0; RECORD: while (<$infh>) { my ($tsec, $obsd, $prn, $antidxm, $obstypem, $interval, $pcktsize, $ca, $cachan, $casnr, $caphase, $capseudo, $p2, $p2chan, $p2snr, $p2phase, $p2pseudo, $p1, $p1chan, $p1snr, $p1phase, $p1pseudo) = split; next unless (defined ($obsd) && $obsd =~ /OBSD/); # skip nav data next unless ($rate == hex($interval)); # skip data of incorrect rate $prn = substr($prn,3,2); # get rid of 'prn' in eg 'prn08' # create opnGps header (data run from -0.5 to +0.48 if 50Hz my $hrhdr = pack ("SSIdCCCC", # see opnGps.h: typedef struct opnGpsHdr $rate, $prn, $tsec, $capseudo, 0, $antin, 0, 0); # $bjdlength, $trkstatus, $flywheel set to zero. # create opnGps one-second body record (convert phases to meters) my $hrrec = pack ("fddSS", $dt, $caphase * $LAMB1, $p2phase * $LAMB2, $casnr*10, $p2snr*10); # Append record onto output file print {$outfh} $hrhdr, $hrrec; } close $infh; close $outfh; return; } #/**---------------------------------------------------------------------- # @sub bin2dec # # Convert an arbitrary length (<= 32 bits long) binary number to a normal # perl value. # Input: '100110' # Output: 38 # # @parameter $bits ASCII string of 0s and 1s <= 32 characters in length # @return decimal value of $bits # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub bin2dec { return unpack("N", pack("B32", substr("0" x 32 . shift, -32))); } #/**---------------------------------------------------------------------- # @sub write_all_opnScn # # Subroutine to write out COSMIC-2 high rate scintillation data to an # opnScn file. Takes care of all the different versions of scintillation # data packets from the TRIG. # # @parameter # @ $outfiles -- hashref to scintillation data from TRIG # @ $date -- Date to put on output file # @ $ndump -- Dump number # @ $leoid -- LEO id # @ $mission -- Mission name # @return Nothing (file written) # @exception Will toss an exception in case of a real problem. # ----------------------------------------------------------------------*/ sub write_all_opnScn { my ($outfiles, $date, $ndump, $leoid, $mission,$opnScn_output,$hr,$min,$groundstation) = @_; #print "$opnScn_output,$hr,$min,$groundstation\n"; #exit; my @scin432 = sort keys %{$outfiles->{'SCNThrat432'}}; # v4.3.2 and older my @scin433 = sort keys %{$outfiles->{'SCNThrat433'}}; # v4.3.3 GLONASS data my @scin434 = sort keys %{$outfiles->{'SCNThrat434'}}; # v4.3.4 GPS data (needs data from @scin432 as well!) # If any v434 obstype 7,8 packets, treat as version 4.3.4 if (@scin434) { print "Processing scintillation v434 packets\n"; my @ants = pdl(@scin432, @scin433, @scin434)->uniq->qsort->list; foreach my $ant (@ants) { #my $outpath = sprintf ("opnScn_$date.%03d.%02d.%02d_bin", $leoid, $ndump, $ant); # Original code # Jun modified 20251208 ---------------------bg my $outpath = sprintf ("opnScn_$date.%02d.%02d.%03d.%02d.${groundstation}_bin", $hr,$min,$leoid,$ant); $outpath="$opnScn_output/$outpath"; print "outpath=$outpath\n"; # Jun modified 20251208 ---------------------ed # Combine v4.3.2, v4.3.3 and v4.3.4 scintillation packets into one opnScn file for this antenna write_opnScn_v434($outfiles, $ant, $outpath); } } elsif (@scin433) { print "Processing scintillation v433 packets\n"; my @ants = pdl(@scin432, @scin433)->uniq->qsort->list; foreach my $ant (@ants) { #my $outpath = sprintf ("opnScn_$date.%03d.%02d.%02d_bin", $leoid, $ndump, $ant); # Original code # Jun modified 20251208 ---------------------bg my $outpath = sprintf ("opnScn_$date.%02d.%02d.%03d.%02d.${groundstation}_bin", $hr,$min,$leoid,$ant); $outpath="$opnScn_output/$outpath"; print "outpath=$outpath\n"; # Jun modified 20251208 ---------------------ed # Combine v4.3.2 and v4.3.3 scintillation packets into one opnScn file for this antenna write_opnScn_v433($outfiles, $ant, $outpath); } } else { print "Processing scintillation v432 packets\n"; foreach my $ant (@scin432) { #my $outpath = sprintf ("opnScn_$date.%03d.%02d.%02d_bin", $leoid, $ndump, $ant); # Original code # Jun modified 20251208 ---------------------bg my $outpath = sprintf ("opnScn_$date.%02d.%02d.%03d.%02d.${groundstation}_bin", $hr,$min,$leoid,$ant); $outpath="$opnScn_output/$outpath"; print "outpath=$outpath\n"; # Jun modified 20251208 ---------------------ed # Process only v4.3.2 scintillation packets into one opnScn file for this antenna write_opnScn_v432($outfiles, $ant, $outpath); } } } #/**---------------------------------------------------------------------- # @sub write_opnScn_v432 # # Write an opnGns format file (file type = opnScn). This version is used # to process scintillation data under firmware revision 4.3.2. # # This version offers GPS and GLONASS data in 10-second, 50 or 100Hz obstype 1,2,3,4 packets. # # These data are not phase-connected! # # @parameter $parser -- Root of the TRIG data parser, containing # @ scintillation data in the SCNThrat, SCNThratII, and # @ SCNThratIII sections # @ $ant -- antenna ID # @ $outpath: Output file name # @output File $outpath written # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub write_opnScn_v432 { my $outfiles = shift; my $ant = shift; my $outpath = shift; my $v432 = $outfiles->{'SCNThrat432'}{$ant}; # obstype 1-4 GPS data # holder for GLONASS and GPS data in the Level1aTools::sort_opnGns_v2 format: # $alldata{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr my %alldata; # GLONASS data are reconstructed solely from the obstype 5,6 data # which are one-second packets with 100Hz I and Q data. # # GPS data are reconstructed from the obstype 7,8 data which # are 1 second phase packets containing 50Hz data, along with the # older SNRs from obstype 1,2 packets. These packets are 10 seconds # long and contain 50Hz SNRs. What a tangled web we weave! # 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 [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 [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 [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 ); my $status = 0x0; # closed loop, no flags # Determine which SATELLIT. file to use (the latest one) #my $Ivs = 'I14'; # hard wired! #my $satinfo = SatInfo->new("$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; foreach my $fullprn (sort { $a cmp $b} keys %$v432) { my ($conid, $prn) = unpack "a a2", $fullprn; $prn += 0; if ($prn == 0) { warn "PRN equals zero in SCNThrat data (very odd...)" if ($prn == 0); next; } # Time fraction fields my $rate = ($conid eq 'R') ? 100 : 50; my $hrt = sequence($rate)/$rate; my $tf = long(($hrt/1e-7)->rint); my $mostsig = byte(($tf & 0x00ff0000) >> 16); my $midsig = byte(($tf & 0x0000ff00) >> 8); my $lowsig = byte ($tf & 0x000000ff); my @times = sort { $a <=> $b} keys %{$v432->{$fullprn}}; my ($svn, $freq_slot, $f1, $f2) = $satinfo->find_svn($fullprn, $times[0]); my $needed_elems = (10*$rate); # 500 or 1000 TIME: foreach my $t10 (@times) { # ## Check to make sure required data for this time step are there # next TIME if ( !defined ($v432->{$fullprn}{$t10}{'1'}[0]) || !defined ($v432->{$fullprn}{$t10}{'2'}[0]) || !defined ($v432->{$fullprn}{$t10}{'3'}[0]) || !defined ($v432->{$fullprn}{$t10}{'4'}[0]) ); if ( ($v432->{$fullprn}{$t10}{'1'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'2'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'3'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'4'}[0]->nelem != $needed_elems) ) { next TIME; } # Now the problem is that we have 10 second epochs from the TRIG, but we need 1 second epochs for the opnGns format for (my $i=0;$i<10;$i++) { my $sidx = $i*$rate; my $eidx = ($i*$rate)+$rate-1; my $time = $t10 + $i; my $l1s = $v432->{$fullprn}{$t10}{'1'}[0]->mslice([$sidx, $eidx]); my $l2s = $v432->{$fullprn}{$t10}{'2'}[0]->mslice([$sidx, $eidx]); my $l1f = $v432->{$fullprn}{$t10}{'3'}[0]->mslice([$sidx, $eidx]); my $l2f = $v432->{$fullprn}{$t10}{'4'}[0]->mslice([$sidx, $eidx]); # Now pack all this together for one second my $rec = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $rate); # Determine GPS L2C vs. L2P my $signal_type2 = $v432->{$fullprn}{$t10}{4}[1]; # L2C vs. L2P for obstype = 4 (L2 phase) # Format Epoch Data Record (v2) Use @hrdesc defined # at a package level for speed Add low rate record count, high # rate record count, then high rate descriptors if ($conid eq 'G' && ($signal_type2 == 0x30)) { # L2P (0x30) instead of L2CL (0x2C) $rec .= pack ('CC', 0, scalar(@hrdesc_l2p_gps)) . join ('', @hrdesc_l2p_gps); } elsif ($conid eq 'R') { $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_glonass)) . join ('', @hrdesc_l2c_glonass); } else { # default to GPS L2C $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_gps)) . join ('', @hrdesc_l2c_gps); } my $hrmat = cat (double($mostsig), double($midsig), double($lowsig), $l1f, $l2f, $l1s, $l2s); my $hrrec = pack ("(CCCddSS)>" x $rate, $hrmat->xchg(0,1)->flat->list); $alldata{$conid}{$prn}{$time} = $rec.$hrrec; } } } Level1aTools::sort_opnGns_v2 (\%alldata, $outpath, {HROFFSET => 0}); print "FILE CREATION: ", Dplib::simplefn($outpath), "\n"; return; } #/**---------------------------------------------------------------------- # @sub write_opnScn_v433 # # Write an opnGns format file (file type = opnScn). This version is used # to process scintillation data under firmware revision 4.3.3. # # This version offers GPS data in 10-second, 50Hz obstype 1,2,3,4 packets # and GLONASS as 1 second, 100Hz obstype 5,6 packets. # # @parameter $outfiles -- Root of the TRIG data parser, containing # @ scintillation data in the SCNThrat, SCNThratII, and # @ SCNThratIII sections # @ $ant -- antenna ID # @ $outpath: Output file name # @output File $outpath written # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub write_opnScn_v433 { my $outfiles = shift; my $ant = shift; my $outpath = shift; my $v432 = $outfiles->{'SCNThrat432'}{$ant}; # obstype 1-4 GPS data my $v433 = $outfiles->{'SCNThrat433'}{$ant}; # obstype 5,6 GLONASS data # holder for GLONASS and GPS data in the Level1aTools::sort_opnGns_v2 format: # $alldata{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr my %alldata; # GLONASS data are reconstructed solely from the obstype 5,6 data # which are one-second packets with 100Hz I and Q data. # # GPS data are reconstructed from the obstype 7,8 data which # are 1 second phase packets containing 50Hz data, along with the # older SNRs from obstype 1,2 packets. These packets are 10 seconds # long and contain 50Hz SNRs. What a tangled web we weave! # 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 [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 [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 [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 ); # ## First process GLONASS data into the intermediate hash format required by ## Level1aTools::sort_opnGns_v2: ## $data{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr # my $status = 0x0; # closed loop, no flags # Determine which SATELLIT. file to use (the latest one) #my $Ivs = 'I14'; #my $satinfo = SatInfo->new("SATELLIT.$Ivs.bern")->parse; # Jun modified, use new file provided by Yong my $Ivs = 'I20'; # Hard-wired to IGS 14 coordinates my $satinfo = SatInfo->new("SATELLIT_$Ivs.SAT")->parse; foreach my $fullprn (sort { $a cmp $b} keys %$v433) { my ($conid, $prn) = unpack "a a2", $fullprn; $prn += 0; my @times = sort { $a <=> $b} keys %{$v433->{$fullprn}}; my ($svn, $freq_slot, $f1, $f2) = $satinfo->find_svn($fullprn, $times[0]); TIME: foreach my $time (@times) { if (!exists($v433->{$fullprn}{$time}{'L1'}) || !exists($v433->{$fullprn}{$time}{'L2'})) { print "Both L1 and L2 data not present for $fullprn at $time, skipping\n"; next TIME; } my ($dt1, $hrfaz1, $mdl1, $snr1, $signal_type1) = eval { @{$v433->{$fullprn}{$time}{'L1'}} }; die "write_opnScnII: Error reading L1 data for $time, $fullprn, quitting" if ($@); my ($dt2, $hrfaz2, $mdl2, $snr2, $signal_type2) = eval { @{$v433->{$fullprn}{$time}{'L2'}} }; die "write_opnScnII: Error reading L2 data for $time, $fullprn, quitting" if ($@); my $rate = $dt1->nelem; if ($dt1->at(0) != $dt2->at(0)) { die "L1 and L2 times do not match: PRN$fullprn T$time"; } my $tf = long(($dt1/1e-7)->rint); my $mostsig = byte(($tf & 0x00ff0000) >> 16); my $midsig = byte(($tf & 0x0000ff00) >> 8); my $lowsig = byte ($tf & 0x000000ff); # Now pack all this together for one second my $rec = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $rate); # Format Epoch Data Record (v2) Use @hrdesc defined # at a package level for speed Add low rate record count, high # rate record count, then high rate descriptors $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_glonass)) . join ('', @hrdesc_l2c_glonass); my $hrmat = cat (double($mostsig), double($midsig), double($lowsig), $hrfaz1, $hrfaz2, $snr1, $snr2); my $hrrec = pack ("(CCCddSS)>" x $rate, $hrmat->xchg(0,1)->flat->list); $alldata{'R'}{$prn}{$time} = $rec.$hrrec; } } # ## Now we have the GLONASS data (obstype 5,6) taken care of, let's look at the GPS data ## in the two packet types at two different cadences # foreach my $fullprn (sort { $a cmp $b} keys %$v432) { my ($conid, $prn) = unpack "a a2", $fullprn; $prn += 0; # Canonicalize if ($prn == 0) { warn "PRN equals zero in SCNThrat data (very odd...)" if ($prn == 0); next; } # Time fraction fields my $rate = ($conid eq 'R') ? 100 : 50; my $hrt = sequence($rate)/$rate; my $tf = long(($hrt/1e-7)->rint); my $mostsig = byte(($tf & 0x00ff0000) >> 16); my $midsig = byte(($tf & 0x0000ff00) >> 8); my $lowsig = byte ($tf & 0x000000ff); my @times = sort { $a <=> $b} keys %{$v432->{$fullprn}}; my ($svn, $freq_slot, $f1, $f2) = $satinfo->find_svn($fullprn, $times[0]); my $needed_elems = (10*$rate); # 500 or 1000 TIME: foreach my $t10 (@times) { # ## Check to make sure required data for this time step are there # next TIME if ( !defined ($v432->{$fullprn}{$t10}{'1'}[0]) || !defined ($v432->{$fullprn}{$t10}{'2'}[0]) || !defined ($v432->{$fullprn}{$t10}{'3'}[0]) || !defined ($v432->{$fullprn}{$t10}{'4'}[0]) ); if ( ($v432->{$fullprn}{$t10}{'1'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'2'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'3'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'4'}[0]->nelem != $needed_elems) ) { next TIME; } # Now the problem is that we have 10 second epochs from the TRIG, but we need 1 second epochs for the opnGns format for (my $i=0;$i<10;$i++) { my $sidx = $i*$rate; my $eidx = ($i*$rate)+$rate-1; my $time = $t10 + $i; my $l1s = $v432->{$fullprn}{$t10}{'1'}[0]->mslice([$sidx, $eidx]); my $l2s = $v432->{$fullprn}{$t10}{'2'}[0]->mslice([$sidx, $eidx]); my $l1f = $v432->{$fullprn}{$t10}{'3'}[0]->mslice([$sidx, $eidx]); my $l2f = $v432->{$fullprn}{$t10}{'4'}[0]->mslice([$sidx, $eidx]); # Now pack all this together for one second my $rec = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $rate); # Determine GPS L2C vs. L2P my $signal_type2 = $v432->{$fullprn}{$t10}{4}[1]; # L2C vs. L2P for obstype = 4 (L2 phase) # Format Epoch Data Record (v2) Use @hrdesc defined # at a package level for speed Add low rate record count, high # rate record count, then high rate descriptors if ($signal_type2 == 0x30) { # L2P (0x30) instead of L2CL (0x2C) $rec .= pack ('CC', 0, scalar(@hrdesc_l2p_gps)) . join ('', @hrdesc_l2p_gps); } else { # default to GPS L2C $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_gps)) . join ('', @hrdesc_l2c_gps); } my $hrmat = cat (double($mostsig), double($midsig), double($lowsig), $l1f, $l2f, $l1s, $l2s); my $hrrec = pack ("(CCCddSS)>" x $rate, $hrmat->xchg(0,1)->flat->list); $alldata{'G'}{$prn}{$time} = $rec.$hrrec; } } } Level1aTools::sort_opnGns_v2 (\%alldata, $outpath, {HROFFSET => 0}); print "FILE CREATION: ", Dplib::simplefn($outpath), "\n"; return; } #/**---------------------------------------------------------------------- # @sub write_opnScn_v434 # # Write an opnGns format file (file type = opnScn). This version is used # to attempt to meld v4.3.2, v4.3.3, and v4.3.4 style packets. # NOTE: The oldest style of scintillation data (the oldest format: # obstype 3 and 4 phases contained in the 'SCNThrat' section of the parser. # # @parameter $outfiles -- Root of the TRIG data parser, containing # @ scintillation data in the SCNThrat, SCNThratII, and # @ SCNThratIII sections # @ $ant -- antenna ID # @ $outpath: Output file name # @output File $outpath written # @exception Will toss an exception in case of real trouble. # ----------------------------------------------------------------------*/ sub write_opnScn_v434 { my $outfiles = shift; my $ant = shift; my $outpath = shift; my $v432 = $outfiles->{'SCNThrat432'}{$ant}; # obstype 1-4 GPS data my $v433 = $outfiles->{'SCNThrat433'}{$ant}; # obstype 5,6 GLONASS data my $v434 = $outfiles->{'SCNThrat434'}{$ant}; # obstype 7,8 GPS data # holder for GLONASS and GPS data in the Level1aTools::sort_opnGns_v2 format: # $alldata{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr my %alldata; # GLONASS data are reconstructed solely from the obstype 5,6 data # which are one-second packets with 100Hz I and Q data. # # GPS data are reconstructed from the obstype 7,8 data which # are 1 second phase packets containing 50Hz data, along with the # older SNRs from obstype 1,2 packets. These packets are 10 seconds # long and contain 50Hz SNRs. What a tangled web we weave! # 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 [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' (L2C6L) bit-free pilot signal [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 [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 ); # ## First process GLONASS data into the intermediate hash format required by ## Level1aTools::sort_opnGns_v2: ## $data{$conid}{$prn}{$tpkt} = $hdr . $lrhdr . $hrhdr . $lrvals . $hr # my $status = 0x0; # closed loop, no flags # Determine which SATELLIT. file to use (the latest one) #my $Ivs = 'I14'; # hard wired! #my $satinfo = SatInfo->new("SATELLIT.$Ivs.bern")->parse; # Jun modified, use new file provided by Yong my $Ivs = 'I20'; # Hard-wired to IGS 14 coordinates my $satinfo = SatInfo->new("SATELLIT_$Ivs.SAT")->parse; foreach my $fullprn (sort { $a cmp $b} keys %$v433) { my ($conid, $prn) = unpack "a a2", $fullprn; $prn += 0; my @times = sort { $a <=> $b} keys %{$v433->{$fullprn}}; my ($svn, $freq_slot, $f1, $f2) = $satinfo->find_svn($fullprn, $times[0]); TIME: foreach my $time (@times) { if (!exists($v433->{$fullprn}{$time}{'L1'}) || !exists($v433->{$fullprn}{$time}{'L2'})) { print "Both L1 and L2 data not present for $fullprn at $time, skipping\n"; next TIME; } my ($dt1, $hrfaz1, $mdl1, $snr1, $signal_type1) = eval { @{$v433->{$fullprn}{$time}{'L1'}} }; die "write_opnScnII: Error reading L1 data for $time, $fullprn, quitting" if ($@); my ($dt2, $hrfaz2, $mdl2, $snr2, $signal_type2) = eval { @{$v433->{$fullprn}{$time}{'L2'}} }; die "write_opnScnII: Error reading L2 data for $time, $fullprn, quitting" if ($@); my $rate = $dt1->nelem; if ($dt1->at(0) != $dt2->at(0)) { die "L1 and L2 times do not match: PRN$fullprn T$time"; } my $tf = long(($dt1/1e-7)->rint); my $mostsig = byte(($tf & 0x00ff0000) >> 16); my $midsig = byte(($tf & 0x0000ff00) >> 8); my $lowsig = byte ($tf & 0x000000ff); # Now pack all this together for one second my $rec = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $rate); # Format Epoch Data Record (v2) Use @hrdesc defined # at a package level for speed Add low rate record count, high # rate record count, then high rate descriptors $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_glonass)) . join ('', @hrdesc_l2c_glonass); my $hrmat = cat (double($mostsig), double($midsig), double($lowsig), $hrfaz1, $hrfaz2, $snr1, $snr2); my $hrrec = pack ("(CCCddSS)>" x $rate, $hrmat->xchg(0,1)->flat->list); $alldata{'R'}{$prn}{$time} = $rec.$hrrec; } } # ## Now we have the GLONASS data (obstype 5,6) taken care of, let's look at the GPS data ## in the two packet types at two different cadences # # Time fraction fields for GPS SNRs my $rate = 50; # GPS my $hrt = sequence($rate)/$rate; # 0, 0.02, ..., 0.98 my $tf = long(($hrt/1e-7)->rint); my $mostsig = byte(($tf & 0x00ff0000) >> 16); my $midsig = byte(($tf & 0x0000ff00) >> 8); my $lowsig = byte ($tf & 0x000000ff); my $needed_elems = (10*$rate); # 500 foreach my $fullprn (sort { $a cmp $b} keys %$v432) { my ($conid, $prn) = unpack "a a2", $fullprn; $prn += 0; if ($prn == 0) { warn "PRN equals zero in SCNThrat data (very odd...)" if ($prn == 0); next; } next if ($conid eq 'R'); # only GPS please! my @times = sort { $a <=> $b} keys %{$v432->{$fullprn}}; my ($svn, $freq_slot, $f1, $f2) = $satinfo->find_svn($fullprn, $times[0]); TIME: foreach my $t10 (@times) { # ## Check to make sure required data for this time step are there # next TIME if ( !defined ($v432->{$fullprn}{$t10}{'1'}) || !defined ($v432->{$fullprn}{$t10}{'2'}) ); if ( ($v432->{$fullprn}{$t10}{'1'}[0]->nelem != $needed_elems) || ($v432->{$fullprn}{$t10}{'2'}[0]->nelem != $needed_elems) ) { next TIME; } # Now the problem is that we have 10 second epochs from the TRIG, but we need 1 second epochs for the opnGns format for (my $i=0;$i<10;$i++) { my $sidx = $i*$rate; my $eidx = ($i*$rate)+$rate-1; my $time = $t10 + $i; my $l1s = $v432->{$fullprn}{$t10}{'1'}[0]->mslice([$sidx, $eidx]); my $l2s = $v432->{$fullprn}{$t10}{'2'}[0]->mslice([$sidx, $eidx]); # Now pack all this together for one second my $rec = pack ("(LCCASCS)>", $time, $status, $ant, $conid, $svn, $prn, $rate); # Determine GPS L2C vs. L2P my $signal_type2 = $v434->{$fullprn}{$time}{'L2'}[3]; if (!defined($signal_type2)) { print "Could not find corresponding v434 GPS scintillation packet for $fullprn at $time\n"; next; } # Format Epoch Data Record (v2) Use @hrdesc defined # at a package level for speed Add low rate record count, high # rate record count, then high rate descriptors if ($signal_type2 == 0x30) { # L2P (0x30) instead of L2CL (0x2C) $rec .= pack ('CC', 0, scalar(@hrdesc_l2p_gps)) . join ('', @hrdesc_l2p_gps); } else { # default to GPS L2C $rec .= pack ('CC', 0, scalar(@hrdesc_l2c_gps)) . join ('', @hrdesc_l2c_gps); } # At this point, we have the 1 second SNR data, but lack the one second phase data. # Pull this from the other hash type (v434) # $v434->{$fullprn}{$time}{$freqid} = [$dt, $hrfaz, $mdl, $signal_type] my $l1f_p1 = $v434->{$fullprn}{$time}{'L1'}[1]; # dt, -0.5 .. 0.48 :-( my $l1f_p2 = $v434->{$fullprn}{$time+1}{'L1'}[1]; next if (!defined($l1f_p1) || !defined($l1f_p2)); # Need to find the values from $l1f_p1 and $l1f_p2, which run from $time-0.5 to $time+1+0.48 # to $hrt, which is $time+0 .. $time+0.98 my $l1f = append($l1f_p1->mslice([25,-1]), $l1f_p2->mslice([0,24])); my $l2f_p1 = $v434->{$fullprn}{$time}{'L2'}[1]; # dt, -0.5 .. 0.48 :-( my $l2f_p2 = $v434->{$fullprn}{$time+1}{'L2'}[1]; next if (!defined($l2f_p1) || !defined($l2f_p2)); my $l2f = append($l2f_p1->mslice([25,-1]), $l2f_p2->mslice([0,24])); my $hrmat = cat (double($mostsig), double($midsig), double($lowsig), $l1f, $l2f, $l1s, $l2s); my $hrrec = pack ("(CCCddSS)>" x $rate, $hrmat->xchg(0,1)->flat->list); $alldata{'G'}{$prn}{$time} = $rec.$hrrec; } } } Level1aTools::sort_opnGns_v2 (\%alldata, $outpath, {HROFFSET => 0}); print "FILE CREATION: ", Dplib::simplefn($outpath), "\n"; return; } 1;