This notebook documents and computes fit coefficients for a simple model that
gives the relative ACA alignment as a linear function of the ACA CCD temperature.
It also includes validation of the implementation of the new model in the
chandra_aca.drift package.
This is based on the origin fit_aimpoint_drift notebook but updated through
2018:310 to include data from after the 2018:283 safe mode normal sun dwell.
Two jumps were added to the model, one corresponding to the 2017:066 NSM
and one from the 2018:283 safe mode.
The ACA alignment is measured accurately for each science observation via the apparent positions of the fid lights. These are referred to by their CXC aspect solution designation as the SIM DY and DZ offsets. This is actually a misnomer based on the pre-launch understanding of what physical mechanism would generate such offsets. We now know via HRMA optical axis measurements that a temperature-dependent change in the ACA boresight alignment is responsible. The HRMA to SIM alignment is quite stable.
The ACA alignment relates directly to the X-ray detector aimpoint that is used in observation planning and analysis. With this model it will be possible to improve the aimpoint accuracy by introducing a dynamic pointing offset based on the predicted ACA CCD temperature for each observation.
The model is
DY/Z = (t_ccd - offset) * scale + (year - 2016.0) * trend + JUMPSwhere
t_ccd : ACA CCD temperature (degF)
scale : scaling in arcsec / degF
offset : ACA CCD temperature corresponding to DY/Z = 0.0 arcsec
trend : Trend in DY/Z (arcsec / year)
year : decimal year
jumpYYYYDDD : step function from 0.0 to jumpYYYYDDD (arcsec) for date > YYYY:DDDThe jumps are persistent step function changes in alignment that have been observed following extended dwells at normal sun where the ACA gets substantially hotter than during normal operations. The exact mechanism is not understood, but could be due to a non-linear stiction release of a stress point that impacts alignment.
Note that the ACA alignment has a direct linear correlation to the ACA housing temperature (AACH1T). However, in this model we use the ACA CCD temperature as the model dependent variable because it is linearly related to housing temperature (AACCDPT = m * AACH1T + b) as long as the TEC is at max drive current. Since there is already an existing Xija model to predict ACA CCD temperature this reduces duplication.
This model was fitted to data from 2012:180 to 2018:310 using Sherpa. The key fit results are:
DY
-----
scale = 2.1 arcsec / degF = 3.9 arcsec / degC
trend = -1.11 arcsec / year
jumps ~ -2 to -13 arcsec
model error = +/- 1.9 arcsec (1st to 99th percentile range)
DZ
-----
scale = 1.0 arcsec / degF = 1.8 arcsec / degC
trend = -0.16 arcsec / year
jumps ~ -0.4 to -6.1 arcsec
model error = +/- 2.8 arcsec (1st to 99th percentile range)The model accuracy will be degraded somewhat when ACA CCD temperature is taken from a predictive Xija model instead of from telemetry.
This notebook lives in the aimpoint_mon project repository
In [1]:
import re
import tables
import matplotlib.pyplot as plt
import numpy as np
from astropy.time import Time
from astropy.table import Table
import Ska.engarchive.fetch_eng as fetch
from Ska.engarchive import fetch_sci
from Chandra.Time import DateTime
from Ska.Numpy import interpolate
from kadi import events
from sherpa import ui
from Ska.Matplotlib import plot_cxctime
In [2]:
%matplotlib inline
In [3]:
SIM_MM_TO_ARCSEC = 20.493
In [4]:
# Discrete jumps after 2012:001. Note also jumps at:
# '2008:293', # IU-reset
# '2010:151', # IU-reset
# '2011:190', # Safe mode
JUMPS = ['2015:006', # IU-reset
'2015:265', # Safe mode 6
'2016:064', # Safe mode 7
'2017:066', # NSM
'2018:285', # Safe mode 8
]
In [5]:
ltt_bads = events.ltt_bads(pad=(0, 200000))
normal_suns = events.normal_suns(pad=(0, 100000))
safe_suns = events.safe_suns(pad=(0, 86400 * 7))
In [6]:
# Aspect camera CCD temperature trend since 2010
t_ccd = fetch.Msid('aacccdpt', start='2010:001', stat='5min')
t_ccd.remove_intervals(ltt_bads | normal_suns | safe_suns)
plt.figure(figsize=(12, 4.5))
t_ccd.plot()
plt.ylabel('T_ccd (degF)')
plt.title('ACA CCD temperature')
plt.ylim(None, 20)
plt.grid()
In [7]:
# Get aspect solution DY and DZ (apparent SIM offsets via fid light positions)
# which are sampled at 1 ksec intervals and updated daily.
if 'adat' not in globals():
h5 = tables.open_file('/proj/sot/ska/data/aimpoint_mon/aimpoint_asol_values.h5')
adat = h5.root.data[:]
h5.close()
adat.sort(order=['time'])
# Filter bad data when asol DY and DZ are both exactly 0.0 (doesn't happen normally)
bad = (adat['dy'] == 0.0) & (adat['dz'] == 0.0)
adat = adat[~bad]
In [8]:
class AcaDriftModel(object):
"""
Class to encapsulate necessary data and compute the model of ACA
alignment drift. The object created from this class is called
by Sherpa as a function during fitting. This gets directed to
the __call__() method.
"""
YEAR0 = 2016.0 # Reference year for linear offset
def __init__(self, adat, start='2012:001', stop=None):
"""
adat is the raw data array containing aspect solution data
sampled at 1 ksec intervals.
"""
# Get the ACA CCD temperature telemetry
t_ccd = fetch.Msid('aacccdpt', stat='5min', start=start, stop=stop)
# Slice the ASOL data corresponding to available ACA CCD temps
i0, i1 = np.searchsorted(adat['time'], [t_ccd.times[0], t_ccd.times[-1]])
self.asol = adat[i0:i1].copy()
# Convert from mm to arcsec for convenience
self.asol['dy'] *= SIM_MM_TO_ARCSEC
self.asol['dz'] *= SIM_MM_TO_ARCSEC
self.times = self.asol['time']
self.years = Time(self.times, format='cxcsec').decimalyear
self.years_0 = self.years - self.YEAR0
# Resample CCD temp. data to the 1 ksec ASOL time stamps
self.t_ccd = interpolate(t_ccd.vals, t_ccd.times, self.asol['time'], method='linear')
# Get indices corresponding to jump times for later model computation
self.jump_times = Time(JUMPS).cxcsec
self.jump_idxs = np.searchsorted(self.times, self.jump_times)
def __call__(self, pars, years=None, t_ccd=None):
"""
Calculate model prediction for DY or DZ. Params are:
scale : scaling in arcsec / degF
offset : ACA CCD temperature corresponding to DY/Z = 0.0 arcsec
trend : Trend in DY/Z (arcsec / year)
jumpYYYYDDD : discrete jump in arcsec at date YYYY:DDD
"""
# Sherpa passes the parameters as a list
scale, offset, trend = pars[0:3]
jumps = pars[3:]
# Allow for passing in a different value for ACA CCD temperature
if t_ccd is None:
t_ccd = self.t_ccd
# Compute linear part of model
out = (t_ccd - offset) * scale + self.years_0 * trend
# Put in the step function jumps
for jump_idx, jump in zip(self.jump_idxs, jumps):
if jump_idx > 10 and jump_idx < len(out) - 10:
out[jump_idx:] += jump
return out
In [9]:
def fit_aimpoint_aca_temp(axis='dy', start='2012:180', stop=None):
"""
Use Sherpa to fit the model parameters
"""
# Create the object used to define the Sherpa user model, then
# load as a model and create parameters
aca_drift = AcaDriftModel(adat, start, stop)
ui.load_user_model(aca_drift, 'aca_drift_model')
parnames = ['scale', 'offset', 'trend']
parnames += ['jump{}'.format(re.sub(':', '', x)) for x in JUMPS]
ui.add_user_pars('aca_drift_model', parnames)
# Sherpa automatically puts 'aca_drift_model' into globals, but
# make this explicit so code linters don't complain.
aca_drift_model = globals()['aca_drift_model']
# Get the DY or DZ values and load as Sherpa data
dyz = aca_drift.asol[axis]
ui.load_arrays(1, aca_drift.years, dyz)
# Set the model and fit using Simplex (Nelder-Mead) minimization
ui.set_model(1, aca_drift_model)
ui.set_method('simplex')
ui.fit(1)
return aca_drift, ui.get_fit_results()
In [10]:
def plot_aimpoint_drift(axis, aca_drift, fit_results, start='2010:001', stop=None, plot_t_ccd=False):
"""
Plot our results
"""
y_start = DateTime(start).frac_year
y_stop = DateTime(stop).frac_year
years = aca_drift.years
ok = (years > y_start) & (years < y_stop)
years = aca_drift.years[ok]
times = aca_drift.times[ok]
# Call model directly with best-fit parameters to get model values
dyz_fit = aca_drift(fit_results.parvals)[ok]
# DY or DZ values from aspect solution
dyz = aca_drift.asol[axis][ok]
dyz_resid = dyz - dyz_fit
if plot_t_ccd:
plt.figure(figsize=(12, 4.5))
plt.subplot(1, 2, 1)
plot_cxctime(times, dyz, label='Data')
plot_cxctime(times, dyz_fit, 'r-', alpha=0.5, label='Fit')
plot_cxctime(times, dyz_resid, 'r-', label='Residual')
plt.title('Fit aspect solution {} to scaled ACA CCD temperature'
.format(axis.upper()))
plt.ylabel('{} (arcsec)'.format(axis.upper()))
plt.grid()
plt.legend(loc='upper left', framealpha=1.0)
if plot_t_ccd:
dat = fetch_sci.Msid('aacccdpt', start, stop, stat='5min')
plt.subplot(1, 2, 2)
dat.plot()
plt.grid()
plt.ylabel('AACCCDPT (degC)')
if isinstance(plot_t_ccd, tuple):
plt.ylim(*plot_t_ccd)
std = dyz_resid.std()
p1, p99 = np.percentile(dyz_resid, [1, 99])
print('Fit residual stddev = {:.2f} arcsec'.format(std))
print('Fit residual 99th - 1st percentile = {:.2f}'.format(p99 - p1))
In [11]:
aca_drift_dy, fit_dy = fit_aimpoint_aca_temp('dy')
In [12]:
plot_aimpoint_drift('dy', aca_drift_dy, fit_dy)
In [13]:
start = '2018:260'
stop = '2018:310'
plot_aimpoint_drift('dy', aca_drift_dy, fit_dy, start=start, stop=stop, plot_t_ccd=(-16, -8))
Temperatures just after the safe mode gap are cool, around -12.5 C. In this plot the offset
gets more positive with lower temperature (i.e. anti-correlated with T_ccd, about -4 arcsec / degC).
Implications:
So effectively right now we have a new constraint that the ACA CCD needs to be warmer than -13.5 C. Probably not a problem in practice.
This just shows in raw form the 13 arcsec jump. The sign is the same, so increasing temperature means more negative Y offset. Focus on the downwarn spikesaround Y offset = -46 for the blue model, looking before and after safe mode. The raw fid locations have shifted by about 13 arcsec.
Y offset around -22 arcsec.Y offset is around -8 arcsec.
In [14]:
dyz_fit = aca_drift_dy(fit_dy.parvals, t_ccd=14) # degF = -10 C
plot_cxctime(aca_drift_dy.times, dyz_fit)
plt.title('DY drift model assuming constant ACA temperature')
plt.grid();
In [15]:
aca_drift_dz, fit_dz = fit_aimpoint_aca_temp('dz')
In [16]:
plot_aimpoint_drift('dz', aca_drift_dz, fit_dz)
In [17]:
start = '2018:260'
stop = '2018:310'
plot_aimpoint_drift('dz', aca_drift_dz, fit_dz, start=start, stop=stop, plot_t_ccd=(-16, -8))
In [18]:
dyz_fit = aca_drift_dz(fit_dz.parvals, t_ccd=14) # degF = -10 C
plot_cxctime(aca_drift_dz.times, dyz_fit)
plt.title('DZ drift model assuming constant ACA temperature')
plt.grid();
Compare the actual flight aca_offset_y/z from the *_dynamical_offsets.txt files
to predictions with the new chandra_aca.drift module.
A key point is to use the observed mean T_ccd with the new model to be able to reproduce
the observed aimpoint shift of about 8 arcsec. The jump was 13 arcsec but we did not
see that directly because of the ~1.4 C error in the temperatures being used to
predict the aimpoint offset.
In [19]:
text = """
obsid detector chipx chipy chip_id aca_offset_y aca_offset_z mean_t_ccd mean_date
----- -------- ------- ------- ------- ------------ ------------ ---------- ---------------------
21152 ACIS-S 210.0 520.0 7 -0.9 -22.67 -11.72 2018:307:18:07:54.816
20332 ACIS-I 970.0 975.0 3 -14.27 -21.89 -11.88 2018:308:04:03:46.816
21718 HRC-I 7590.0 7745.0 0 -13.39 -22.8 -11.53 2018:313:03:14:10.816
21955 HRC-S 2195.0 8915.0 2 -12.50 -22.57 -11.53 2018:305:16:28:34.816
"""
obss = Table.read(text, format='ascii.fixed_width_two_line')
In [20]:
import sys
import os
sys.path.insert(0, os.path.join(os.environ['HOME'], 'git', 'chandra_aca'))
import chandra_aca
from chandra_aca import drift
from kadi import events
In [21]:
chandra_aca.test(get_version=True)
Out[21]:
In [22]:
for obs in obss:
dwell = events.dwells.filter(obsid=21152)[0]
t_ccd = fetch_sci.Msid('aacccdpt', dwell.start, dwell.stop, stat='5min')
mean_t_ccd = np.mean(t_ccd.vals)
offsets = drift.get_aca_offsets(obs['detector'], chip_id=obs['chip_id'],
chipx=obs['chipx'], chipy=obs['chipy'],
time=obs['mean_date'], t_ccd=mean_t_ccd)
print(obs)
print('T_ccd:', mean_t_ccd, ' Delta offsets Y Z:',
'%.2f' % (obs['aca_offset_y'] - offsets[0]),
'%.2f' % (obs['aca_offset_z'] - offsets[1]))
print()
In [23]:
from chandra_aca.tests.test_all import simple_test_aca_drift
In [24]:
dy, dz, times = simple_test_aca_drift()
In [25]:
plt.figure(figsize=(12, 4.5))
plt.subplot(1, 2, 1)
dy_fit = aca_drift_dy(fit_dy.parvals, t_ccd=14) # degF = -10 C
plot_cxctime(aca_drift_dy.times, dy_fit)
plt.title('DY drift model assuming constant ACA temperature')
plt.grid();
plt.subplot(1, 2, 2)
plot_cxctime(times, dy);
plt.grid()
plt.ylabel('DY (arcsec)');
plt.title('DY drift model from chandra_aca');
In [26]:
plt.figure(figsize=(12, 4.5))
plt.subplot(1, 2, 1)
dz_fit = aca_drift_dz(fit_dz.parvals, t_ccd=14) # degF = -10 C
plot_cxctime(aca_drift_dz.times, dz_fit)
plt.title('DZ drift model assuming constant ACA temperature')
plt.grid();
plt.subplot(1, 2, 2)
plot_cxctime(times, dz);
plt.grid()
plt.ylabel('DZ (arcsec)');
plt.title('DZ drift model from chandra_aca');