Add optika.radiometry with an effective-area model

optika/__init__.py

@@ -19,6 +19,7 @@
 from . import sensors
 from . import distortion
 from . import vignetting
+from . import radiometry
 from . import systems
 
 __all__ = [
@@ -41,5 +42,6 @@
     "sensors",
     "distortion",
     "vignetting",
+    "radiometry",
     "systems",
 ]

optika/radiometry/__init__.py

@@ -0,0 +1,11 @@
+"""Model the radiometry of an optical system."""
+
+from ._effective_area import (
+    AbstractEffectiveAreaModel,
+    InterpolatedEffectiveAreaModel,
+)
+
+__all__ = [
+    "AbstractEffectiveAreaModel",
+    "InterpolatedEffectiveAreaModel",
+]

optika/radiometry/_effective_area.py

@@ -0,0 +1,102 @@
+import abc
+import dataclasses
+import named_arrays as na
+import optika
+
+__all__ = [
+    "AbstractEffectiveAreaModel",
+    "InterpolatedEffectiveAreaModel",
+]
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class AbstractEffectiveAreaModel(
+    optika.mixins.Printable,
+):
+    """
+    An interface describing the effective area of an optical system as a
+    function of wavelength.
+    """
+
+    @abc.abstractmethod
+    def __call__(
+        self,
+        wavelength: na.AbstractScalar,
+    ) -> na.AbstractScalar:
+        """
+        The effective area of the optical system at the given wavelength.
+
+        Parameters
+        ----------
+        wavelength
+            The wavelength at which to evaluate the effective area.
+        """
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class InterpolatedEffectiveAreaModel(
+    AbstractEffectiveAreaModel,
+):
+    """
+    An effective area model which linearly interpolates between measured
+    calibration points.
+
+    Linear interpolation is used (rather than a polynomial fit) because the
+    effective area of a real system has sharp features --- absorption edges,
+    filter cutoffs, and multilayer reflectivity peaks --- that a polynomial
+    would fit poorly.
+
+    Examples
+    --------
+
+    Interpolate a measured effective area curve and plot the result.
+
+    .. jupyter-execute::
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        import astropy.units as u
+        import named_arrays as na
+        import optika
+
+        # A coarse set of calibration measurements
+        wavelength = na.linspace(100, 1000, axis="wavelength", num=10) * u.AA
+        area = 10 * np.exp(-(((wavelength - 500 * u.AA) / (150 * u.AA)) ** 2))
+        area = area * u.cm**2
+
+        model = optika.radiometry.InterpolatedEffectiveAreaModel(
+            wavelength=wavelength,
+            area=area,
+            axis_wavelength="wavelength",
+        )
+
+        # Evaluate the model on a finer grid
+        wavelength_fit = na.linspace(100, 1000, axis="wavelength", num=201) * u.AA
+
+        fig, ax = plt.subplots(constrained_layout=True)
+        na.plt.scatter(wavelength, area, ax=ax, label="calibration")
+        na.plt.plot(wavelength_fit, model(wavelength_fit), ax=ax, label="interpolated")
+        ax.set_xlabel(f"wavelength ({na.unit(wavelength):latex_inline})")
+        ax.set_ylabel(f"effective area ({na.unit(area):latex_inline})")
+        ax.legend();
+    """
+
+    wavelength: na.AbstractScalar = dataclasses.MISSING
+    """The wavelength of each calibration point."""
+
+    area: na.AbstractScalar = dataclasses.MISSING
+    """The measured effective area at each calibration point."""
+
+    axis_wavelength: str = dataclasses.MISSING
+    """The logical axis corresponding to changing wavelength."""
+
+    def __call__(
+        self,
+        wavelength: na.AbstractScalar,
+    ) -> na.AbstractScalar:
+        return na.interp(
+            wavelength,
+            self.wavelength,
+            self.area,
+            axis=self.axis_wavelength,
+        )

optika/radiometry/_effective_area_test.py

@@ -0,0 +1,58 @@
+import pytest
+import numpy as np
+import astropy.units as u
+import named_arrays as na
+import optika
+from .._tests import test_mixins
+
+
+def _wavelength() -> na.AbstractScalar:
+    return na.linspace(100, 1000, axis="wavelength", num=10) * u.AA
+
+
+def _area() -> na.AbstractScalar:
+    return na.linspace(1, 5, axis="wavelength", num=10) * u.cm**2
+
+
+class AbstractTestAbstractEffectiveAreaModel(
+    test_mixins.AbstractTestPrintable,
+):
+    def test__call__(self, a: optika.radiometry.AbstractEffectiveAreaModel):
+        wavelength = na.linspace(200, 900, axis="wavelength", num=5) * u.AA
+        result = a(wavelength)
+        assert isinstance(result, na.AbstractScalar)
+        assert na.unit_normalized(result).is_equivalent(u.cm**2)
+
+
+@pytest.mark.parametrize(
+    argnames="a",
+    argvalues=[
+        optika.radiometry.InterpolatedEffectiveAreaModel(
+            wavelength=_wavelength(),
+            area=_area(),
+            axis_wavelength="wavelength",
+        ),
+    ],
+)
+class TestInterpolatedEffectiveAreaModel(
+    AbstractTestAbstractEffectiveAreaModel,
+):
+    def test_wavelength(self, a: optika.radiometry.InterpolatedEffectiveAreaModel):
+        assert isinstance(a.wavelength, na.AbstractScalar)
+
+    def test_area(self, a: optika.radiometry.InterpolatedEffectiveAreaModel):
+        assert isinstance(a.area, na.AbstractScalar)
+
+    def test_axis_wavelength(
+        self,
+        a: optika.radiometry.InterpolatedEffectiveAreaModel,
+    ):
+        assert isinstance(a.axis_wavelength, str)
+
+    def test_interpolates_calibration(
+        self,
+        a: optika.radiometry.InterpolatedEffectiveAreaModel,
+    ):
+        # evaluating at the calibration wavelengths returns the calibration areas
+        result = a(a.wavelength)
+        assert np.all(result == a.area)