Adding functionality to represent non linear spectral dispersion

.github/workflows/build_workshop.yml

@@ -31,13 +31,13 @@ jobs:
         if: ${{ !startsWith(github.ref_name, 'v') }}
         run: echo 'BRANCHNAME=latest' >> "$GITHUB_ENV"
 
-      - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
+      - uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0
         with:
           fetch-depth: 0
           persist-credentials: false
 
       - name: Setup Python 3.13
-        uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0
+        uses: actions/setup-python@a309ff8b426b58ec0e2a45f0f869d46889d02405 # v6.2.0
         with:
           python-version: "3.13"
 
@@ -50,7 +50,7 @@ jobs:
         run: pip install -U pipx
 
       - name: Clone workshop repository
-        uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
+        uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0
         with:
           persist-credentials: true
           repository: 'DKISTDC/DKIST-Workshop'

changelog/713.feature.rst

@@ -0,0 +1 @@
+Added `dkist.wcs.models.build_grating_spectral_transform()` to build FITS ``-GRA``/``-GRI`` spectral transforms by composing a constant incident-angle term, `~gwcs.spectroscopy.RefractedAngleSineModel`, and `~gwcs.spectroscopy.WavelengthFromGrismEquation` using the FITS-facing grating parameters.

dkist/wcs/models.py

@@ -14,6 +14,7 @@
 import astropy.units as u
 from astropy.modeling import CompoundModel, Model, Parameter, separable
 from astropy.utils.decorators import deprecated_renamed_argument
+from gwcs.spectroscopy import RefractedAngleSineModel, WavelengthFromGrismEquation
 
 from dkist.utils.decorators import deprecated
 from dkist.utils.exceptions import DKISTDeprecationWarning
@@ -30,11 +31,63 @@
     "VaryingCelestialTransform",
     "VaryingCelestialTransform2D",
     "VaryingCelestialTransform3D",
+    "build_grating_spectral_transform",
     "generate_celestial_transform",
     "varying_celestial_transform_from_tables",
 ]
 
 
+def build_grating_spectral_transform(
+    reference_pixel: float,
+    reference_wavelength: u.Quantity,
+    dispersion: u.Quantity,
+    groove_density: u.Quantity,
+    spectral_order: u.Quantity,
+    incident_angle: u.Quantity,
+    refractive_index: u.Quantity = 1 * u.one,
+    refractive_index_derivative: u.Quantity = 0 / u.m,
+    out_of_plane_angle: u.Quantity = 0 * u.deg,
+    camera_angle: u.Quantity = 0 * u.deg,
+) -> CompoundModel:
+    """
+    Build a FITS grating spectral transform from header-derived parameters.
+
+    Composes a constant incident-angle sine term (`~astropy.modeling.models.Const1D`),
+    a pixel-dependent refracted-angle sine model
+    (`~gwcs.spectroscopy.RefractedAngleSineModel`), and
+    `~gwcs.spectroscopy.WavelengthFromGrismEquation`, following the FITS
+    grating/grism spectral-coordinate formalism described by Greisen et al.
+    (2006). The input and output angles are computed from the Greisen
+    relations within the component models:
+    https://scixplorer.org/abs/2006A%26A...446..747G/abstract
+    """
+    model = WavelengthFromGrismEquation(
+        groove_density=groove_density,
+        spectral_order=spectral_order,
+        reference_wavelength=reference_wavelength,
+        refractive_index=refractive_index,
+        refractive_index_derivative=refractive_index_derivative,
+        out_of_plane_angle=out_of_plane_angle,
+    )
+
+    alpha_in = m.Const1D(amplitude=np.sin(incident_angle))
+
+    alpha_out = RefractedAngleSineModel(
+        reference_pixel=reference_pixel,
+        reference_wavelength=reference_wavelength,
+        dispersion=dispersion,
+        groove_density=groove_density,
+        spectral_order=spectral_order,
+        incident_angle=incident_angle,
+        refractive_index=refractive_index,
+        refractive_index_derivative=refractive_index_derivative,
+        out_of_plane_angle=out_of_plane_angle,
+        camera_angle=camera_angle,
+    )
+
+    return m.Mapping((0, 0)) | (alpha_in & alpha_out) | model
+
+
 def generate_celestial_transform(
         crpix: Iterable[float] | u.Quantity,
         cdelt: Iterable[float] | u.Quantity,

dkist/wcs/tests/test_models.py

@@ -9,10 +9,13 @@
 from astropy.coordinates.matrix_utilities import rotation_matrix
 from astropy.modeling import CompoundModel
 from astropy.modeling.models import Tabular1D
+from astropy.wcs import WCS
+from gwcs.spectroscopy import RefractedAngleSineModel
 
 from dkist.wcs.models import (AsymmetricMapping, Ravel, Unravel, VaryingCelestialTransform,
                               VaryingCelestialTransform2D, VaryingCelestialTransform3D,
-                              generate_celestial_transform, update_celestial_transform_parameters,
+                              build_grating_spectral_transform, generate_celestial_transform,
+                              update_celestial_transform_parameters,
                               varying_celestial_transform_from_tables)
 
 
@@ -52,6 +55,49 @@ def test_generate_celestial_unitless():
     assert u.allclose(shift1.offset, 0)
 
 
+def test_build_grating_spectral_transform() -> None:
+    header = {
+        "CTYPE1": "AWAV-GRA",
+        "CUNIT1": "nm",
+        "CRPIX1": 218,
+        "CRVAL1": 854.1738582455826,
+        "CDELT1": 0.0022975580183395555,
+        "PV1_0": 23000.0,
+        "PV1_1": 90,
+        "PV1_2": 65.696,
+        "PV1_3": 1.25,
+        "PV1_4": 1000.0,
+        "PV1_5": 1.5,
+        "PV1_6": 0.8,
+    }
+    transform = build_grating_spectral_transform(
+        reference_pixel=header["CRPIX1"] - 1,
+        reference_wavelength=header["CRVAL1"] * u.nm,
+        dispersion=header["CDELT1"] * u.nm / u.pix,
+        groove_density=header["PV1_0"] / u.m,
+        spectral_order=header["PV1_1"] * u.one,
+        incident_angle=header["PV1_2"] * u.deg,
+        refractive_index=header["PV1_3"] * u.one,
+        refractive_index_derivative=header["PV1_4"] / u.m,
+        out_of_plane_angle=header["PV1_5"] * u.deg,
+        camera_angle=header["PV1_6"] * u.deg,
+    )
+
+    pixels = np.array([0, 100, 217, 300, 511], dtype=float)
+    expected = WCS(header).spectral.pixel_to_world(pixels)
+    result = transform(pixels)
+
+    assert isinstance(transform, CompoundModel)
+    # The transform should contain the gwcs refracted-angle model.
+    assert any(
+        isinstance(sm, RefractedAngleSineModel)
+        for sm in transform.traverse_postorder()
+    )
+    np.testing.assert_allclose(
+        result.to_value(u.nm), expected.to_value(u.nm), rtol=1e-10, atol=1e-10
+    )
+
+
 def test_update_celestial():
     trsfm  = generate_celestial_transform(
               crpix=[0, 0] * u.pix,