Subs::ExtractTiming (version 0.0)


package Subs::ExtractTiming;
##############################################################################
#
# DESCRIPTION:
#
# HISTORY: $Log: ExtractTiming.pm,v $
# HISTORY: Revision 1.22  2017/01/04 16:35:43  apsop
# HISTORY: Update for end of 2016 leap second (3.17.09)
# HISTORY:
# HISTORY: Revision 1.21  2015/06/04 20:16:42  apsop
# HISTORY: Updates for 3.16.09:
# HISTORY: * Address upcoming leap second end of June 2015.
# HISTORY: * UVOT shift+add exposure maps now fixed for dead time correction.
# HISTORY:
# HISTORY: Revision 1.20  2012/07/11 06:18:18  apsop
# HISTORY: Commit of 2012-07-01 leap second correction.
# HISTORY:
# HISTORY:
# HISTORY: JRG as apsop on sdc 2012-06-28: Correct for 2012-07-01 leap
# HISTORY:    second, based on proc3.12.9 version.
# HISTORY:
# HISTORY: Revision 1.19  2009/12/18 14:37:35  apsop
# HISTORY: This release contains several new items:
# HISTORY:   1. call to Subs::UvotGraspCorr for apsect correction
# HISTORY:   2. call to Subs::UvotProduct, which generates gifs light curves
# HISTORY:   3. call to Subs::XrtGrbLc, which generates gifs light curves for GRB based on XRT data
# HISTORY:
# HISTORY: Revision 1.18  2006/09/20 20:28:01  apsop
# HISTORY: Do not extract timing information from the uvot engineering packets, ApID 852.
# HISTORY:
# HISTORY: Revision 1.17  2006/05/06 21:14:06  apsop
# HISTORY: Remove column specific keywords copied from first engineering file extension
# HISTORY:
# HISTORY: Revision 1.16  2006/01/04 19:08:52  apsop
# HISTORY: Filter out UTCF values with abs(UTCF) > 100.
# HISTORY:
# HISTORY: Revision 1.15  2005/12/22 13:54:17  apsop
# HISTORY: Leap second correction for first bit of 2006, before MOC jams UTCF.
# HISTORY:
# HISTORY: Revision 1.14  2005/06/02 12:43:48  apsop
# HISTORY: Fix bug in getting tstart and tstop for second extension.
# HISTORY:
# HISTORY: Revision 1.13  2005/05/31 15:28:55  apsop
# HISTORY: Speed up by using CFITSIO/FTOOL capabilities (row filtering, min/max, heap
# HISTORY: sort).
# HISTORY:
# HISTORY: Revision 1.12  2005/03/15 20:02:59  apsop
# HISTORY: Fix bug in call for getting sc eng hk file name. Remove columns for which UTCF value is set to zero.
# HISTORY:
# HISTORY: Revision 1.12  2005/03/15 19:03:26  apsop
# HISTORY: Check for presence of xrt no position message when making tdrss cat file.
# HISTORY:
# HISTORY: Revision 1.11  2004/11/02 21:17:49  apsop
# HISTORY: Write DATE keyword into timing file.
# HISTORY:
# HISTORY: Revision 1.10  2004/10/25 14:50:01  apsop
# HISTORY: Write standard timing keywords into timing file.
# HISTORY:
# HISTORY: Revision 1.9  2004/08/13 15:41:06  apsop
# HISTORY: Fix to algorithm for thinning time info extensions.
# HISTORY:
# HISTORY: Revision 1.8  2004/07/23 19:50:10  apsop
# HISTORY: Increase allowed string length for deleting rows to 511.
# HISTORY:
# HISTORY: Revision 1.7  2004/07/19 16:00:37  apsop
# HISTORY: Check if string for deleting rows in timing file is too long.
# HISTORY:
# HISTORY: Revision 1.6  2004/07/11 20:39:38  apsop
# HISTORY: Lots of changes to account for bugs in first version
# HISTORY:
# HISTORY: Revision 1.5  2004/06/11 19:51:08  apsop
# HISTORY: Handle case where there are 3 or less attitude entries
# HISTORY:
# HISTORY: Revision 1.4  2004/06/08 00:08:20  apsop
# HISTORY: Fix for handling case with no attitude info.
# HISTORY:
# HISTORY: Revision 1.3  2004/06/07 17:01:27  apsop
# HISTORY: Use ftmerge over fmerge, cause it does not have a file number limit.
# HISTORY:
# HISTORY: Revision 1.2  2004/05/28 19:40:57  apsop
# HISTORY: Bug fixes and enhancement to make first working version
# HISTORY:
# HISTORY: Revision 1.1  2004/05/06 20:02:25  dah
# HISTORY: New module to extract timing info (like UTCF) from the telemetry.
# HISTORY:
# HISTORY:
#
# VERSION: 0.0
#
#
##############################################################################

use Subs::Sub;

@ISA = ("Subs::Sub");
use strict;

sub new {
    my $proto=shift;
    my $self=$proto->SUPER::new();

    $self->{DESCRIPTION}="Extracting Timing Information from Telemetry";

    return $self;
}

##################
# METHODS:
##################

sub body {
    my $self=shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();

    ########################################################
    # Extract UTC adjustment parameters from the hk data
    ########################################################
    my $params_tmp = 'utcf_params.tmp';
    my $tmp_file = 'time_files.tmp';

    my $schk = $filename->get('scenhk', 'swift', '', 0);

    if( -f $schk ){
      my $specs = "[col EXTNAME='UTCF_PARAMS'; " .
	           'TIME; ENABLED_A = (INT)(SUDLUPSTAT1_A/2**14)%2==1;' .
		     'ENABLED_B = (INT)(SUDLUPSTAT1_B/2**14)%2==1;' .
		       'DIRECTION_A = (INT)(SUDLUPSTAT1_A/2**13)%2==1;' .
			 'DIRECTION_B = (INT)(SUDLUPSTAT1_B/2**13)%2==1;' .
			   ' INTERVAL == SUDLADJINTCNT ';

      my $fstruct = Util::Ftool->new("fstruct")
	                       ->verbose(1)
                               ->params({infile  => $schk,
					 outfile => 'STDOUT',
					 colinfo => 'no'})
                               ->run();

      my @hdus = $fstruct->stdout() =~ /(hk00[\d,a-c]x00\d)/gs;

      if(@hdus){
	###############################################################################
	# Remove column specific keywords copied from first engineering file extension
	###############################################################################
	my $tmpfits = Util::FITSfile->new($schk, $hdus[0]);
	my %keywords = $tmpfits->keywords();
	my @badkeys = grep /^T(YE|RD)(LO|HI)/, keys(%keywords);
	@badkeys = map '-'.$_ , @badkeys;

	open TMP, ">$tmp_file";
	print TMP $schk .'['.  shift(@hdus) .']'. $specs . join('; ', @badkeys) .']', "\n";
	foreach (@hdus){
	  print TMP $schk . "\[$_\]" . $specs .']', "\n";
	}
	close TMP;

	$log->entry('Extract utcf parameters from house keeping: '. join(' ',@hdus));

	Util::HEAdas->new('ftmerge')
	            ->params({infile => '@'.$tmp_file,
			      outfile => $params_tmp})
		    ->run();

	unlink $tmp_file;
      }
    }

    ##############################################################
    # Extract the value of UTCF as a function of MET from LPD and
    # teriary header files
    ##############################################################

    my $extract = Util::Tool->new($procpar->read('fitspackets'), 'extract_utcf');
    my $i = 0;
    my $utcf_file = $filename->get('timedata', 'swift', '', 0);

    open TMP, ">$tmp_file";
    # list of ApIDs not to use.  This includes bat cpio packets and uvot engineering packets.
    my %exclude_apid = map { $_ => 1 } (299, 417, 481, 484, 568, 570, 852);
    foreach my $tel_file ( $filename->get('telemetry', '*', 'ldp', '*'),
			   $filename->get('telemetry', '*', 'head3', '*') ){

      my ($inst, $mode, $apid) = $filename->parse($tel_file, 'telemetry');
      next if $exclude_apid{$apid};

      my $tmp_utcf = 'utcf_' . ++$i . '.tmp';

      $log->entry("Extracting UTCF values from file $tel_file.");

      $extract->command_line("-infile $tel_file",
			     "-apid_list /aps/lists/swift_apids.list",
			     "-outfile $tmp_utcf")
              ->run();

      print TMP $tmp_utcf."[1][abs(UTCF) > 1E-9 && abs(UTCF) < 100]\n"  unless $extract->had_error();
    }
    close TMP;

    if($i){
      Util::HEAdas->new('ftmerge')
	          ->params({infile => '@'.$tmp_file,
			    outfile => $utcf_file})
		  ->run();
      unlink glob('utcf_[0-9]*.tmp'), $tmp_file;

      if( -f $params_tmp ){
	$log->entry('Appending utcf parameters extension to timing file.');

	Util::FITSfile->new($params_tmp, 'UTCF_PARAMS')
	              ->append_to($utcf_file);
	unlink $params_tmp;
      }

      ######################################################################
      # Trim the tables so that only the rows where changes occur are kept.
      # Always keep the first and last row.
      ######################################################################

      my $now = Util::Date->new();
      my $Tnow = $now->date() .'T'. $now->time();
      my $utcf_fits = Util::FITSfile->new($utcf_file, 0);
      $utcf_fits->keyword('DATE', $Tnow);

      $utcf_fits->ext(1);
      $utcf_fits->cols('TIME')->sort('heap', 'unique');

      my ($sum,$mean,$sigma,$min,$max) = $utcf_fits->stats;

      ##################################################
      # Check if we need to do a leap second correction
      # The work is done by ftcopying the file using CFITSIO
      # column-filtering syntax.
      ##################################################
      my $correct_infile;

      # @CORR defines the times when UTCF needs to be adjusted for leap seconds.
      # START is the time of the leap second. STOP is the end time of the
      # 'correction'.
      # The entries with UTCF_TEST allow correcting for the leap second without
      # knowing precisely when the MOC updates UTCF. This lets us update @CORR
      # in advance of the leap second and the filtering will appropriately
      # adjust the UTCF from the time of the leap second to when the change
      # in UTCF appears in the telemetry. STOP should be chosen after the
      # expected time of the MOC update which is normally a FOT monitored
      # pass soon after the leap second, but may be on a later day for instance
      # in the event of a holiday.
      my @CORR = (
          { START => '2006-01-01T00:00:00',
            STOP  => '2006-01-01T02:58:59',
          },
          { START => '2009-01-01T00:00:00',
            STOP  => '2009-01-02T17:33:48',
          },

          { # As of the leap second on 2012-07-01, UTCF changes from ~-8.3 to ~-9.3,
            # but the MOC doesn't upload the new value right away.  So if TIME is
            # within a 4-day window and still has the old value, adjust it ourselves.
            START => '2012-07-01T00:00:00',
            STOP  => '2012-07-04T00:00:00',
            UTCF_TEST => -9.0,
          },

          { # As of the leap second on 2015-07-01, UTCF changes from ~-13.6 to ~-14.6
            START => '2015-07-01T00:00:00',
            STOP  => '2015-07-03T00:00:00',
            UTCF_TEST => -14.0,
          },

          { # As of the leap second on 2017-01-01, UTCF changes from ~-17.4 to ~-18.4
            START => '2017-01-01T00:00:00',
            STOP  => '2017-01-03T00:00:00',
            UTCF_TEST => -18.0,
          },
      );

      foreach my $corr (@CORR) {

          my $corr_start = Util::Date->new($corr->{START});
          my $corr_stop  = Util::Date->new($corr->{STOP});
          my $met_start = $corr_start->seconds();
          my $met_stop  = $corr_stop->seconds();
          if( $met_start < $max && $met_stop > $min ){
	     $correct_infile = $utcf_file . "[UTCF]" .
	              "[col TIME; UTCF = TIME > $met_start && TIME < $met_stop";
             if ($corr->{UTCF_TEST}) {
                 $correct_infile .= " && UTCF > $corr->{UTCF_TEST}";
             }
             $correct_infile .= ' ? UTCF-1 : UTCF; SOURCE; APID]';
             $log->entry("ExtractTiming: plan to correct using $correct_infile");
          }
      }

      if ($correct_infile){
	  $log->entry("ExtractTiming.pm: Correcting ${utcf_file}[UTCF] for leap second");
        my $tmpfile = 'leap_sec_corr.tmp';
        my $filter = Util::HEAdas->new('ftcopy')
                        ->params({
                            infile => $correct_infile,
                            outfile => $tmpfile,
                        })
                        ->run;
         rename($tmpfile, $utcf_file)
             or $log->error(1, "unable to rename $tmpfile to $utcf_file [$!]");
      }

      {
        # use row filtering instead of dumping and parsing table and code
        my $firstLast = '(#ROW==1).or.(#ROW==#NAXIS2)';
        my $changed = '(UTCF!=UTCF{-1}).or.(UTCF!=UTCF{1})';
        my $infile = $utcf_file . "[UTCF][$firstLast.or.$changed]";
        my $tmpfile = 'filtered.tmp';
        my $filter = Util::HEAdas->new('ftcopy')
                        ->params({
                            infile => $infile,
                            outfile => $tmpfile,
                        })
                        ->run;
         rename($tmpfile, $utcf_file)
             or $log->error(1, "unable to rename $tmpfile to $utcf_file [$!]");
      }

      $utcf_fits->begin_many_keywords();
      $utcf_fits->keyword('DATE', $Tnow);
      $utcf_fits->keyword('EXTNAME', 'UTCF');
      $utcf_fits->keyword('TIMESYS', 'TT', 'time measured from');
      $utcf_fits->keyword('MJDREFI', 51910, 'MJD reference day');
      $utcf_fits->keyword('MJDREFF', 7.428703700000000E-04,
			  'MJD reference (fraction of day)');
      $utcf_fits->keyword('CLOCKAPP', 'F', 'default');
      $utcf_fits->keyword('TIMEUNIT', 's', 'unit for time keywords');

      $utcf_fits->keyword('TSTART', $min);
      $utcf_fits->keyword('TSTOP', $max);
      $utcf_fits->end_many_keywords();

      if( $utcf_fits->nhdus() > 2 ){
	$utcf_fits->ext(2);

	$utcf_fits->begin_many_keywords();
	$utcf_fits->keyword('DATE', $Tnow);
	$utcf_fits->keyword('TIMESYS', 'TT', 'time measured from');
	$utcf_fits->keyword('MJDREFI', 51910, 'MJD reference day');
	$utcf_fits->keyword('MJDREFF', 7.428703700000000E-04,
			    'MJD reference (fraction of day)');
	$utcf_fits->keyword('CLOCKAPP', 'F', 'default');
	$utcf_fits->keyword('TIMEUNIT', 's', 'unit for time keywords');

	##########################################################################
	# The TSCALE keyword results in the wrong format for the INTERVAL column
	# so delete it if it has the default value
	##########################################################################
	my $tscal2 = $utcf_fits->keyword('TSCAL2');
	if( $tscal2 && $tscal2 == 1 ){
	  $utcf_fits->keyword('-TSCAL2' => 0);
	}

	my $nrows = $utcf_fits->nrows()-1;
	if($nrows>1){
	  $utcf_fits->cols('TIME')->sort('heap');

          my $rowFilter = '(#ROW==1).or.(#ROW==#NAXIS2)'
                    . '.or.(INTERVAL!=INTERVAL{1})';
          my $pre = '{-1}';
          foreach my $c (qw(INTERVAL ENABLED_A ENABLED_B
                            DIRECTION_A DIRECTION_B)) {
             $rowFilter .= ".or.($c.ne.$c$pre)";
          }
          my $filter = Util::HEAdas->new('ftcopy')
                 ->params({
                      infile => $utcf_fits->fullname . "[$rowFilter]",
                      outfile => $params_tmp,
                 })->run;
          rename($params_tmp, $utcf_file)
             or $log->error(1, "unable to rename $params_tmp to $utcf_file [$!]");
        }

	{
	  my ($sum,$mean,$sigma,$min,$max) = $utcf_fits->stats;
	  $utcf_fits->keyword('TSTART', $min);
	  $utcf_fits->keyword('TSTOP', $max);
	}
	$utcf_fits->end_many_keywords();
      }
    }
    else{
      unlink $params_tmp if -f $params_tmp;
      unlink $tmp_file if -f $tmp_file;
    }
}

1;