From 413c8c147e4db2484a72d88dd02f9345fc164a90 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 09:09:47 -0600
Subject: [PATCH 1/8] Add `optika.distortion`, a polynomial model of optical
 distortion.

Introduce the `optika.distortion` subpackage, which maps scene
coordinates (wavelength + field position) to sensor coordinates and
back:

- `AbstractDistortionModel` defines the `distort`/`undistort` interface.
- `AbstractInterpolatedDistortionModel` adds the calibration-point
  interface (`coordinates_scene`, `coordinates_sensor`, axes).
- `PolynomialDistortionModel` fits forward and inverse polynomials (with
  a mean-centered reference for conditioning) and exposes them as the
  public `fit`/`fit_inverse` cached properties.

Register the subpackage in `optika/__init__.py` and add round-trip tests.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/__init__.py                    |   2 +
 optika/distortion/__init__.py         |  13 ++
 optika/distortion/_distortion.py      | 227 ++++++++++++++++++++++++++
 optika/distortion/_distortion_test.py |  99 +++++++++++
 4 files changed, 341 insertions(+)
 create mode 100644 optika/distortion/__init__.py
 create mode 100644 optika/distortion/_distortion.py
 create mode 100644 optika/distortion/_distortion_test.py

diff --git a/optika/__init__.py b/optika/__init__.py
index 1f4d05a..bbb8ea4 100644
--- a/optika/__init__.py
+++ b/optika/__init__.py
@@ -17,6 +17,7 @@
 from . import rulings
 from . import surfaces
 from . import sensors
+from . import distortion
 from . import systems
 
 __all__ = [
@@ -37,5 +38,6 @@
     "rulings",
     "surfaces",
     "sensors",
+    "distortion",
     "systems",
 ]
diff --git a/optika/distortion/__init__.py b/optika/distortion/__init__.py
new file mode 100644
index 0000000..01f41c1
--- /dev/null
+++ b/optika/distortion/__init__.py
@@ -0,0 +1,13 @@
+"""Model the distortion of a scene observed by an optical system."""
+
+from ._distortion import (
+    AbstractDistortionModel,
+    AbstractInterpolatedDistortionModel,
+    PolynomialDistortionModel,
+)
+
+__all__ = [
+    "AbstractDistortionModel",
+    "AbstractInterpolatedDistortionModel",
+    "PolynomialDistortionModel",
+]
diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
new file mode 100644
index 0000000..d45f076
--- /dev/null
+++ b/optika/distortion/_distortion.py
@@ -0,0 +1,227 @@
+import abc
+import dataclasses
+import functools
+import named_arrays as na
+import optika
+
+__all__ = [
+    "AbstractDistortionModel",
+    "AbstractInterpolatedDistortionModel",
+    "PolynomialDistortionModel",
+]
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class AbstractDistortionModel(
+    optika.mixins.Printable,
+):
+    """
+    An interface describing an arbitrary distortion model,
+    which maps scene coordinates to sensor coordinates (and vice versa).
+
+    A distortion model carries the wavelength along with the position, since
+    the mapping from a point in the scene to a point on the sensor generally
+    depends on wavelength (for example, the dispersion of a spectrograph).
+    As a result :meth:`distort` and :meth:`undistort` are inverses of one
+    another only up to the accuracy of the model.
+    """
+
+    @abc.abstractmethod
+    def distort(
+        self,
+        coordinates: optika.vectors.SceneVectorArray,
+    ) -> na.SpectralPositionalVectorArray:
+        """
+        Convert scene coordinates to sensor coordinates.
+
+        Parameters
+        ----------
+        coordinates
+            The wavelength and field position of each point in the scene.
+        """
+
+    @abc.abstractmethod
+    def undistort(
+        self,
+        coordinates: na.SpectralPositionalVectorArray,
+    ) -> optika.vectors.SceneVectorArray:
+        """
+        Convert sensor coordinates to scene coordinates.
+
+        Parameters
+        ----------
+        coordinates
+            The wavelength and sensor position of each point.
+        """
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class AbstractInterpolatedDistortionModel(
+    AbstractDistortionModel,
+):
+    """
+    A distortion model defined by interpolating between known scene/sensor
+    coordinates.
+
+    This class has two main members, :attr:`coordinates_scene` and
+    :attr:`coordinates_sensor`, the calibration points between which subclasses
+    interpolate.
+    """
+
+    @property
+    @abc.abstractmethod
+    def coordinates_scene(self) -> optika.vectors.SceneVectorArray:
+        """
+        The wavelength and field position of each calibration point in the scene.
+        """
+
+    @property
+    @abc.abstractmethod
+    def coordinates_sensor(self) -> na.Cartesian2dVectorArray:
+        """
+        The position of each calibration point mapped onto the sensor.
+        """
+
+    @property
+    @abc.abstractmethod
+    def axis_wavelength(self) -> str:
+        """The logical axis corresponding to changing wavelength."""
+
+    @property
+    @abc.abstractmethod
+    def axis_field(self) -> tuple[str, str]:
+        """The logical axes corresponding to changing position in the scene."""
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class PolynomialDistortionModel(
+    AbstractInterpolatedDistortionModel,
+):
+    """
+    A distortion model which fits a polynomial to known scene/sensor coordinates.
+
+    The forward model (:meth:`distort`) is a polynomial fit mapping scene field
+    position to sensor position as a function of wavelength.
+    The inverse model (:meth:`undistort`) is a *separate* polynomial fit in the
+    opposite direction, so the round trip is exact only to the accuracy of the
+    two fits.
+
+    Examples
+    --------
+
+    Build a distortion model from a grid of scene/sensor calibration points and
+    confirm that the round trip recovers the original scene coordinates.
+
+    .. jupyter-execute::
+
+        import numpy as np
+        import astropy.units as u
+        import named_arrays as na
+        import optika
+
+        # Define a grid of known field positions and wavelengths
+        scene = optika.vectors.SceneVectorArray(
+            wavelength=na.linspace(500, 600, axis="wavelength", num=3) * u.nm,
+            field=na.Cartesian2dVectorLinearSpace(
+                start=-1 * u.deg,
+                stop=+1 * u.deg,
+                axis=na.Cartesian2dVectorArray("field_x", "field_y"),
+                num=5,
+            ),
+        )
+
+        # Map them onto the sensor with a simple linear distortion
+        sensor = na.Cartesian2dVectorArray(
+            x=scene.field.x * (10 * u.mm / u.deg),
+            y=scene.field.y * (10 * u.mm / u.deg),
+        )
+
+        model = optika.distortion.PolynomialDistortionModel(
+            coordinates_scene=scene,
+            coordinates_sensor=sensor,
+            axis_wavelength="wavelength",
+            axis_field=("field_x", "field_y"),
+            degree=1,
+        )
+
+        scene_roundtrip = model.undistort(model.distort(scene))
+    """
+
+    coordinates_scene: optika.vectors.SceneVectorArray = dataclasses.MISSING
+    """The wavelength and field position of each calibration point in the scene."""
+
+    coordinates_sensor: na.Cartesian2dVectorArray = dataclasses.MISSING
+    """The position of each calibration point mapped onto the sensor."""
+
+    axis_wavelength: str = dataclasses.MISSING
+    """The logical axis corresponding to changing wavelength."""
+
+    axis_field: tuple[str, str] = dataclasses.MISSING
+    """The logical axes corresponding to changing position in the scene."""
+
+    degree: int = 1
+    """The degree of the polynomial used to model the distortion."""
+
+    where: bool | na.AbstractScalar = True
+    """A boolean mask selecting which calibration points to use for fitting."""
+
+    @property
+    def _axis_scene(self) -> tuple[str, ...]:
+        """The logical axes over which the calibration points are distributed."""
+        return (self.axis_wavelength, *self.axis_field)
+
+    @functools.cached_property
+    def fit(self) -> na.PolynomialFitFunctionArray:
+        """The polynomial fit mapping scene field position to sensor position."""
+        scene = self.coordinates_scene
+        inputs = na.SpectralPositionalVectorArray(
+            wavelength=scene.wavelength,
+            position=scene.field,
+        )
+        return na.PolynomialFitFunctionArray(
+            inputs=inputs,
+            outputs=self.coordinates_sensor,
+            center=inputs.mean(self._axis_scene),
+            degree=self.degree,
+            where_polynomial=self.where,
+        )
+
+    @functools.cached_property
+    def fit_inverse(self) -> na.PolynomialFitFunctionArray:
+        """The polynomial fit mapping sensor position back to scene field position."""
+        scene = self.coordinates_scene
+        inputs = na.SpectralPositionalVectorArray(
+            wavelength=scene.wavelength,
+            position=self.coordinates_sensor,
+        )
+        return na.PolynomialFitFunctionArray(
+            inputs=inputs,
+            outputs=scene.field,
+            center=inputs.mean(self._axis_scene),
+            degree=self.degree,
+            where_polynomial=self.where,
+        )
+
+    def distort(
+        self,
+        coordinates: optika.vectors.SceneVectorArray,
+    ) -> na.SpectralPositionalVectorArray:
+        inputs = na.SpectralPositionalVectorArray(
+            wavelength=coordinates.wavelength,
+            position=coordinates.field,
+        )
+        position = self.fit(inputs).outputs
+        return na.SpectralPositionalVectorArray(
+            wavelength=coordinates.wavelength,
+            position=position,
+        )
+
+    def undistort(
+        self,
+        coordinates: na.SpectralPositionalVectorArray,
+    ) -> optika.vectors.SceneVectorArray:
+        field = self.fit_inverse(coordinates).outputs
+        return optika.vectors.SceneVectorArray(
+            wavelength=coordinates.wavelength,
+            field=field,
+        )
diff --git a/optika/distortion/_distortion_test.py b/optika/distortion/_distortion_test.py
new file mode 100644
index 0000000..c6debaf
--- /dev/null
+++ b/optika/distortion/_distortion_test.py
@@ -0,0 +1,99 @@
+import pytest
+import numpy as np
+import astropy.units as u
+import named_arrays as na
+import optika
+from .._tests import test_mixins
+
+
+def _scene() -> optika.vectors.SceneVectorArray:
+    return optika.vectors.SceneVectorArray(
+        wavelength=na.linspace(500, 600, axis="wavelength", num=3) * u.nm,
+        field=na.Cartesian2dVectorLinearSpace(
+            start=-1 * u.deg,
+            stop=+1 * u.deg,
+            axis=na.Cartesian2dVectorArray("field_x", "field_y"),
+            num=5,
+        ),
+    )
+
+
+class AbstractTestAbstractDistortionModel(
+    test_mixins.AbstractTestPrintable,
+):
+    def test_distort(self, a: optika.distortion.AbstractDistortionModel):
+        coordinates = _scene()
+        result = a.distort(coordinates)
+        assert isinstance(result, na.SpectralPositionalVectorArray)
+        assert isinstance(result.position, na.AbstractCartesian2dVectorArray)
+        # the wavelength is carried through unchanged
+        assert np.all(result.wavelength == coordinates.wavelength)
+
+    def test_undistort(self, a: optika.distortion.AbstractDistortionModel):
+        coordinates = a.distort(_scene())
+        result = a.undistort(coordinates)
+        assert isinstance(result, optika.vectors.SceneVectorArray)
+        assert np.all(result.wavelength == coordinates.wavelength)
+
+    def test_roundtrip(self, a: optika.distortion.AbstractDistortionModel):
+        scene = _scene()
+        result = a.undistort(a.distort(scene))
+        error = (result.field - scene.field).length
+        assert np.all(error < 1e-9 * u.deg)
+
+
+class AbstractTestAbstractInterpolatedDistortionModel(
+    AbstractTestAbstractDistortionModel,
+):
+    def test_coordinates_scene(
+        self,
+        a: optika.distortion.AbstractInterpolatedDistortionModel,
+    ):
+        assert isinstance(a.coordinates_scene, optika.vectors.SceneVectorArray)
+
+    def test_coordinates_sensor(
+        self,
+        a: optika.distortion.AbstractInterpolatedDistortionModel,
+    ):
+        assert isinstance(a.coordinates_sensor, na.AbstractCartesian2dVectorArray)
+
+    def test_axis_wavelength(
+        self,
+        a: optika.distortion.AbstractInterpolatedDistortionModel,
+    ):
+        assert isinstance(a.axis_wavelength, str)
+
+    def test_axis_field(
+        self,
+        a: optika.distortion.AbstractInterpolatedDistortionModel,
+    ):
+        assert isinstance(a.axis_field, tuple)
+        assert all(isinstance(ax, str) for ax in a.axis_field)
+
+
+@pytest.mark.parametrize(
+    argnames="a",
+    argvalues=[
+        optika.distortion.PolynomialDistortionModel(
+            coordinates_scene=_scene(),
+            coordinates_sensor=na.Cartesian2dVectorArray(
+                x=_scene().field.x * (10 * u.mm / u.deg),
+                y=_scene().field.y * (10 * u.mm / u.deg),
+            ),
+            axis_wavelength="wavelength",
+            axis_field=("field_x", "field_y"),
+            degree=degree,
+        )
+        for degree in [1, 2]
+    ],
+)
+class TestPolynomialDistortionModel(
+    AbstractTestAbstractInterpolatedDistortionModel,
+):
+    def test_fit(self, a: optika.distortion.PolynomialDistortionModel):
+        assert isinstance(a.fit, na.PolynomialFitFunctionArray)
+        assert a.fit.degree == a.degree
+
+    def test_fit_inverse(self, a: optika.distortion.PolynomialDistortionModel):
+        assert isinstance(a.fit_inverse, na.PolynomialFitFunctionArray)
+        assert a.fit_inverse.degree == a.degree

From 01bb9c8baeba3e150a5c7fbd84b4378604b1b0c7 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 12:11:14 -0600
Subject: [PATCH 2/8] Require named-arrays~=1.5 for
 PolynomialFitFunctionArray.center.

The distortion model relies on named-arrays features (the polynomial fit
`center` field) that are only available in named-arrays 1.5+.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 850dafb..165a955 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,7 +15,7 @@ classifiers = [
 ]
 dependencies = [
     "astropy!=6.1.5",
-    "named-arrays~=1.0",
+    "named-arrays~=1.5",
     "pymupdf",
     "joblib",
     "ezdxf",

From c7560369badc983e46016b3441d9eb0807e54575 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 13:14:20 -0600
Subject: [PATCH 3/8] Add `PolynomialDistortionModel.plot_residual()`.

Plot the magnitude of the forward-fit residual (calibration sensor
positions minus fit predictions) as a function of field angle using
`na.plt.pcolormesh`, with a separate subplot for each wavelength and a
shared colorbar. Update the class example to demonstrate it.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/distortion/_distortion.py      | 92 ++++++++++++++++++++++++---
 optika/distortion/_distortion_test.py |  9 +++
 2 files changed, 93 insertions(+), 8 deletions(-)

diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
index d45f076..d447f86 100644
--- a/optika/distortion/_distortion.py
+++ b/optika/distortion/_distortion.py
@@ -1,6 +1,12 @@
 import abc
 import dataclasses
 import functools
+import matplotlib.axes
+import matplotlib.cm
+import matplotlib.colors
+import matplotlib.figure
+import matplotlib.pyplot as plt
+import astropy.visualization
 import named_arrays as na
 import optika
 
@@ -109,8 +115,8 @@ class PolynomialDistortionModel(
     Examples
     --------
 
-    Build a distortion model from a grid of scene/sensor calibration points and
-    confirm that the round trip recovers the original scene coordinates.
+    Plot the fit residual of a distortion model with a deliberately
+    underfit (linear) polynomial.
 
     .. jupyter-execute::
 
@@ -119,20 +125,18 @@ class PolynomialDistortionModel(
         import named_arrays as na
         import optika
 
-        # Define a grid of known field positions and wavelengths
         scene = optika.vectors.SceneVectorArray(
             wavelength=na.linspace(500, 600, axis="wavelength", num=3) * u.nm,
             field=na.Cartesian2dVectorLinearSpace(
                 start=-1 * u.deg,
                 stop=+1 * u.deg,
                 axis=na.Cartesian2dVectorArray("field_x", "field_y"),
-                num=5,
+                num=11,
             ),
         )
-
-        # Map them onto the sensor with a simple linear distortion
         sensor = na.Cartesian2dVectorArray(
-            x=scene.field.x * (10 * u.mm / u.deg),
+            x=scene.field.x * (10 * u.mm / u.deg)
+            + scene.field.x**2 * (1 * u.mm / u.deg**2),
             y=scene.field.y * (10 * u.mm / u.deg),
         )
 
@@ -144,7 +148,7 @@ class PolynomialDistortionModel(
             degree=1,
         )
 
-        scene_roundtrip = model.undistort(model.distort(scene))
+        fig, ax = model.plot_residual()
     """
 
     coordinates_scene: optika.vectors.SceneVectorArray = dataclasses.MISSING
@@ -225,3 +229,75 @@ def undistort(
             wavelength=coordinates.wavelength,
             field=field,
         )
+
+    def plot_residual(
+        self,
+        cmap: None | str | matplotlib.colors.Colormap = None,
+        **kwargs,
+    ) -> tuple[matplotlib.figure.Figure, na.ScalarArray]:
+        """
+        Plot the residual of the forward :attr:`fit` as a function of field
+        angle, with a separate subplot for each wavelength.
+
+        The residual is the magnitude of the difference between the calibration
+        sensor positions, :attr:`coordinates_sensor`, and the positions
+        predicted by the forward polynomial fit.
+
+        Parameters
+        ----------
+        cmap
+            The colormap used to map the residual magnitude to colors.
+        kwargs
+            Additional keyword arguments passed to
+            :func:`named_arrays.plt.pcolormesh`.
+        """
+        scene = self.coordinates_scene
+        field = scene.field
+        wavelength = na.as_named_array(scene.wavelength)
+        axis_wavelength = self.axis_wavelength
+
+        residual = (self.coordinates_sensor - self.fit.predictions).length
+        unit = na.unit(residual)
+
+        if axis_wavelength in na.shape(wavelength):
+            ncols = wavelength.shape[axis_wavelength]
+        else:
+            ncols = 1
+
+        with astropy.visualization.quantity_support():
+            fig, ax = na.plt.subplots(
+                axis_cols=axis_wavelength,
+                ncols=ncols,
+                sharex=True,
+                sharey=True,
+                squeeze=False,
+                constrained_layout=True,
+            )
+
+            colorizer = plt.Colorizer(
+                cmap=cmap,
+                norm=plt.Normalize(vmin=0, vmax=residual.ndarray.max().value),
+            )
+
+            na.plt.pcolormesh(
+                field,
+                C=residual,
+                ax=ax,
+                colorizer=colorizer,
+                **kwargs,
+            )
+
+            na.plt.set_xlabel(f"field $x$ ({na.unit(field.x):latex_inline})", ax=ax)
+            na.plt.set_ylabel(
+                f"field $y$ ({na.unit(field.y):latex_inline})",
+                ax=ax[{axis_wavelength: 0}],
+            )
+            na.plt.set_title(wavelength.to_string_array(), ax=ax)
+
+            plt.colorbar(
+                mappable=matplotlib.cm.ScalarMappable(colorizer=colorizer),
+                ax=ax.ndarray,
+                label=f"residual ({unit:latex_inline})",
+            )
+
+        return fig, ax
diff --git a/optika/distortion/_distortion_test.py b/optika/distortion/_distortion_test.py
index c6debaf..be6d109 100644
--- a/optika/distortion/_distortion_test.py
+++ b/optika/distortion/_distortion_test.py
@@ -1,6 +1,8 @@
 import pytest
 import numpy as np
 import astropy.units as u
+import matplotlib.figure
+import matplotlib.pyplot as plt
 import named_arrays as na
 import optika
 from .._tests import test_mixins
@@ -97,3 +99,10 @@ def test_fit(self, a: optika.distortion.PolynomialDistortionModel):
     def test_fit_inverse(self, a: optika.distortion.PolynomialDistortionModel):
         assert isinstance(a.fit_inverse, na.PolynomialFitFunctionArray)
         assert a.fit_inverse.degree == a.degree
+
+    def test_plot_residual(self, a: optika.distortion.PolynomialDistortionModel):
+        fig, ax = a.plot_residual()
+        assert isinstance(fig, matplotlib.figure.Figure)
+        assert isinstance(ax, na.ScalarArray)
+        assert a.axis_wavelength in na.shape(ax)
+        plt.close(fig)

From a094ffa62d08b956a5544337d352d1851793c9c7 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 13:47:09 -0600
Subject: [PATCH 4/8] Simplify plot_residual column count for full test
 coverage.

Replace the if/else branch selecting the subplot count with
`na.shape(wavelength).get(axis_wavelength, 1)`, which the existing tests
fully exercise (the scalar-wavelength else-branch was never hit).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/distortion/_distortion.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
index d447f86..c305065 100644
--- a/optika/distortion/_distortion.py
+++ b/optika/distortion/_distortion.py
@@ -259,10 +259,7 @@ def plot_residual(
         residual = (self.coordinates_sensor - self.fit.predictions).length
         unit = na.unit(residual)
 
-        if axis_wavelength in na.shape(wavelength):
-            ncols = wavelength.shape[axis_wavelength]
-        else:
-            ncols = 1
+        ncols = na.shape(wavelength).get(axis_wavelength, 1)
 
         with astropy.visualization.quantity_support():
             fig, ax = na.plt.subplots(

From 389943c92f9ba304474846e2628912eae053ea96 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 13:47:09 -0600
Subject: [PATCH 5/8] Set MPLBACKEND=agg in the tests workflow.

Use a non-interactive matplotlib backend in CI so plotting tests (e.g.
`PolynomialDistortionModel.plot_residual`) run headless, matching the
named-arrays workflow.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 .github/workflows/tests.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 471e50e..9597ea4 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -32,6 +32,8 @@ jobs:
           pip install setuptools wheel
           pip install -e .[test]
       - name: Test with pytest
+        env:
+          MPLBACKEND: "agg"
         run: |
           pip install pytest pytest-cov
           pytest --doctest-modules --junitxml=junit/test-results.xml --cov=. --cov-report=xml --cov-report=html

From 5d393b7843d7c3a1a8c82b0ef08a6ea70e16171f Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 14:21:52 -0600
Subject: [PATCH 6/8] Improve plot_residual example: quadratic in both axes,
 equal aspect.

Make the example sensor distortion quadratic in both field components,
increase the field sampling, and set the subplot aspect ratio to equal.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/distortion/_distortion.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
index c305065..330e798 100644
--- a/optika/distortion/_distortion.py
+++ b/optika/distortion/_distortion.py
@@ -131,13 +131,14 @@ class PolynomialDistortionModel(
                 start=-1 * u.deg,
                 stop=+1 * u.deg,
                 axis=na.Cartesian2dVectorArray("field_x", "field_y"),
-                num=11,
+                num=13,
             ),
         )
         sensor = na.Cartesian2dVectorArray(
             x=scene.field.x * (10 * u.mm / u.deg)
             + scene.field.x**2 * (1 * u.mm / u.deg**2),
-            y=scene.field.y * (10 * u.mm / u.deg),
+            y=scene.field.y * (10 * u.mm / u.deg)
+            + scene.field.y**2 * (1 * u.mm / u.deg**2),
         )
 
         model = optika.distortion.PolynomialDistortionModel(
@@ -149,6 +150,7 @@ class PolynomialDistortionModel(
         )
 
         fig, ax = model.plot_residual()
+        na.plt.set_aspect("equal", ax=ax)
     """
 
     coordinates_scene: optika.vectors.SceneVectorArray = dataclasses.MISSING

From 142e9f5c86286aec302a9acc40711a0b1e1a5789 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 15:16:56 -0600
Subject: [PATCH 7/8] Add figsize, vmin, and vmax keyword arguments to
 plot_residual().

`vmin`/`vmax` set the residual color limits (defaulting to zero and the
maximum residual). `figsize` overrides the figure size; when omitted it
is now derived automatically from the number of wavelength subplots and
the field-of-view aspect ratio (via `field.x.ptp() / field.y.ptp()`).
Extend the test to exercise the new arguments.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/distortion/_distortion.py      | 34 ++++++++++++++++++++++++++-
 optika/distortion/_distortion_test.py | 15 ++++++++++--
 2 files changed, 46 insertions(+), 3 deletions(-)

diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
index 330e798..0d4aaf3 100644
--- a/optika/distortion/_distortion.py
+++ b/optika/distortion/_distortion.py
@@ -234,7 +234,10 @@ def undistort(
 
     def plot_residual(
         self,
+        figsize: None | tuple[float, float] = None,
         cmap: None | str | matplotlib.colors.Colormap = None,
+        vmin: None | na.ArrayLike = None,
+        vmax: None | na.ArrayLike = None,
         **kwargs,
     ) -> tuple[matplotlib.figure.Figure, na.ScalarArray]:
         """
@@ -247,8 +250,18 @@ def plot_residual(
 
         Parameters
         ----------
+        figsize
+            The size of the returned figure in inches.
+            If :obj:`None`, the size is chosen automatically from the number
+            of wavelengths and the aspect ratio of the field of view.
         cmap
             The colormap used to map the residual magnitude to colors.
+        vmin
+            The residual value mapped to the lowest color.
+            If :obj:`None`, defaults to zero.
+        vmax
+            The residual value mapped to the highest color.
+            If :obj:`None`, defaults to the maximum residual.
         kwargs
             Additional keyword arguments passed to
             :func:`named_arrays.plt.pcolormesh`.
@@ -261,8 +274,23 @@ def plot_residual(
         residual = (self.coordinates_sensor - self.fit.predictions).length
         unit = na.unit(residual)
 
+        if vmin is None:
+            vmin = 0 * unit
+        if vmax is None:
+            vmax = residual.max()
+
         ncols = na.shape(wavelength).get(axis_wavelength, 1)
 
+        if figsize is None:
+            # shape each subplot to the field-of-view aspect ratio, and widen
+            # the figure to fit one subplot per wavelength
+            height_subplot = 3
+            aspect = (field.x.ptp() / field.y.ptp()).ndarray.value
+            figsize = (
+                ncols * height_subplot * aspect + 1.5,
+                height_subplot + 1,
+            )
+
         with astropy.visualization.quantity_support():
             fig, ax = na.plt.subplots(
                 axis_cols=axis_wavelength,
@@ -270,12 +298,16 @@ def plot_residual(
                 sharex=True,
                 sharey=True,
                 squeeze=False,
+                figsize=figsize,
                 constrained_layout=True,
             )
 
             colorizer = plt.Colorizer(
                 cmap=cmap,
-                norm=plt.Normalize(vmin=0, vmax=residual.ndarray.max().value),
+                norm=plt.Normalize(
+                    vmin=na.as_named_array(vmin).ndarray.to_value(unit),
+                    vmax=na.as_named_array(vmax).ndarray.to_value(unit),
+                ),
             )
 
             na.plt.pcolormesh(
diff --git a/optika/distortion/_distortion_test.py b/optika/distortion/_distortion_test.py
index be6d109..832532b 100644
--- a/optika/distortion/_distortion_test.py
+++ b/optika/distortion/_distortion_test.py
@@ -100,8 +100,19 @@ def test_fit_inverse(self, a: optika.distortion.PolynomialDistortionModel):
         assert isinstance(a.fit_inverse, na.PolynomialFitFunctionArray)
         assert a.fit_inverse.degree == a.degree
 
-    def test_plot_residual(self, a: optika.distortion.PolynomialDistortionModel):
-        fig, ax = a.plot_residual()
+    @pytest.mark.parametrize(
+        argnames="kwargs",
+        argvalues=[
+            dict(),
+            dict(figsize=(8, 4), cmap="viridis", vmin=0 * u.um, vmax=5 * u.um),
+        ],
+    )
+    def test_plot_residual(
+        self,
+        a: optika.distortion.PolynomialDistortionModel,
+        kwargs: dict,
+    ):
+        fig, ax = a.plot_residual(**kwargs)
         assert isinstance(fig, matplotlib.figure.Figure)
         assert isinstance(ax, na.ScalarArray)
         assert a.axis_wavelength in na.shape(ax)

From ccc1d2aa0fc623e2d6b49f9a5fc55c92c5cd2522 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 30 May 2026 16:13:27 -0600
Subject: [PATCH 8/8] Suppress plot_residual example output with a trailing
 semicolon.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 optika/distortion/_distortion.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py
index 0d4aaf3..e29da3a 100644
--- a/optika/distortion/_distortion.py
+++ b/optika/distortion/_distortion.py
@@ -150,7 +150,7 @@ class PolynomialDistortionModel(
         )
 
         fig, ax = model.plot_residual()
-        na.plt.set_aspect("equal", ax=ax)
+        na.plt.set_aspect("equal", ax=ax);
     """
 
     coordinates_scene: optika.vectors.SceneVectorArray = dataclasses.MISSING