This notebook performs functional testing of test load products created starting with
existing OFLS products and processing with Matlab tools 2016_210.
This release introduces dynamical aimpoint offsets that
compensate for temperature-dependent alignment changes of the ACA. The Matlab
tools code in turn depends a new version 0.7 of the chandra_aca package which
provides the core function to compute the temperature-dependent dynamical aimpoint offsets.
Overall test requirements are defined in the Aimpoint Transition Plan for Cycle 18. The functional testing in this notebook consists of the following:
For post-NOV0215 products with a nominal zero-offsets aimpoint table
Check dynamical offsets ACA attitude matches FLIGHT maneuver summary for ACIS observations within ~10 arcsec.
Check that predicted CHIPX / CHIPY matches expectation to within 10 arcsec.
dmcoords to compute predicted CHIPX / CHIPY.For pre-NOV1615 products with an empty zero-offsets aimpoint table
ALL TESTS PASS
In [1]:
# Required imports
from __future__ import division, print_function
import glob
import shelve
from astropy.table import Table
import numpy as np
from Quaternion import Quat
import chandra_aca
from chandra_aca import calc_aca_from_targ, calc_targ_from_aca
import parse_cm
from Ska.engarchive import fetch_sci
from Chandra.Time import DateTime
import Ska.Shell
import functools
import mica.archive.obspar
import Ska.arc5gl
import os
import parse_cm.maneuver
%matplotlib inline
These are nominal flight-like products created using the existing MAR0516 loads.
One expected difference is that Matlab uses the OFLS CHARACTERISTICS file
CHARACTERIS_07JUL16 to back out the OFLS ODB_SI_ALIGN transform.
MAR0516 was planned using CHARACTERIS_28FEB16. There is a (1.77, 0.33) arcsec difference
for both ACIS-S and ACIS-I in the ODB_SI_ALIGN offsets. This 1.8 arcsec offset
is seen in comparisons while in newly generated products there will be no such offset.
In [2]:
PRODUCTS = 'MAR0516O'
TEST_DIR = '/proj/sot/ska/ops/SFE'
FLIGHT_DIR = '/data/mpcrit1/mplogs/2016' # Must match PRODUCTS
In [3]:
# SI_ALIGN from Matlab code
SI_ALIGN = chandra_aca.ODB_SI_ALIGN
SI_ALIGN
Out[3]:
In [4]:
def print_dq(q1, q2):
"""
Print the difference between two quaternions in a nice formatted way.
"""
dq = q1.inv() * q2
dr, dp, dy, _ = np.degrees(dq.q) * 2 * 3600
print('droll={:6.2f}, dpitch={:6.2f}, dyaw={:6.2f} arcsec'.format(dr, dp, dy))
In [5]:
def check_obs(obs):
"""
Check `obs` (which is a row out of the dynamic offsets table) for consistency
between target and aca coordinates given the target and aca offsets and the
SI_ALIGN alignment matrix.
"""
y_off = (obs['target_offset_y'] + obs['aca_offset_y']) / 3600
z_off = (obs['target_offset_z'] + obs['aca_offset_z']) / 3600
q_targ = Quat([obs['target_ra'], obs['target_dec'], obs['target_roll']])
q_aca = Quat([obs['aca_ra'], obs['aca_dec'], obs['aca_roll']])
q_aca_out = calc_aca_from_targ(q_targ, y_off, z_off, SI_ALIGN)
print('{} {:6s} '.format(obs['obsid'], obs['detector']), end='')
print_dq(q_aca, q_aca_out)
In [6]:
def check_obs_vs_manvr(obs, manvr):
"""
Check against attitude from actual flight products (from maneuver summary file)
"""
mf = manvr['final']
q_flight = Quat([mf['q1'], mf['q2'], mf['q3'], mf['q4']])
q_aca = Quat([obs['aca_ra'], obs['aca_dec'], obs['aca_roll']])
print('{} {:6s} chipx={:8.2f} chipy={:8.2f} '
.format(obs['obsid'], obs['detector'], obs['chipx'], obs['chipy']), end='')
print_dq(q_aca, q_flight)
In [7]:
filename = os.path.join(TEST_DIR, PRODUCTS, 'ofls', 'output', '{}_dynamical_offsets.txt'.format(PRODUCTS))
dat = Table.read(filename, format='ascii')
In [8]:
dat[:5]
Out[8]:
This checks that applying the sum of target and ACA offsets along with the nominal SI_ALIGN to the target attitude produces the ACA attitude.
All the ACIS observations have a similar offset due to the CHARACTERISTICS mismatch noted earlier, while the HRC observations show the expected 0.0 offset.
In [9]:
for obs in dat:
check_obs(obs)
In [10]:
filename = glob.glob(os.path.join(TEST_DIR, PRODUCTS, 'ofls', 'mps', 'mm*.sum'))[0]
print('Reading', filename)
mm = parse_cm.maneuver.read_maneuver_summary(filename, structured=True)
mm = {m['final']['id']: m for m in mm} # Turn into a dict
In [11]:
for obs in dat:
check_obs_vs_manvr(obs, mm[obs['obsid'] * 100])
It is assumed that the maneuver summary matches the load product attitudes.
There are three discrepancies below: obsids 18725, 18790, and 18800. All of these are DDT observations that are configured in the OR list to use Cycle 18 aimpoint values, but without changing the target offsets from cycle 17 values. This is an artifact of testing and would not occur in flight planning.
In [12]:
os.path.join(FLIGHT_DIR, PRODUCTS[:-1], 'ofls', 'mps', 'mm*.sum')
Out[12]:
In [13]:
filename = glob.glob(os.path.join(FLIGHT_DIR, PRODUCTS[:-1], 'ofls', 'mps', 'mm*.sum'))[0]
print('Reading', filename)
mm = parse_cm.maneuver.read_maneuver_summary(filename, structured=True)
mm = {m['final']['id']: m for m in mm} # Turn into a dict
In [14]:
for obs in dat:
check_obs_vs_manvr(obs, mm[obs['obsid'] * 100])
dmcoords to compute predicted CHIPX / CHIPY.In the results below there are two discrepancies, obsids 18168 and 18091. These are very hot observations coming right after safe mode recovery. In this case the thermal model is inaccurate and the commanded pointing would have been offset by up to 15 arcsec. Future improvements in thermal modeling could reduce this offset, but it should be understood that pointing accuracy will be degraded in such a situation.
In [15]:
ciaoenv = Ska.Shell.getenv('source /soft/ciao/bin/ciao.sh')
ciaorun = functools.partial(Ska.Shell.bash, env=ciaoenv)
In [16]:
dmcoords_cmd = ['dmcoords', 'none',
'asolfile=none',
'detector="{detector}"',
'fpsys="{fpsys}"',
'opt=cel',
'ra={ra_targ}',
'dec={dec_targ}',
'celfmt=deg',
'ra_nom={ra_nom}',
'dec_nom={dec_nom}',
'roll_nom={roll_nom}',
'ra_asp=")ra_nom"',
'dec_asp=")dec_nom"',
'roll_asp=")roll_nom"',
'sim="{sim_x} 0 {sim_z}"',
'displace="0 {dy} {dz} 0 0 0"',
'verbose=0']
dmcoords_cmd = ' '.join(dmcoords_cmd)
In [17]:
def dmcoords_chipx_chipy(keys, verbose=False):
"""
Get the dmcoords-computed chipx and chipy for given event file
header keyword params. NOTE: the ``dy`` and ``dz`` inputs
to dmcoords are flipped in sign from the ASOL values. Generally the
ASOL DY/DZ are positive and dmcoord input values are negative. This
sign flip is handled *here*, so input to this is ASOL DY/DZ.
:param keys: dict of event file keywords
"""
# See the absolute_pointing_uncertainty notebook in this repo for the
# detailed derivation of this -15.5, 6.0 arcsec offset factor. See the
# cell below for the summary version.
ciaorun('punlearn dmcoords')
fpsys_map = {'HRC-I': 'HI1',
'HRC-S': 'HS2',
'ACIS': 'ACIS'}
keys = {key.lower(): val for key, val in keys.items()}
det = keys['detnam']
keys['detector'] = (det if det.startswith('HRC') else 'ACIS')
keys['dy'] = -keys['dy_avg']
keys['dz'] = -keys['dz_avg']
keys['fpsys'] = fpsys_map[keys['detector']]
cmd = dmcoords_cmd.format(**keys)
ciaorun(cmd)
if verbose:
print(cmd)
return [float(x) for x in ciaorun('pget dmcoords chipx chipy chip_id')]
In [18]:
def get_evt_meta(obsid, detector):
"""
Get event file metadata (FITS keywords) for ``obsid`` and ``detector`` and cache for later use.
Returns a dict of key=value pairs.
"""
evts = shelve.open('event_meta.shelf')
sobsid = str(obsid)
if sobsid not in evts:
det = 'hrc' if detector.startswith('HRC') else 'acis'
arc5gl = Ska.arc5gl.Arc5gl()
arc5gl.sendline('obsid={}'.format(obsid))
arc5gl.sendline('get {}2'.format(det) + '{evt2}')
del arc5gl
files = glob.glob('{}f{}*_evt2.fits.gz'.format(det, obsid))
if len(files) != 1:
raise ValueError('Wrong number of files {}'.format(files))
evt2 = Table.read(files[0])
os.unlink(files[0])
evts[sobsid] = {k.lower(): v for k, v in evt2.meta.items()}
out = evts[sobsid]
evts.close()
return out
In [19]:
def check_predicted_chipxy(obs):
"""
Compare the predicted CHIPX/Y values with planned using observed event file
data on actual ACA alignment.
"""
obsid = obs['obsid']
detector = obs['detector']
try:
evt = get_evt_meta(obsid, detector)
except ValueError as err:
print('Obsid={} detector={}: fail {}'.format(obsid, detector, err))
return
f_chipx, f_chipy, f_chip_id = dmcoords_chipx_chipy(evt)
q_nom_flight = Quat([evt['ra_nom'], evt['dec_nom'], evt['roll_nom']])
q_aca = Quat([obs['aca_ra'], obs['aca_dec'], obs['aca_roll']])
mf = mm[obsid * 100]['final']
q_flight = Quat([mf['q1'], mf['q2'], mf['q3'], mf['q4']])
dq = q_flight.dq(q_aca)
q_nom_test = q_nom_flight * dq
evt_test = dict(evt)
evt_test['ra_nom'] = q_nom_test.ra
evt_test['dec_nom'] = q_nom_test.dec
evt_test['roll_nom'] = q_nom_test.roll
scale = 0.13175 if detector.startswith('HRC') else 0.492
aim_chipx = obs['chipx']
aim_chipy = obs['chipy']
if detector == 'ACIS-S':
aim_chipx += -obs['target_offset_y'] / scale
aim_chipy += -obs['target_offset_z'] / scale + 20.5 / 0.492 * (-190.14 - evt['sim_z'])
elif detector == 'ACIS-I':
aim_chipx += -obs['target_offset_z'] / scale + 20.5 / 0.492 * (-233.59 - evt['sim_z'])
aim_chipy += +obs['target_offset_y'] / scale
chipx, chipy, chip_id = dmcoords_chipx_chipy(evt_test)
print('{} {:6s} aimpoint:{:6.1f},{:6.1f} test:{:6.1f},{:6.1f} '
'flight:{:6.1f},{:6.1f} delta: {:.1f} arcsec'
.format(obsid, detector, aim_chipx, aim_chipy, chipx, chipy, f_chipx, f_chipy,
np.hypot(aim_chipx - chipx, aim_chipy - chipy) * scale))
In [20]:
for obs in dat:
check_predicted_chipxy(obs)
Check that TEST (JUL0415M) and FLIGHT attitudes from maneuver summary match to within 0.1 arcsec.
JUL0415M was constructed with an OR-list zero-offset aimpoint table which exists but has no row entries. This has the effect of telling SAUSAGE to run through the attitude replacement machinery but use 0.0 for the ACA offset y/z values. This should output attitudes that are precisely the same as the FLIGHT attitudes.
Results: pass
In [21]:
PRODUCTS = 'JUL0415M'
FLIGHT_DIR = '/data/mpcrit1/mplogs/2015' # Must match PRODUCTS
In [22]:
# Get FLIGHT maneuver summary
filename = glob.glob(os.path.join(FLIGHT_DIR, PRODUCTS[:-1], 'ofls', 'mps', 'mm*.sum'))[0]
print('Reading', filename)
mmf = parse_cm.maneuver.read_maneuver_summary(filename, structured=True)
mmf = {m['final']['id']: m for m in mmf} # Turn into a dict
In [23]:
# Get TEST maneuver summary
filename = glob.glob(os.path.join(TEST_DIR, PRODUCTS, 'ofls', 'mps', 'mm*.sum'))[0]
print('Reading', filename)
mmt = parse_cm.maneuver.read_maneuver_summary(filename, structured=True)
mmt = {m['final']['id']: m for m in mmt} # Turn into a dict
In [24]:
# Make sure set of obsids are the same
set(mmf) == set(mmt)
Out[24]:
In [25]:
# Now do the actual attitude comparison
for trace_id, mf in mmf.items():
mt = mmt[trace_id]['final']
mf = mf['final']
qt = Quat([mt['q1'], mt['q2'], mt['q3'], mt['q4']])
qf = Quat([mf['q1'], mf['q2'], mf['q3'], mf['q4']])
print(trace_id, ' ', end='')
print_dq(qt, qf)
In [ ]: