package Subs::Attitude; ############################################################################## # # DESCRIPTION: Extract attitude data. Swift attitude data appears # DESCRIPTION: in several places in the telemetry. This subroutine # DESCRIPTION: collects all the attitude data into a single quaternion- # DESCRIPTION: based attitude file and removed redundant information. # DESCRIPTION: # DESCRIPTION: The main source of attitude data is the "ACS packets" # DESCRIPTION: (APID 404). The BAT also places ACS records in the # DESCRIPTION: headers of each LDP and in the body of a special # DESCRIPTION: LDP generated suring slews. # DESCRIPTION: The XRT places ACS records in the header for each CCD # DESCRIPTION: frame, however these are represented by "floats" instead # DESCRIPTION: of "doubles", so they are not used if higher accuracy # DESCRIPTION: data are available. # DESCRIPTION: # DESCRIPTION: The sub checkDatabaseForNominalPointing reads the database # DESCRIPTION: file "/sw/nominal_pnt.rdb" from the directory indicated by # DESCRIPTION: env var $RDB_TABLES_DIR (normally set in aps.config; # DESCRIPTION: defaults to /aps/db/rdb/tables). If the sequence is found # DESCRIPTION: in this table, those values are used for the nominal # DESCRIPTION: pointing; this allows for manual override when the pointing # DESCRIPTION: cannot be reliably determined automatically. Otherwise, the # DESCRIPTION: mean pointing is calculated and compared with the requested # DESCRIPTION: pointing (SACSREQ, read from the engineering housekeeping # DESCRIPTION: file), if these disagree an E2 error is logged and an email # DESCRIPTION: notification is sent. Currently the tolerance is 10 arcmin. # DESCRIPTION: Someone should then review the sequence and add the correct # DESCRIPTION: position to nominal_pnt.rdb, then rerun the sequence. # DESCRIPTION: # DESCRIPTION: Here is a pseudocode outline of the algorithm used to # DESCRIPTION: determine the nominal attitude: # DESCRIPTION: # DESCRIPTION: IF obs in ADB [nominal pointing database] # DESCRIPTION: pnt_db = ADB[obs] # DESCRIPTION: ELSE # DESCRIPTION: pnt_sacs = determineSACSREQ(sw<obs>sen.hk) # DESCRIPTION: # 1-2% of the time, determineSACSREQ fails... # DESCRIPTION: pnt_aspect = aspect(sw<obs>[pu]at.fits) # DESCRIPTION: IF pnt_sacs AND pnt_aspect THEN # DESCRIPTION: IF separation(pnt_sacs, pnt_aspect) > 10 arcmin THEN # DESCRIPTION: E2 ERROR SACSREQ and aspect positions inconsistent # DESCRIPTION: SEND EMAIL NOTIFICATION # DESCRIPTION: ENDIF # DESCRIPTION: IF minRollDelta(pnt_sacs, pnt_aspect) > 2 deg THEN # DESCRIPTION: E1 ERROR # DESCRIPTION: ENDIF # DESCRIPTION: ELSEIF pnt_aspect # DESCRIPTION: IF pnt_aspect has ERROR # DESCRIPTION: E2 ERROR Unable to determine nominal pointing # DESCRIPTION: SEND EMAIL NOTIFICATION # DESCRIPTION: ENDIF # DESCRIPTION: ENDIF # DESCRIPTION: ENDIF # DESCRIPTION: # # HISTORY: $Log: Attitude.pm,v $ # HISTORY: Revision 1.63 2016/11/04 13:04:49 apsop # HISTORY: Add UTCFINIT to attitude files so prefilter can use it to address clock drift. Pass prefilter jump corrected (.pat) instead of uncorrected (.sat) attitude if available. # HISTORY: # HISTORY: Revision 1.62 2014/08/14 09:19:47 apsop # HISTORY: Send an email notification to watchers if nominal pointing cannot be # HISTORY: determined by any method. Include job title in both nominal-pointing # HISTORY: emails. # HISTORY: # HISTORY: Revision 1.61 2014/02/27 08:26:31 apsop # HISTORY: ATTFLAG comment explains value. # HISTORY: Sequence and seqprocnum included on subject line of # HISTORY: aspect/SACSREQ Inconsistent email. # HISTORY: # HISTORY: Revision 1.60 2013/04/29 22:08:50 apsop # HISTORY: Add seqprocnum to the email sent out if aspect & sacsreq are inconsistent # HISTORY: # HISTORY: Revision 1.59 2013/03/28 09:00:58 apsop # HISTORY: Added more comments. # HISTORY: # HISTORY: Revision 1.58 2013/03/14 10:02:11 apsop # HISTORY: Math.pm has been added to Util, so use Util::Math # HISTORY: # HISTORY: Revision 1.57 2012/09/08 07:03:56 apsop # HISTORY: Added object, sequence, target to email alert for inconsitent # HISTORY: SACSREQ/aspect positions. # HISTORY: # HISTORY: Revision 1.56 2012/09/05 03:44:12 apsop # HISTORY: Modified command that reads nominal_pnt.rdb to not return header. # HISTORY: # HISTORY: Revision 1.55 2012/08/15 06:53:37 apsop # HISTORY: Check database nominal_pnt.rdb for manual pointing override; # HISTORY: compare mean aspect solution with requested pointing (SACSREQ); # HISTORY: add time keywords (TSTART, TSTOP, DATE, etc.) to all extensions. # HISTORY: [Bob W.] # HISTORY: # HISTORY: Revision 1.54 2007/10/10 20:32:29 apsop # HISTORY: Put tolerance for target-pointing GTI back to 1.0 . # HISTORY: # HISTORY: Revision 1.53 2007/10/08 21:03:56 apsop # HISTORY: For targ pointing GTI, reduced tolerance to 0.1 degrees, and rename STDGTI extension to GTI. # HISTORY: # HISTORY: Revision 1.52 2007/09/28 17:14:48 apsop # HISTORY: Change writing of ATTFLAG keyword so that it is a string rather than an integer. # HISTORY: # HISTORY: Revision 1.51 2007/09/17 15:32:15 apsop # HISTORY: Set chatter for attjumpcorr to 3. # HISTORY: # HISTORY: Revision 1.50 2007/09/11 18:22:23 apsop # HISTORY: Turn up chatter on attjumpcorr tasked to 5. # HISTORY: # HISTORY: Revision 1.49 2007/04/01 19:08:48 apsop # HISTORY: Set chatter for attjumpcorr to 3. # HISTORY: # HISTORY: Revision 1.48 2007/02/22 18:35:33 apsop # HISTORY: Fix subtle, intemittent bug in running acs2fits w/ setting $_ . # HISTORY: # HISTORY: Revision 1.47 2007/02/08 21:33:12 apsop # HISTORY: Write ATTFLAG keyword to all extensions. # HISTORY: # HISTORY: Revision 1.46 2007/01/31 20:22:14 apsop # HISTORY: Do attitude correction with attjumpcorr and set attitude correction flags. # HISTORY: # HISTORY: Revision 1.45 2006/10/01 18:55:57 apsop # HISTORY: Insert date keywords into the attitude file so that new version of aspect task will work properly. # HISTORY: # HISTORY: Revision 1.44 2006/06/28 18:54:42 apsop # HISTORY: Create new target pointing GTI file. # HISTORY: # HISTORY: Revision 1.43 2006/06/15 22:11:52 apsop # HISTORY: Make pointing and slew GTIs more restrictive, so that times beyond end measurements are excluded. # HISTORY: # HISTORY: Revision 1.42 2006/04/27 16:28:12 apsop # HISTORY: Make a lack of any attitude file a fatal error. # HISTORY: # HISTORY: Revision 1.41 2006/03/13 17:39:17 apsop # HISTORY: Use XRT attitude info in 991 seqs. Fix bug in choosing GTI for mean pointing determination. # HISTORY: # HISTORY: Revision 1.40 2006/01/29 19:37:26 apsop # HISTORY: Add 60 secs of buffer at each end of the range for the att/orbit file. # HISTORY: # HISTORY: Revision 1.39 2005/11/20 21:26:10 apsop # HISTORY: Unexpected OO behaviour for $inter_file requires seperate call to append_to. # HISTORY: # HISTORY: Revision 1.38 2005/11/20 20:30:03 apsop # HISTORY: Check for presence of BUS_V column before making and using ACS_DATA extension in attitude file. # HISTORY: # HISTORY: Revision 1.37 2005/11/08 16:58:18 apsop # HISTORY: Change fdiff to ftdiff and set reltol to 10e-8. Use caldb to get alignment file in prefilter and abberator. # HISTORY: # HISTORY: Revision 1.36 2005/09/26 21:32:35 apsop # HISTORY: Make not_pointing GTI file. # HISTORY: # HISTORY: Revision 1.35 2005/07/29 14:04:03 apsop # HISTORY: look for ACS packets that are labeled as head2 # HISTORY: # HISTORY: Revision 1.34 2005/06/01 17:38:51 apsop # HISTORY: More robust algorithm for selecting GTI file to use for determining the mean pointing. # HISTORY: # HISTORY: Revision 1.33 2005/04/19 16:05:48 apsop # HISTORY: Check for existence of gti files before setting keywords. # HISTORY: # HISTORY: Revision 1.32 2005/03/25 21:25:34 apsop # HISTORY: Fix bug in setting TSTOP for GTIs. Set gti extname to GTI. # HISTORY: # HISTORY: Revision 1.31 2005/03/25 20:18:04 apsop # HISTORY: Set TSTART and TSTOP for GTI files. # HISTORY: # HISTORY: Revision 1.30 2005/03/15 20:04:42 apsop # HISTORY: Fix bug in call for getting xrt eng hk file name. # HISTORY: # HISTORY: Revision 1.30 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.29 2005/02/18 01:53:09 apsop # HISTORY: Only use secondary attitude data if primary data is not available. Fix up xrt modes. # HISTORY: # HISTORY: Revision 1.28 2005/02/14 19:23:39 apsop # HISTORY: Remove batatt call. Instead pickup bat attitude info from bat2fits output. # HISTORY: # HISTORY: Revision 1.27 2005/02/10 02:52:14 apsop # HISTORY: Get start and stop times from attitude file, as StartEndTimes has not been run yet. # HISTORY: # HISTORY: Revision 1.26 2005/02/08 18:24:17 apsop # HISTORY: Abberation correction using abberator tool with time test. Move calculation of attorb file to this class as it is needed for abberator. # HISTORY: # HISTORY: Revision 1.25 2005/01/12 17:19:13 apsop # HISTORY: Exclude safeholds from settling and pointing GTIs. Include 10arcmin flag in pointing GTI. # HISTORY: # HISTORY: Revision 1.24 2004/12/10 02:16:09 apsop # HISTORY: Do a diff of attitude files after normalization, so we have a record of what changed. # HISTORY: # HISTORY: Revision 1.23 2004/12/02 18:55:50 apsop # HISTORY: Normalize quaternions before passing to aspect. # HISTORY: # HISTORY: Revision 1.22 2004/11/02 21:12:15 apsop # HISTORY: Set DATE keyword in attitude file. # HISTORY: # HISTORY: Revision 1.21 2004/09/03 00:26:01 apsop # HISTORY: Temporarily remove bat lc attitude information from attitude file. # HISTORY: # HISTORY: Revision 1.20 2004/08/16 15:23:09 apsop # HISTORY: Turn on writing of HISTORY keywords. # HISTORY: # HISTORY: Revision 1.19 2004/07/06 19:59:32 apsop # HISTORY: Add test for existence of nfi GTI file. # HISTORY: # HISTORY: Revision 1.18 2004/06/29 14:34:02 apsop # HISTORY: New ontarget gti is pointing ANDed with nfi. Use as input to aspect tool. # HISTORY: # HISTORY: Revision 1.17 2004/06/14 14:26:23 apsop # HISTORY: Write in name of alignment file into attitude file. # HISTORY: # HISTORY: Revision 1.16 2004/06/10 19:55:03 apsop # HISTORY: Add updating of TSTART/TSTOP keywords in second extension. # HISTORY: # HISTORY: Revision 1.15 2004/05/28 19:16:07 apsop # HISTORY: Changes to attitude file format # HISTORY: # HISTORY: Revision 1.14 2004/05/06 20:02:33 dah # HISTORY: Add version number back into the header comments. # HISTORY: # HISTORY: Revision 1.13 2004/05/04 16:30:51 dah # HISTORY: Set ra/dec to 0/0 if not attitude data. Issue error. # HISTORY: # HISTORY: Revision 1.12 2004/04/28 13:45:50 dah # HISTORY: Change xrt header and engineering file modes to fit new standard. # HISTORY: # HISTORY: Revision 1.11 2004/04/16 20:21:18 dah # HISTORY: Begin using embedded history records # HISTORY: # # VERSION: $Revision: 1.63 $ # ############################################################################## use Subs::SwiftSub; use Subs::NominalPointing; use Util::AttTool; use Scalar::Util; use Util::Email; # Math.pm is in the Util directory (copied from heasoft), # but isn't part of the Util package. use Util::Math; @ISA = ('Subs::NominalPointing', 'Subs::SwiftSub'); use strict; use Util::SwiftTags; sub new { my $proto=shift; my $self=$proto->SUPER::new(); my $procpar =$self->procpar(); $self->{DESCRIPTION}="Extracting attitude data"; 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 attitude files from various sources ############################################## my @acs = $self->extract_acs("acs_att_tmp"); my @bat = $filename->get('attitude', 'bat', '', '*'); my @xrt; my @batlc; unless( @acs || @bat ){ $log->entry("No primary attitude data, using secondary sources."); push @xrt, $self->extract_xrt("xrt_att_tmp"); push @bat, $filename->get('attlpd', 'bat', '', '*'); }else{ unlink $filename->get('attlpd', 'bat', '', '*'); my $seq = $jobpar->read('sequence'); if( $seq%1000==991 ){ $log->entry("Using xrt attitude data in 991 sequence."); push @xrt, $self->extract_xrt("xrt_att_tmp"); } } ############################################################# # merge just the attitude files from the ACS type packets ########################################################### my $acslist = Util::FITSlist->new(@acs, @bat); my @acs_merged; unless($acslist->count() == 0) { $log->entry("Merging ". join(" ", $acslist->files()) ); push @acs_merged, 'acs_merged.tmp'; my $merge_out = $acslist->merge($acs_merged[0]); if($merge_out ne $acs_merged[0]) { ############################ # there was only one file ########################### rename $merge_out, $acs_merged[0]; }else{ ########################## # delete the files ########################## foreach ($acslist->files() ) { $log->entry("Deleting $_"); unlink $_; } } }else{ my $index = 0; foreach my $lcfile ($filename->get('batlcatt', 'bat', '', '*'), $filename->get('tdrsslcatt', 'bat', '', '*')){ push @batlc, "bat_lcatt_${index}.tmp"; Util::FITSfile->new($lcfile, 'ATTITUDE', '[col TIME; QPARAM; POINTING; SOURCE(B)=4]') ->copy($batlc[-1]); unlink $lcfile; $index++; } } ####################################### # merge all the attitude data ####################################### my $attlist = Util::FITSlist->new(@acs_merged, @batlc, @xrt); if($attlist->count() ==0) { $log->error(2, 'No attitude files produced, setting ra,dec to 0,0'); $jobpar->set({ra => 0.0, dec => 0.0}); return; } else { $log->entry("Merging ". join(" ", $attlist->files()) ); } my $attitude = $filename->get('attitude', 's'); my $merged=$attlist->merge($attitude, '', 'TIME', 'QPARAM', 'POINTING', 'SOURCE'); ############################################################################### # Make a seperate extension with just the extra columns from the ACS packets ############################################################################### if($merged ne $attitude) { ############################ # there was only one file ########################### rename $merged, $attitude; my $inter_file = Util::FITSfile->new($attitude, 'ATTITUDE'); if($inter_file->find_column('BUS_V')){ $inter_file->specs('[col TIME; POSITION; FLAGS; BUS_V; SOURCE]'); $inter_file->append_to($attitude); my $delcol = Util::Ftool->new('fdelcol') ->params({infile => $attitude.'[1]', confirm => 'no', proceed => 'yes'}); $delcol->params({colname => 'POSITION'})->run(); $delcol->params({colname => 'FLAGS'})->run(); $delcol->params({colname => 'BUS_V'})->run(); } }else{ if(@acs_merged){ my $inter_file = Util::FITSfile->new($acs_merged[0], 'ATTITUDE'); if($inter_file->find_column('BUS_V')){ $inter_file->specs('[col TIME; POSITION; FLAGS; BUS_V; SOURCE]'); $inter_file->append_to($attitude); } } ########################## # delete the files ########################## foreach ($attlist->files() ) { $log->entry("Deleting $_"); unlink $_; } } ########################################## # TIME sort the merged attitude file and # remove overdetermined values ########################################## $log->entry("Sorting and uniqing $attitude"); my $now = Util::Date->new(); my $Tnow = $now->date() .'T'. $now->time(); my $fits = Util::FITSfile->new($attitude, 0); $fits->keyword('DATE', $Tnow); $fits->keyword('MJDREFI', 51910, 'MJD reference day'); $fits->keyword('MJDREFF', 7.428703700000000E-04, 'MJD reference (fraction of day)'); $fits->ext(1); $fits->sort(); if($fits->nhdus() > 2){ $fits->ext(2); $fits->cols('TIME'); $fits->sort('unique') ->keyword('EXTNAME', 'ACS_DATA'); } Util::AttTool->new("att_thin") ->command_line($attitude) ->run(); $log->entry("normalizing quaternions"); my $colfilter = '[col *; QPARAM = QPARAM / sqrt(sum(QPARAM*QPARAM))]'; Util::HEAdas->new('ftcopy') ->params({ infile => $attitude . $colfilter, outfile => 'attitude.tmp', }) ->run; ############################################################### # Run diff so that we have a record in the log of what changed ############################################################### my $diff = Util::HEAdas->new('ftdiff') ->params({infile1 => "$attitude\[1]", infile2 => "attitude.tmp\[1]", reltol => 10e-8}); $diff->run(); rename('attitude.tmp', $attitude); ########################################## # set TSTART, TSTOP, and UTCFINIT in the merged file ########################################## $fits = Util::FITSfile->new($attitude, 1) ->cols("TIME"); my $nrows = $fits->nrows(); my $tstart = $fits->rows(1 )->table(); my $start_date = Util::Date->new($tstart); my $start_str = $start_date->date() . 'T' . $start_date->time(); my $tstop = $fits->rows($nrows)->table(); my $stop_date = Util::Date->new($tstop); my $stop_str = $stop_date->date() . 'T' . $stop_date->time(); my $utcfinit = undef; { # determine UTCFINIT at TSTART my $time_file = $filename->get('timedata', 'swift', '', 0); my (@utcf_times, @utcf); if (-f $time_file) { my $time_fits = Util::FITSfile->new($time_file, 'UTCF'); @utcf_times = $time_fits->cols('TIME')->table(); @utcf = $time_fits->cols('UTCF')->table(); } if (@utcf) { my $itime = 0; while ($tstart > $utcf_times[$itime] && $itime < $#utcf_times) { $itime++ } $utcfinit = $utcf[$itime]; } if (defined($utcfinit)) { $log->entry("found UTCFINIT=$utcfinit at $tstart"); } else { $log->error(1, "unknown UTCFINIT at $tstart"); } } my $nhdus = $fits->nhdus; for (my $i = 0; $i < $nhdus; ++$i) { $fits->ext($i); $fits->keyword('DATE', $Tnow); $fits->keyword('MJDREFI', 51910, 'MJD reference day'); $fits->keyword('MJDREFF', 7.428703700000000E-04, 'MJD reference (fraction of day)'); $fits->keyword('TSTART', $tstart, 'Time of first attitude record'); $fits->keyword('TSTOP' , $tstop, 'Time of last attitude record'); $fits->keyword('DATE-OBS', $start_str, 'Date corresponding to TSTART'); $fits->keyword('DATE-END', $stop_str, 'Date corresponding to TSTOP'); if (defined($utcfinit)) { $fits->keyword('UTCFINIT', $utcfinit, '[s] UTCF at TSTART'); } } $fits->ext(1); my $align_fits = Util::FITSfile->new($filename->fetch_cal('alignment'), 0); $fits->keyword('ALGN_NAM', $align_fits->keyword('FILENAME'), 'Name of alignment teldef file used.'); ################################################# # Make settling and pointing gtis ################################################# my $settling = $filename->get('gti', 's', 'st', 0); my $pointing = $filename->get('gti', 's', 'po', 0); my $not_pointing = $filename->get('gti', 's', 'np', 0); my $ontarget = $filename->get('gti', 's', 'ot', 0); my $nfis = $filename->get('gti', 's', 'nf', 0); my $target_pointing = $filename->get('gti', 's', 'tp', 0); if( (grep /ACS_DATA/, $fits->list_hdus()) ){ my $tempfile = 'acs.tmp'; Util::HEAdas->new('ftsort') ->params({infile => $attitude.'[ACS_DATA]', outfile => $tempfile, unique => 'yes', columns => 'TIME, FLAGS'}) ->run(); my $maketime = Util::Ftool->new('maketime') ->params({infile => $tempfile.'[ACS_DATA]', compact => 'no'}); $log->entry("Producing GTIs for settling and pointing."); $maketime->params({outfile => $settling, prefr => 0, postfr => 1, expr => 'FLAGS == b10x0xxxx'}) ->run(); $maketime->params({outfile => $pointing, prefr => 0, postfr => 0, expr => 'FLAGS == b11x0xxxx'}) ->run(); $maketime->params({outfile => $not_pointing, prefr => 1, postfr => 1, expr => 'FLAGS != b11xxxxxx || FLAGS == bxxx1xxxx'}) ->run(); if( -s $nfis ){ Util::Ftool->new('mgtime') ->params({ingtis => $pointing .' '. $nfis, outgti => $ontarget, merge => 'AND'}) ->run(); }else{ Util::FITSfile->new($pointing) ->copy($ontarget); } unlink $tempfile if -f $tempfile; } foreach my $gtifile ($settling, $pointing, $ontarget, $nfis){ next unless -s $gtifile; my $gtifits = Util::FITSfile->new($gtifile, 1); unless( $gtifits->nrows() ){ $log->error([ 1, ATT_NO_GTI ], "GTI file $gtifile is empty. Deleting."); unlink $gtifile; next; } $gtifits->keyword('EXTNAME', 'GTI'); $gtifits->keyword('TSTART', ($gtifits->cols('START')->table())[0] ); $gtifits->keyword('TSTOP', ($gtifits->cols('STOP')->table())[-1] ); } ########################################## # determine the nominal pointing # using method inherited from superclass ########################################## unless( -f $ontarget ){ if( -f $pointing ) { $ontarget = $pointing; }else{ $ontarget = $self->no_gap_gtis($attitude, 100); } } # check if the nominal pointing database has an entry for the sequence my $pnt_db = $self->checkDatabaseForNominalPointing; if ($pnt_db) { $self->storeNominalPointing($pnt_db); } else { my $pnt_sacs = $self->determineRequestedPointing; my $pnt_aspect = $self->determineAspect($ontarget); if ($pnt_sacs and $pnt_aspect) { my $sacs_str = formatPNT($pnt_sacs); my $aspect_str = formatPNT($pnt_aspect); my $details = "SACSREQ $sacs_str and aspect $aspect_str"; $log->entry("comparing $details"); my $ten_arcmin_deg = sprintf('%.2f', 10 / 60); if (separation_deg($pnt_sacs, $pnt_aspect) > $ten_arcmin_deg) { my $text = "$details positions are inconsistent by more than $ten_arcmin_deg deg"; my $sequence = $jobpar->read('sequence'); my $seqprocnum = $jobpar->read('seqprocnum'); $log->error(2, $text); # (NB: Can't include obsdate because that parameter isn't # set until a later module) $text = $text . "\n\nobject= " . $jobpar->read('object') . "\nsequence= $sequence" . "\ntarget= " . $jobpar->read('target') . "\nprocessing number= ${seqprocnum}" . "\njob title= " . $jobpar->read('job_title') . "\n" ; Util::Email::sendEmail(SUB => $self, SUBJECT => "aspect and SACSREQ inconsistent in ${sequence}.${seqprocnum}", TEXT => $text); } my $rolls = $pnt_sacs->{ROLLS_deg} || [ ]; my $roll_a = $pnt_aspect->{ROLL_deg}; my $matched = undef; my $nonzero = 0; my $max_roll_deg = 2; foreach my $roll (@$rolls) { if ($roll) { ++$nonzero; } my $delta_d = deltaRoll_d($roll, $roll_a); if ($delta_d < 0) { $log->error(1, "bad input: deltaRoll_d($roll, $roll_a)"); } elsif ($delta_d < $max_roll_deg) { $matched = 1; } } if (not $nonzero) { $log->error(1, "no non-zero SACSREQROLL"); } if (not $matched) { $log->error(1, "aspect roll not within $max_roll_deg deg of SACSREQROLL"); } } elsif ($pnt_aspect) { if ($pnt_aspect->{ERROR}) { $log->error(2, 'Unable to determine nominal pointing'); # (NB: Can't include obsdate because that parameter isn't # set until a later module) my $sequence = $jobpar->read('sequence'); my $seqprocnum = $jobpar->read('seqprocnum'); my $text = "Unable to determine nominal pointing by either aspect or SAQSREQ" . "\n\nobject= " . $jobpar->read('object') . "\nsequence= $sequence" . "\ntarget= " . $jobpar->read('target') . "\nprocessing number= ${seqprocnum}" . "\njob title= " . $jobpar->read('job_title') . "\n" ; Util::Email::sendEmail(SUB => $self, SUBJECT => "Unable to determine nominal pointing in ${sequence}.${seqprocnum}", TEXT => $text); } } } if($ontarget =~ /\.tmp$/ ) { #################################### # delete no-gap temporary GTI file #################################### unlink $ontarget; } ################################################################# # GTI of pointings where we are actually pointing in the nominal # direction. Call this 'tp' for 'target pointing. ################################################################# if( -f $pointing ){ my $maketime = Util::Ftool->new('maketime') ->params({infile => $attitude.'[ATTITUDE]', compact => 'no', outfile => 'gti.tmp', prefr => 1, postfr => 1, expr => 'near(POINTING[1],' .$jobpar->read('ra') .',1.0) ' . ' && near(POINTING[2],' .$jobpar->read('dec') .',1.0)' }) ->run(); Util::Ftool->new('mgtime') ->params({ingtis => $pointing .' gti.tmp', outgti => $target_pointing, merge => 'AND'}) ->run(); unlink 'gti.tmp'; my $tpfits = Util::FITSfile->new($target_pointing, 'STDGTI'); $tpfits->keyword('EXTNAME', 'GTI'); } ################################################################# # Only do abberation correction in velocity adding was turned on ################################################################# my $velocity_add = Util::Date->new('2005-01-31T21:57:00'); $self->correct_aberration() if $tstart < $velocity_add->seconds(); my $att_proc = $filename->get('attcorr', 'p'); my $att_corr = Util::HEAdas->new('attjumpcorr')->is_script(1); $att_corr->params({infile => $attitude, outfile => $att_proc, chatter => 3}); $att_corr->run(); my $attstatus; if( $att_corr->had_error() ){ $attstatus = 120; unlink $att_proc; }else{ $attstatus = 110; my $diff = Util::HEAdas->new('ftdiff') ->params({infile1 => $attitude, infile2 => $att_proc}); $diff->seriousness(0); $diff->run(); my $ndiffs = $diff->parfile()->read('numdiffs'); if($ndiffs==0){ $attstatus = 100; unlink $att_proc; } } $fits->ext(0); $fits->keyword('ATTSTATU', $attstatus, 'Status of corrected attitude files'); for(my $i=0; $i<3; $i++){ $fits->ext($i); $fits->keyword('ATTFLAG', "'100'", 'Attitude origin: 100=sat/spacecraft'); } if( -f $att_proc){ my $pattfits = Util::FITSfile->new($att_proc); for(my $i=0; $i<3; $i++){ $pattfits->ext($i); $pattfits->keyword('ATTFLAG', "'110'", 'Attitude origin: 110=pat/jump corrected'); } } ############################################################ # calculate attitude/orbit calculated filtering quantities # 2016-11-02 moved after attjumpcorr so prefilter can use pat (if available) ############################################################ $self->calculate_att_orb(); } # end of body method ############################################################################# # Extract an attitude file from the ACS packets ############################################################################# sub extract_acs { my $self=shift; my $base=shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); my @list=(); ############################## # get the ACS packets file ############################## my @acs = ( $filename->get("telemetry", "*", "acs", "*"), $filename->get("telemetry", "*", "head2", 485) ); unless(@acs && -f $acs[0]) { $log->entry("No ACS packets available"); return @list; } ########################## # run acs2fits ########################## my $index=0; foreach my $ccsds_file (@acs) { ####################### # get the file name ####################### my $file = "$base.$index"; $log->entry("Creating $file from $ccsds_file"); unlink $file; ############################# # set up the extraction tool ############################# my $teldef = $filename->fetch_cal('alignment'); my $acs2fits=Util::AttTool->new('acs2fits'); $acs2fits->command_line("-infile $ccsds_file", "-alignfile $teldef", "-outfile $file"); ##################################### # run the tool and check for errors ##################################### $acs2fits->run(); if($acs2fits->had_error() ) { ############# # error ############# if(-f $file ) { $log->entry("Removing $file"); unlink $file; } } else { ############# # success ############# push @list, ($file); } } # end of loop over input files return @list; } # end of extract_acs method ############################################################################# # Extract an attitude file from the XRT FITS files ############################################################################# sub extract_xrt { my $self=shift; my $base=shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); my @list=(); ########################################################### # types of XRT HK files with attitude data in them ########################################################### my @modes = ('48a', '4e0', '4e1', '4e2', '4f0', '500', '534', '536', '53a', '540'); ### "48Ah", "4E0h", "4E1h", "4E2h", ### "4F0h", "500h", "534h", "536h", "53Ah", "540h"); my @files=(); my @head_files = ( $filename->get('hk', 'xrt', 'hd', '*') ); push @files, map $_ .= "[1]", @head_files; my $eng_hk = $filename->get('enhk', 'xrt', '', 0); if( -f $eng_hk ){ my $fitsEng = Util::FITSfile->new($eng_hk); my $nhdus = $fitsEng->nhdus(); for(my $iext=1; $iext < $nhdus; $iext++){ $fitsEng->ext($iext); my $name = $fitsEng->keyword('EXTNAME'); if( grep(/^hk${name}x$/, @modes) ){ push @files, $eng_hk . "[$name]"; } } } unless(@files) { $log->entry("No attitude data in the XRT FITS files"); return @list; } ################################### # extract the attitude files ################################### my $index=0; foreach my $xrt_file (@files) { ####################### # get the file name ####################### my $file = "$base.".$index++; $log->entry("Creating $file from $xrt_file"); unlink $file; ############################# # set up the extraction tool ############################# my $teldef = $filename->fetch_cal('alignment'); my $tool=Util::AttTool->new('xrtatt'); $tool->command_line("-infile '$xrt_file'", "-alignfile $teldef", "-outfile $file"); ##################################### # run the tool and check for errors ##################################### $tool->run(); if($tool->had_error() ) { ############# # error ############# if(-f $file ) { $log->entry("Removing $file"); unlink $file; } } else { ############# # success ############# push @list, ($file); } } # end of loop over input files return @list; } # end of extract_xrt method ############################################################################# # create a GTI file from an attitude file which excludes gaps in the # attitude records ############################################################################# sub no_gap_gtis { my $self = shift; my $att = shift; my $gap = shift; my $log =$self->log(); my $filename=$self->filename(); $log->entry("Creating GTI file which excludes attitude ". "record gaps > $gap s"); ############################################# # read the time records in the attitude file ############################################# my $fits = Util::FITSfile->new($att, 1); my @time = $fits->cols("TIME")->table(); unless(@time) { $log->entry("No attitude records"); return ""; } ####################################### # open an ASCII file to hold the data ####################################### my $data = "att_gtis.tmp"; unlink $data; open DATA, ">$data"; my $start=$time[0]; my $i; my $nrows = $fits->nrows(); for($i=0; $i<$nrows; $i++) { $start = $time[$i]; ############################################################### # skip over consecutive rows to find the start of the next gap ############################################################### while ($i < $nrows && $time[$i] - $time[$i-1] < $gap ) { $i++ } #################### # record the GTI #################### print DATA "$start $time[$i-1]\n"; ############################################################### # skip over the gap ############################################################### while ($i < $nrows && $time[$i] - $time[$i-1] >= $gap ) { $i++ } } # end of loop over attitude records close DATA; my $header = "att_gtis_header.tmp"; open HEADER, ">$header"; print HEADER "START 1D\n"; print HEADER "STOP 1D\n"; close HEADER; # my $gti = "att_gtis_fits.tmp"; my $gti = $filename->get('gti', 's', 'ng', 0); my $fcreate = Util::Ftool->new("fcreate"); $fcreate->params({cdfile => $header, datafile => $data, outfile => $gti, headfile => " ", tbltype => "binary", nskip => 0, nrows => 0, history => "yes", morehdr => 0, extname => "GTI", anull => " ", inull => 0, clobber => "yes"}) ->run(); unlink $data; unlink $header; unless($fcreate->had_error() ) { return $gti } else { return "" } } # end of no_gap_gtis; ############################################################################# # Calculate attitude/orbit derived filter quantities ############################################################################# sub calculate_att_orb { my $self=shift; my @columns=@_; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); $log->entry("Calculating attitude/orbit filtering quantities"); ####################################### # make sure we have an attitude file # othewise there's no sense bothering ####################################### my $attitude = $filename->get('attcorr', 'p'); unless (-f $attitude) { $attitude = $filename->get('attitude', 's'); unless (-f $attitude) { $log->entry("No attitude file available"); return; } } $log->entry("Running prefilter with $attitude attitude"); ################################################################## # assemble the mission zero time into a string to feed prefilter ################################################################## my $epoch = $procpar->read("refdate")."T".$procpar->read("reftime"); my $fits = Util::FITSfile->new($attitude, 1); ############################################## # Add a buffer of 60 seconds to start and stop ############################################## my $tstart = $fits->keyword("TSTART") - 60.0; my $tstop = $fits->keyword("TSTOP") + 60.0; my $attorb = $filename->get('attorb', 's'); unlink $attorb; my $cols = 'ALL'; $cols = join(' ', @columns) if @columns; my $orbname = $filename->fetch_orbit(); ##################################### # set up and run prefilter ##################################### my $prefilter = Util::HEAdas->new("prefilter"); $prefilter->params({outname => $attorb, columns => $cols, orbmode => "TLE_TEXT2", orbname => $orbname, attname => $attitude, alignfile => 'CALDB', leapname => $filename->fetch_cal("leapsec"), rigname => $filename->fetch_cal("rigidity"), start => $tstart, end => $tstop, interval => 1.0, # was 30.0 attextrap => '32', origin => 'GSFC', ranom => $jobpar->read("ra"), decnom => $jobpar->read("dec"), missepoch => $epoch, compcols => "TIME", compapplyquat => "no", chatter => 4, clobber => "yes" }); $log->entry("Running ". $prefilter->name() ); $prefilter->run(); } # end of calculate_att_orb ############################################################################# # Calculate attitude/orbit derived filter quantities ############################################################################# sub correct_aberration { my $self=shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); $log->entry('Doing aberration correction of attitude file'); ####################################### # make sure we have an attitude file # othewise there's no sense bothering ####################################### my $attitude = $filename->get('attitude', 's'); unless( -f $attitude ) { $log->entry('No attitude data available'); return; } $log->entry('Preliminary prefilter run'); $self->calculate_att_orb('TIME', 'POSITION', 'VELOCITY'); my $aberration = Util::HEAdas->new('aberrator') ->params({infile => $attitude, orbfile => $filename->get('attorb', 's'), alignfile => 'CALDB' }); $log->entry('Run aberrator'); $aberration->run(); } # end of correct_aberration sub determineAspect { my ($self, $ontarget) = @_; my $problem = 0; my $level = 1; # determine_pointing is found in comN.N/Subs/NominalPointing.pm. # It stores the aspect in the ra, dec, and roll parameters of the # jobpar file. $self->determine_pointing($ontarget, ERROR_LEVEL => $level, PROBLEM_REF => \$problem); my $jobpar = $self->jobpar; my $ra = $jobpar->read('ra'); my $dec = $jobpar->read('dec'); my $roll = $jobpar->read('roll'); my %aspect = ( TAG => 'aspect', RA_deg => $ra, DEC_deg => $dec, ROLL_deg => $roll, ERROR => $problem, ); return \%aspect; } sub checkDatabaseForNominalPointing { my ($self) = @_; my $procpar = $self->procpar; my $jobpar = $self->jobpar; my $log = $self->log; my $sequence = $jobpar->read('sequence'); my $pnt = undef; my $rdbbin = $ENV{RDB_BIN} || '/aps/tools/rdb'; my $rdbtables = $ENV{RDB_TABLES_DIR} || '/aps/db/rdb/tables'; my $rdbpath = $rdbtables . '/sw/nominal_pnt.rdb'; if (not -f $rdbpath) { $log->error(1, "missing nominal pointing table"); } else { my @out = qx($rdbbin/row < $rdbpath sequence eq $sequence | $rdbbin/headchg -del); if (@out == 1) { my $record = $out[0]; chomp($record); $log->entry("nominal pointing override: $record"); my ($sequence, $ra, $dec, $roll) = split(/\t/, $record); my $ok = 0; foreach my $x ($sequence, $ra, $dec, $roll) { if (Scalar::Util::looks_like_number($x)) { ++$ok; } else { $log->error(1, "bad value '$x' for $sequence"); } } if ($ok == 4) { $pnt = { TAG => 'RDB', RA_deg => $ra, DEC_deg => $dec, ROLL_deg => $roll, }; } } else { $log->entry("no nominal pointing in DB for $sequence"); } } return $pnt; } sub formatPNT { my ($p) = @_; my $roll; if (ref($p->{ROLLS_deg})) { my @rolls = map { sprintf('%.1f', $_) } @{ $p->{ROLLS_deg} }; $roll = join('/', @rolls); } else { $roll = sprintf('%.1f', $p->{ROLL_deg}); } my $str = sprintf('RA=%.3f Dec=%.3f Roll=%s [deg]', $p->{RA_deg}, $p->{DEC_deg}, $roll); return $str; } sub d2r { my ($deg) = @_; my $rad = Math::toRadians($deg); return $rad; } sub separation_deg { my ($p1, $p2) = @_; my $u1 = Math::rd2unit(d2r($p1->{RA_deg}), d2r($p1->{DEC_deg})); my $u2 = Math::rd2unit(d2r($p2->{RA_deg}), d2r($p2->{DEC_deg})); my $rad = Math::u3angle($u1, $u2); my $deg = Math::toDegrees($rad); return $deg; } sub storeNominalPointing { my ($self, $pnt) = @_; my $log = $self->log; my $jobpar = $self->jobpar; my $ra = $pnt->{RA_deg} || 0; my $dec = $pnt->{DEC_deg} || 0; my $roll = $pnt->{ROLL_deg} || 0; my ($euler1, $euler2, $euler3) = convertRADecRollToSwiftEuler($ra, $dec, $roll); ############################################# # convert to galactic coordinates ############################################# my $ftcoco = Util::HEAdas->new('ftcoco')->is_script(1)->params({ infile => 'NONE', outfile => '', incoord => 'R', outcoord => 'G', lon1 => $ra, lat1 => $dec, radecsys => 'FK5', chatter => 0, clobber => 'no', history => 'no' })->run(); # Get the output Lii/Bii params if no error my $glon = 0.0; my $glat = 0.0; if ($ftcoco->had_error()) { $log->error( 2, "Error in ftcoco!" ); } else { $glon = $ftcoco->parfile()->read('lon2', 1e-6); $glat = $ftcoco->parfile()->read('lat2', 1e-6); } #################################################### # log the results and write them to the jobpar file #################################################### $log->entry("Nominal spacecraft Euler angles: ". "Phi=$euler1 Theta=$euler2 Psi=$euler3"); $log->entry("Nominal telescope pointing: ". "R.A.=$ra Dec=$dec Roll Angle=$roll"); $log->entry("Galactic Coordinates: ". "Lii=$glon Bii=$glat"); $jobpar->set({ ra => $ra, dec => $dec, roll => $roll, euler1 => $euler1, euler2 => $euler2, euler3 => $euler3, glon => $glon, glat => $glat }); } sub determineRequestedPointing { my ($self) = @_; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; my $senhk = $filename->get('scenhk', 'swift', '', 0); if (not -f $senhk) { $log->error(1, "unable to determine SACSREQ without engineering housekeeping"); return; } my $outfile = 'sacsreq.fits'; my $obsid = $jobpar->read('sequence'); Util::FITSfile->new($senhk, 0)->keyword(OBS_ID => $obsid); my $swacsreq = Util::HEAdas->new('swacsreq') ->is_script(1) ->params({ infile => $senhk, outfile => $outfile, chatter => 3, clobber => 'yes', history => 'no' })->run(); if (not -f $outfile) { $log->error(1, "unable to determine requested pointing from $senhk"); return; } my %keys = Util::FITSfile->new($outfile, 0)->keywords; my ($ra, $dec); my @roll; foreach my $key (keys(%keys)) { next if ($key !~ /^SW/); my $value = $keys{$key}; if ($key eq 'SWRA') { $ra = $value; } elsif ($key eq 'SWDEC') { $dec = $value; } elsif ($key =~ /^SWROLL/) { push(@roll, $value); } } my $sacsreq = undef; if (defined($ra) and defined($dec) and @roll) { $sacsreq = { RA_deg => $ra, DEC_deg => $dec, ROLL_deg => $roll[0], ROLLS_deg => \@roll, }; $log->entry("determined requested pointing: RA=$ra Dec=$dec roll=$roll[0]"); } else { $log->error(1, "unable to determine requested pointing from $senhk: missing RA/Dec/roll keys"); } return $sacsreq; } sub testCelestialToEuler { testOne(325.4736, 2.669, 62.4474); # => 56.8655, 62.47986, -93.01 testOne(299.5535, 35.1857, 307.5312); # => -174.3233, 60.1397, 131.639 testOne(186.44, 16.321, 293.61) # => 89.4374, 67.3954, 107.722 } sub testOne { my @p = @_; print join("\n\t", 'pointing:', @p), "\n"; my @e = convertRADecRollToSwiftEuler(@p); print join("\n\t", 'euler:', @e), "\n"; } sub convertRADecRollToSwiftEuler { my ($ra_d, $dec_d, $roll_d) = @_; $roll_d *= -1; my $d2r = Math::toRadians(1); my $r2d = 1 / $d2r; my @e_d = ($ra_d, 90 - $dec_d, $roll_d + 90); my @e_r = map { $_ * $d2r } @e_d; my $cos_phi = cos($e_r[0]); my $sin_phi = sin($e_r[0]); my $cos_theta = cos($e_r[1]); my $sin_theta = sin($e_r[1]); my $cos_psi = cos($e_r[2]); my $sin_psi = sin($e_r[2]); # celestial unit vectors my @celux = ( $cos_phi*$cos_theta*$cos_psi - $sin_phi*$sin_psi, $sin_phi*$cos_theta*$cos_psi + $cos_phi*$sin_psi, -$sin_theta*$cos_psi ); my @celuy = ( -$cos_phi*$cos_theta*$sin_psi - $sin_phi*$cos_psi, -$sin_phi*$cos_theta*$sin_psi + $cos_phi*$cos_psi, $sin_theta*$sin_psi ); my @celuz = ( $cos_phi*$sin_theta, $sin_phi*$sin_theta, $cos_theta, ); # spacecraft unit vectors my @scx = @celuz; my @scy = @celux; my @scz = @celuy; # spacecraft euler angles my $scz2 = $scz[2]; # range check and circumvent exception if (abs($scz2) > 1) { if (abs($scz2) - 1 > 1e-6) { print "invalid unit vector: scz2=$scz2\n"; } $scz2 = $scz2 < 0 ? -1 : 1; } my @sce_r = ( POSIX::atan2($scz[1], $scz[0]), POSIX::acos($scz2), POSIX::atan2($scy[2], -$scx[2]) ); my @sce_d = map { $r2d * $_ } @sce_r; return @sce_d; } sub deltaRoll_d { my ($x, $y) = @_; if ($x < 0 || $x > 360) { # carp("invalid roll $x"); return -1; } if ($y < 0 || $y > 360) { # carp("invalid roll $x"); return -1; } my $deg = undef; my $q = abs($x - $y); if ($q > 180) { if ($x < $y) { $deg = $x + 360 - $y; } else { $deg = $y + 360 - $x; } } else { $deg = $q; } return $deg; } 1;