RCAL-853: Apply velocity aberration correction to the WFI WCS (#1354)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

changes/1354.assign_wcs.rst

@@ -0,0 +1 @@
+Apply velocity aberration correction to the WFI WCS

docs/roman/assign_wcs/main.rst

@@ -89,6 +89,16 @@ coordinate system using a package called `PySIAF <https://github.com/spacetelesc
   >>> ap = siaf['WFI01_FULL']
   >>> V2_REF, V3_REF = ap.get_reference_point('tel')
 
+Corrections Due to Spacecraft Motion
+------------------------------------
+
+The WCS transforms contain two corrections due to motion of the observatory.
+
+Absolute velocity aberration is calculated onboard when acquiring the guide star, but
+differential velocity aberration effects are calculated during the ``assign_wcs`` step.
+This introduces corrections in the conversion from sky coordinates to observatory
+V2/V3 coordinates, and is stored in the WCS under the ``v2v3vacorr`` frame.
+
 
 .. rubric:: Footnotes
 

romancal/assign_wcs/assign_wcs_step.py

@@ -96,15 +96,29 @@ def load_wcs(input_model, reference_files=None):
         axes_names=("v2", "v3"),
         unit=(u.arcsec, u.arcsec),
     )
+    v2v3vacorr = cf.Frame2D(
+        name="v2v3vacorr",
+        axes_order=(0, 1),
+        axes_names=("v2", "v3"),
+        unit=(u.arcsec, u.arcsec),
+    )
     world = cf.CelestialFrame(reference_frame=coord.ICRS(), name="world")
 
     # Transforms between frames
     distortion = wfi_distortion(output_model, reference_files)
     tel2sky = pointing.v23tosky(output_model)
 
+    # Compute differential velocity aberration (DVA) correction:
+    va_corr = pointing.dva_corr_model(
+        va_scale=input_model.meta.velocity_aberration.scale_factor,
+        v2_ref=input_model.meta.wcsinfo.v2_ref,
+        v3_ref=input_model.meta.wcsinfo.v3_ref,
+    )
+
     pipeline = [
         Step(detector, distortion),
-        Step(v2v3, tel2sky),
+        Step(v2v3, va_corr),
+        Step(v2v3vacorr, tel2sky),
         Step(world, None),
     ]
     wcs = WCS(pipeline)

romancal/assign_wcs/pointing.py

@@ -1,7 +1,12 @@
+import logging
+
 import numpy as np
-from astropy.modeling.models import RotationSequence3D, Scale
+from astropy.modeling.models import Identity, RotationSequence3D, Scale, Shift
 from gwcs.geometry import CartesianToSpherical, SphericalToCartesian
 
+log = logging.getLogger(__name__)
+log.addHandler(logging.NullHandler())
+
 
 def v23tosky(input_model, wrap_v2_at=180, wrap_lon_at=360):
     """Create the transform from telescope to sky.
@@ -45,3 +50,66 @@ def v23tosky(input_model, wrap_v2_at=180, wrap_lon_at=360):
     )
     model.name = "v23tosky"
     return model
+
+
+def dva_corr_model(va_scale, v2_ref, v3_ref):
+    """
+    Create transformation that accounts for differential velocity aberration
+    (scale).
+
+    Parameters
+    ----------
+    va_scale : float, None
+        Ratio of the apparent plate scale to the true plate scale. When
+        ``va_scale`` is `None`, it is assumed to be identical to ``1`` and
+        an ``astropy.modeling.models.Identity`` model will be returned.
+
+    v2_ref : float, None
+        Telescope ``v2`` coordinate of the reference point in ``arcsec``. When
+        ``v2_ref`` is `None`, it is assumed to be identical to ``0``.
+
+    v3_ref : float, None
+        Telescope ``v3`` coordinate of the reference point in ``arcsec``. When
+        ``v3_ref`` is `None`, it is assumed to be identical to ``0``.
+
+    Returns
+    -------
+    va_corr : astropy.modeling.CompoundModel, astropy.modeling.models.Identity
+        A 2D compound model that corrects DVA. If ``va_scale`` is `None` or 1
+        then `astropy.modeling.models.Identity` will be returned.
+
+    """
+    if va_scale is None or va_scale == 1:
+        return Identity(2)
+
+    if va_scale <= 0:
+        log.warning("Given velocity aberration scale %s", va_scale)
+        log.warning(
+            "Velocity aberration scale must be a positive number. Setting to 1.0"
+        )
+        va_scale = 1.0
+
+    va_corr = Scale(va_scale, name="dva_scale_v2") & Scale(
+        va_scale, name="dva_scale_v3"
+    )
+
+    if v2_ref is None:
+        v2_ref = 0
+
+    if v3_ref is None:
+        v3_ref = 0
+
+    if v2_ref == 0 and v3_ref == 0:
+        return va_corr
+
+    # NOTE: it is assumed that v2, v3 angles and va scale are small enough
+    # so that for expected scale factors the issue of angle wrapping
+    # (180 degrees) can be neglected.
+    v2_shift = (1 - va_scale) * v2_ref
+    v3_shift = (1 - va_scale) * v3_ref
+
+    va_corr |= Shift(v2_shift, name="dva_v2_shift") & Shift(
+        v3_shift, name="dva_v3_shift"
+    )
+    va_corr.name = "DVA_Correction"
+    return va_corr