# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Function(s) related to ACA alignment.
In particular compute the dynamical pointing offset required to achieve
a desired zero-offset target aimpoint.
A key element of this module is the fitting analysis here:
https://github.com/sot/aimpoint_mon/blob/master/fit_aimpoint_drift-2018-11.ipynb
"""
import functools
import json
import os
from pathlib import Path
import numpy as np
from astropy.table import Table
from astropy.utils.data import download_file
from Chandra.Time import DateTime
from cxotime import CxoTime, CxoTimeLike
from ska_helpers import chandra_models
from ska_helpers.utils import LazyDict
def load_drift_pars():
pars_txt, info = chandra_models.get_data(
Path("chandra_models") / "aca_drift" / "aca_drift_model.json"
)
pars = json.loads(pars_txt)
pars["info"] = info
return pars
DRIFT_PARS = LazyDict(load_drift_pars)
# Define transform from aspect solution DY, DZ (mm) to CHIPX, CHIPY for
# nominal target position. This uses the transform:
# CHIPX, CHIPY = c0 + cyz * [DY, DZ] (dot product)
# These values come from the aimpoint_mon/asol_to_chip_transforms notebook:
# https://github.com/sot/aimpoint_mon/blob/cebe5010/asol_to_chip_transforms.ipynb
ASOL_TO_CHIP = {
("ACIS-I", 0): {
"c0": [1100.806, 1110.299],
"cyz": [[0.001, 41.742], [-41.742, 0.105]],
},
("ACIS-I", 1): {
"c0": [-76.208, 962.413],
"cyz": [[0.001, -41.742], [41.741, 0.104]],
},
("ACIS-I", 2): {
"c0": [54.468, 1110.197],
"cyz": [[-0.001, 41.742], [-41.742, -0.105]],
},
("ACIS-I", 3): {
"c0": [970.795, 961.857],
"cyz": [[-0.001, -41.742], [41.742, -0.104]],
},
("ACIS-S", 4): {
"c0": [3382.364, 520.024],
"cyz": [[-41.695, -0.021], [0.021, -41.689]],
},
("ACIS-S", 5): {
"c0": [2338.96, 519.15],
"cyz": [[-41.691, -0.031], [0.037, -41.689]],
},
("ACIS-S", 6): {
"c0": [1296.05, 519.532],
"cyz": [[-41.69, -0.031], [0.036, -41.689]],
},
("ACIS-S", 7): {
"c0": [251.607, 519.818],
"cyz": [[-41.689, -0.022], [0.022, -41.689]],
},
("ACIS-S", 8): {
"c0": [-790.204, 519.784],
"cyz": [[-41.69, -0.02], [0.02, -41.689]],
},
("ACIS-S", 9): {
"c0": [-1832.423, 519.546],
"cyz": [[-41.693, -0.017], [0.017, -41.689]],
},
("HRC-I", 0): {
"c0": [7635.787, 7697.627],
"cyz": [[-109.98, 109.98], [109.98, 109.98]],
},
("HRC-S", 1): {
"c0": [2196.931, 25300.796],
"cyz": [[0.188, -155.535], [155.584, 0.19]],
},
("HRC-S", 2): {
"c0": [2197.753, 8842.306],
"cyz": [[0.2, -155.535], [155.535, 0.184]],
},
("HRC-S", 3): {
"c0": [2196.638, -7615.26],
"cyz": [[0.2, -155.535], [155.571, 0.184]],
},
}
SIM_MM_TO_ARCSEC = 20.493
[docs]
class AcaDriftModel(object):
"""
AcaDriftModel class
Define a drift model for aspect solution SIM DY/DZ values as a function of
time and ACA CCD temperature. This expresses the model which is defined
and fitted in the fit_aimpoint_drift notebook in this repo.
"""
def __init__(self, scale, offset, trend, jumps, year0):
self.scale = scale
self.offset = offset
self.trend = trend
self.jumps = jumps
self.year0 = year0
[docs]
def calc(self, times, t_ccd):
"""
Calculate the drift model
Calculate the drift model for aspect solution SIM DY/DZ values for input
``times`` and ``t_ccd``. The two arrays are broadcasted to match.
The returned drifts are in arcsec and provide the expected aspect solution
SIM DY or DZ values in arcsec. This can be converted to a drift in mm
(corresponding to units in an ASOL file) via the scale factor 20.493 arcsec/mm.
Parameters
----------
times
array of times (CXC secs)
t_ccd
CCD temperatures (degC)
Returns
-------
array of ASOL SIM DY/DZ (arcsec)
"""
# The drift model is calibrated assuming t_ccd is in degF, but we want inputs
# in degC, so convert at this point.
t_ccd_degF = t_ccd * 1.8 + 32.0
times, t_ccd_degF = np.broadcast_arrays(times, t_ccd_degF)
is_scalar = times.ndim == 0 and t_ccd_degF.ndim == 0
times = DateTime(np.atleast_1d(times)).secs
t_ccd_degF = np.atleast_1d(t_ccd_degF)
if times.shape != t_ccd_degF.shape:
raise ValueError("times and t_ccd args must match in shape")
if np.any(np.diff(times) < 0):
raise ValueError("times arg must be monotonically increasing")
if times[0] < DateTime("2012:001:12:00:00").secs:
raise ValueError("model is not applicable before 2012")
# Years from model `year0`
dyears = DateTime(times, format="secs").frac_year - self.year0
# Raw offsets without jumps
out = (t_ccd_degF - self.offset) * self.scale + dyears * self.trend
# Put in the step function jumps
for jump_date, jump in self.jumps:
jump_idx = np.searchsorted(times, DateTime(jump_date).secs)
out[jump_idx:] += jump
return out[0] if is_scalar else out
[docs]
def get_fid_offset(time: CxoTimeLike, t_ccd: float) -> tuple:
"""
Compute the fid light offset values for a given time and temperature.
The ``time`` and ``t_ccd`` inputs can be either scalars or arrays.
Parameters
----------
time : CxoTimeLike format
Time for offset calculation.
t_ccd : float
ACA CCD temperature in degrees Celsius.
Returns
-------
tuple
A tuple containing the y-angle and z-angle offsets (in arcseconds) to apply
additively to the nominal (FEB07) fid positions.
Notes
-----
The apparent fid light positions change in accordance with the ACA alignment drift as a
function of time and temperature. This is captured in the ACA aimpoint drift model. This
function uses that model to provide the offsets in y-angle and z-angle (arcsec) to apply
additively to the nominal fid positions.
The y_offset and z_offset values in this function were calibrated using the
2022-11 aimpoint drift model and the FEB07 fid characteristics.
See https://github.com/sot/fid_drift_mon/blob/master/fid_offset_coeff.ipynb
"""
# Clip the time to the minimum time in the drift model
time = CxoTime(time).secs.clip(CxoTime("2012:001:12:00:00.000").secs, None)
# Define model instances using calibrated parameters
drift_y = AcaDriftModel(**DRIFT_PARS["dy"])
drift_z = AcaDriftModel(**DRIFT_PARS["dz"])
# Compute the predicted asol DY/DZ based on time and ACA CCD temperature
# via the predictive model calibrated in the fit_aimpoint_drift notebook
# in this repo. And flip the signs.
dy_pred = -1.0 * drift_y.calc(time, t_ccd)
dz_pred = -1.0 * drift_z.calc(time, t_ccd)
# Apply internal offset that places the fid lights at ~zero position
# offset during the 2022:094 to 2023:044.
y_offset = 19.6
z_offset = 20.1
return dy_pred + y_offset, dz_pred + z_offset
[docs]
def get_aca_offsets(detector, chip_id, chipx, chipy, time, t_ccd):
"""
Compute the dynamical ACA offset values for the provided inputs.
The ``time`` and ``t_ccd`` inputs can be either scalars or arrays.
Parameters
----------
detector
one of ACIS-I, ACIS-S, HRC-I, HRC-S
chipx
zero-offset aimpoint CHIPX
chipy
zero-offset aimpoint CHIPY
chip_id
zero-offset aimpoint CHIP ID
time
time(s) of observation (any Chandra.Time compatible format)
t_ccd
ACA CCD temperature(s) (degC)
Returns
-------
aca_offset_y, aca_offset_z (arcsec)
"""
# Define model instances using calibrated parameters
drift_y = AcaDriftModel(**DRIFT_PARS["dy"])
drift_z = AcaDriftModel(**DRIFT_PARS["dz"])
try:
asol_to_chip = ASOL_TO_CHIP[detector, chip_id]
except KeyError:
raise KeyError(
"Detector and chip combination {} not in allow values: {}".format(
(detector, chip_id), sorted(ASOL_TO_CHIP.keys())
)
) from None
# Compute the asol DY/DZ that would be required for the aimpoint to be
# exactly at the desired CHIPX/Y values. Uses the geometrical transform
# computed via dmcoords in the asol_to_chip_transform notebook in this repo.
# CHIPX, CHIPY = c0 + cyz * [DY, DZ] (dot product)
chip_xy = np.array([chipx, chipy])
cyz_inv = np.linalg.inv(asol_to_chip["cyz"])
dy_chip, dz_chip = (
cyz_inv.dot(chip_xy - asol_to_chip["c0"]) * SIM_MM_TO_ARCSEC
) # arcsec
# Compute the predicted asol DY/DZ based on time and ACA CCD temperature
# via the predictive model calibrated in the fit_aimpoint_drift notebook
# in this repo.
dy_pred = drift_y.calc(time, t_ccd)
dz_pred = drift_z.calc(time, t_ccd)
# The difference is the dynamic ACA offset that must be applied to the attitude.
# This has the same sign convention as the user-supplied TARGET OFFSET in the
# ObsCat / OR-list.
ddy = dy_chip - dy_pred
ddz = dz_chip - dz_pred
return ddy, ddz
[docs]
@functools.lru_cache
def get_default_zero_offset_table():
"""
Get official SOT MP zero offset aimpoint table.
First try ``/data/mpcrit1/aimpoint_table/zero_offset_aimpoints.txt``.
If that is not available use:
https://cxc.harvard.edu/mta/ASPECT/drift/zero_offset_aimpoints.txt.
The web version is updated weekly on Sunday via a Ska cron job.
Note the definitive source of this file is:
https://icxc.harvard.edu/mp/html/aimpoint_table/zero_offset_aimpoints.txt.
Returns
-------
zero offset aimpoint table as astropy.Table
"""
try:
path = (
Path(os.environ["SKA"])
/ "data"
/ "mpcrit1"
/ "aimpoint_table"
/ "zero_offset_aimpoints.txt"
)
out = Table.read(str(path), format="ascii")
except FileNotFoundError:
url = "https://cxc.harvard.edu/mta/ASPECT/drift/zero_offset_aimpoints.txt"
path = download_file(url, show_progress=False, timeout=10)
out = Table.read(path, format="ascii")
return out
[docs]
def get_target_aimpoint(date, cycle, detector, too=False, zero_offset_table=None):
"""
Given date, proposal cycle, and detector, return aimpoint chipx, chipy, chip_id
Parameters
----------
date
observation date
cycle
proposal cycle of observation
detector
target detector
too
boolean. If target is TOO use current cycle not proposal cycle.
zero_offset_able : table (astropy or numpy) of zero offset aimpoint table
defaults to official SOT MP version if not supplied.
Returns
-------
tuple of chipx, chipy, chip_id
"""
if zero_offset_table is None:
zero_offset_table = get_default_zero_offset_table()
zero_offset_table.sort(["date_effective", "cycle_effective"])
date = DateTime(date).iso[:10]
# Entries for this detector before the 'date' given
ok = (zero_offset_table["detector"] == detector) & (
zero_offset_table["date_effective"] <= date
)
# If a regular observation, the entry must also be before or equal to proposal cycle
if not too:
ok = ok & (zero_offset_table["cycle_effective"] <= cycle)
filtered_table = zero_offset_table[ok]
# Return the desired keys in the most recent [-1] row that matches
return tuple(filtered_table[["chipx", "chipy", "chip_id"]][-1])