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 my $file ($acslist->files() ) {
$log->entry("Deleting $file");
unlink $file;
}
}
} 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++;
}
}
my $attitude = $filename->get('attitude', 's');
my $tstart = 0;
# check for attitude override;
my $att_force = undef;
{
my $attdir = $procpar->read('att_force_dir');
my $seq = $jobpar->read('sequence');
my $path = "${attdir}/sw${seq}sat.fits";
if (-f $path) {
$att_force = $path;
$log->entry("Found attitude force file: $path");
}
}
#######################################
# 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 if not $att_force;
} else {
$log->entry("Merging ". join(" ", $attlist->files()) );
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 my $file ($attlist->files() ) {
$log->entry("Deleting $file");
unlink $file;
}
}
$tstart = $self->normalizeAttitude($attitude);
}
if ($att_force) {
# override attitude- append raw attitude to this file
my $tmpfile = 'attitude.force';
Util::FITSfile->new($att_force)
->copy($tmpfile);
$tstart = $self->normalizeAttitude($tmpfile);
if (-f $attitude) {
$log->entry('Appending raw attitude to forced attitude file');
Util::FITSfile->new($attitude, 'ATTITUDE', q([col EXTNAME='RAW_ATTITUDE']))
->append_to($tmpfile);
Util::FITSfile->new($attitude, 'ACS_DATA', q([col EXTNAME='RAW_ACS_DATA']))
->append_to($tmpfile);
}
rename($tmpfile, $attitude);
# An E2 error is introduced to ensure operator
# intervention for the FFA processing
$log->error(2, "Overrode attitude with $att_force");
}
return if not $tstart;
my $fits = Util::FITSfile->new($attitude, 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 $snapshots = $filename->get('gti', 's', 'ss');
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 $pointing0 = $filename->get('gti', 's', 'po0', 0);
my $not_pointing0 = $filename->get('gti', 's', 'np0', 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);
my $gtitmp = 'gti.tmp';
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 => $pointing0,
prefr => 0,
postfr => 0,
expr => 'FLAGS == b11x0xxxx'})
->run();
if (-f $snapshots) {
# this for loop is here in case we desire to mask other GTIs (e.g., settling)
foreach my $gtifile ($pointing) {
next if not -f $gtifile;
unlink $gtitmp if -f $gtitmp;
$log->entry("masking $gtifile with $snapshots");
Util::Ftool->new('mgtime')
->params({ingtis => "$gtifile,$snapshots",
outgti => $gtitmp,
merge => 'AND'})
->run();
if (-f $gtitmp) {
$log->entry("comparing original and corrected $gtifile");
my $diff = Util::HEAdas->new('ftdiff')
->params({infile1 => $gtifile,
infile2 => $gtitmp});
$diff->seriousness(0);
$diff->run();
unlink($gtifile);
rename($gtitmp, $gtifile);
}
}
}
{
# not pointing = inverse(POINTING) AND extended(SNAPSHOTS)
my $delta = 200;
my $prestart = $jobpar->{TIMELIST}{start} - $delta;
my $poststop = $jobpar->{TIMELIST}{stop} + $delta;
if (-f $pointing) {
my @np;
my $pofits = Util::FITSfile->new($pointing);
my @start = $pofits->cols('START')->table();
my @stop = $pofits->cols('STOP')->table();
for (my $i = 0; $i < @start; ++$i) {
if (not @np) {
push(@np, $prestart, $start[$i]);
}
else {
push(@np, $stop[$i-1], $start[$i]);
}
}
push(@np, $stop[-1], $poststop);
Util::GTIfile->new($not_pointing, @np);
}
else {
$log->entry("no pointing GTI so creating $not_pointing covering extended snapshots");
Util::GTIfile->new($not_pointing, $prestart, $poststop);
}
}
{
# the old not pointing GTI:
$maketime->params({outfile => $not_pointing0,
prefr => 1,
postfr => 1,
expr => 'FLAGS != b11xxxxxx || FLAGS == bxxx1xxxx'})
->run();
if (-f $not_pointing and -f $not_pointing0) {
$log->entry("comparing old and new $not_pointing");
my $diff = Util::HEAdas->new('ftdiff')
->params({infile1 => $not_pointing0,
infile2 => $not_pointing});
$diff->seriousness(0);
$diff->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, $pointing0, $not_pointing0){
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');
}
if (not $att_force) {
#################################################################
# Only do aberration correction before 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
sub normalizeAttitude
{
my ($self, $attitude) = @_;
my $log = $self->log;
##########################################
# TIME sort the merged attitude file and
# remove overdetermined values
##########################################
$log->entry("Normalizing $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 $filename = $self->filename;
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');
}
}
return $tstart;
}
#############################################################################
# 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;