diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index a7d43b9..a7bfe48 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 --cov=. --cov-report=xml
diff --git a/ctis/__init__.py b/ctis/__init__.py
index fe5117f..2708fa6 100644
--- a/ctis/__init__.py
+++ b/ctis/__init__.py
@@ -5,8 +5,10 @@
 
 from . import scenes
 from . import instruments
+from . import inverters
 
 __all__ = [
     "scenes",
     "instruments",
+    "inverters",
 ]
diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py
index 748a579..f65ae17 100644
--- a/ctis/instruments/_instruments.py
+++ b/ctis/instruments/_instruments.py
@@ -114,12 +114,19 @@ def uncertainty(self) -> Callable[[na.ScalarArray], na.ScalarArray]:
         for a given number of photons.
         """
 
+    @property
+    @abc.abstractmethod
+    def channel(self):
+        """
+        Human-readable name of each independent CTIS channel.
+        """
+
     @property
     @abc.abstractmethod
     def axis_channel(self) -> str | tuple[str, ...]:
         """
         The logical axis or axes of this instrument corresponding to
-        the different dispersion magnitudes and angles.
+        the different CTIS channels.
         """
 
     @property
@@ -391,6 +398,11 @@ class IdealInstrument(
     A grid of wavelength and position coordinates on the sensor plane.
     """
 
+    channel: str | na.AbstractScalar = dataclasses.MISSING
+    """
+    Human-readable name of each independent CTIS channel.
+    """
+
     axis_channel: str | tuple[str, ...] = dataclasses.MISSING
     """
     The logical axis or axes of this instrument corresponding to
@@ -492,10 +504,6 @@ def _coordinates_output(self) -> na.AbstractSpectralPositionalVectorArray:
 
         coordinates_output = coordinates_output.cell_centers(self.axis_wavelength)
 
-        p = coordinates_output.position
-        coordinates_output.position.x = na.random.normal(p.x, 1e-3 * u.pix)
-        coordinates_output.position.y = na.random.normal(p.y, 1e-3 * u.pix)
-
         return coordinates_output
 
     @functools.cached_property
@@ -518,12 +526,12 @@ def weights_transpose(self):
         coordinates_input = self._coordinates_input
         coordinates_output = self._coordinates_output
 
-        return na.regridding.weights(
-            coordinates_input=coordinates_output.position,
-            coordinates_output=coordinates_input.position,
-            axis_input=self.axis_sensor_xy,
-            axis_output=self.axis_scene_xy,
-            method="conservative",
+        return na.regridding.transpose_weights_conservative(
+            weights=self.weights,
+            coordinates_input=coordinates_input.position,
+            coordinates_output=coordinates_output.position,
+            axis_input=self.axis_scene_xy,
+            axis_output=self.axis_sensor_xy,
         )
 
     def image(
diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py
index 0e06ff1..9f1ad33 100644
--- a/ctis/instruments/_instruments_test.py
+++ b/ctis/instruments/_instruments_test.py
@@ -7,6 +7,8 @@
 
 velocity = na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s
 
+wavelength_rest = 171 * u.AA
+
 position_scene = na.Cartesian2dVectorLinearSpace(
     start=-20 * u.arcsec,
     stop=20 * u.arcsec,
@@ -19,15 +21,18 @@
     y=na.arange(0, 64, axis="sensor_y") * u.pix,
 )
 
-coordinates_scene = na.SpectralPositionalVectorArray(velocity, position_scene)
-coordinates_sensor = na.SpectralPositionalVectorArray(velocity, position_sensor)
-
-gaussians = ctis.scenes.gaussians(
-    inputs=coordinates_scene,
-    width=na.SpectralPositionalVectorArray(30 * u.km / u.s, 1 * u.arcsec),
+coordinates_scene = na.DopplerPositionalVectorArray.from_velocity(
+    velocity=velocity,
+    wavelength_rest=wavelength_rest,
+    position=position_scene,
+)
+coordinates_sensor = na.DopplerPositionalVectorArray.from_velocity(
+    velocity=velocity,
+    wavelength_rest=wavelength_rest,
+    position=position_sensor,
 )
 
-wavelength_rest = 171 * u.AA
+gaussians = ctis.scenes.gaussians(coordinates_scene)
 
 AA = dict(
     unit=u.AA,
@@ -40,16 +45,19 @@
 dispersion = 200 * u.km / u.s
 dispersion = (dispersion.to(**AA) - wavelength_rest) / u.pix
 
+angle = na.linspace(0, 360, axis="channel", num=3, endpoint=False)
+
 instrument_ideal = ctis.instruments.IdealInstrument(
     area_effective=1 * u.cm**2,
     timedelta_exposure=10 * u.s,
     plate_scale=2 * u.arcsec / u.pix,
     dispersion=dispersion,
-    angle=na.linspace(0, 360, axis="channel", num=3, endpoint=False),
+    angle=angle,
     wavelength_ref=wavelength_rest,
     position_ref=32 * u.pix,
     coordinates_scene=coordinates_scene,
     coordinates_sensor=coordinates_sensor,
+    channel=angle,
     axis_channel="channel",
     axis_wavelength="wavelength",
     axis_scene_xy=("scene_x", "scene_y"),
diff --git a/ctis/inverters/__init__.py b/ctis/inverters/__init__.py
new file mode 100644
index 0000000..f2f8514
--- /dev/null
+++ b/ctis/inverters/__init__.py
@@ -0,0 +1,20 @@
+"""Inversion algorithms which can reconstruct scenes from observed images."""
+
+from . import merit
+from ._results import AbstractInversionResult, InversionResult
+from ._inverters import AbstractInverter
+from ._iterative import (
+    AbstractIterativeInverter,
+    MartInverter,
+    IterativeInversionResult,
+)
+
+__all__ = [
+    "merit",
+    "AbstractInverter",
+    "AbstractIterativeInverter",
+    "MartInverter",
+    "AbstractInversionResult",
+    "InversionResult",
+    "IterativeInversionResult",
+]
diff --git a/ctis/inverters/_inverters.py b/ctis/inverters/_inverters.py
new file mode 100644
index 0000000..0331b51
--- /dev/null
+++ b/ctis/inverters/_inverters.py
@@ -0,0 +1,48 @@
+import abc
+import dataclasses
+import named_arrays as na
+import ctis
+from ._results import InversionResult
+
+__all__ = [
+    "AbstractInverter",
+]
+
+
+@dataclasses.dataclass
+class AbstractInverter(
+    abc.ABC,
+):
+    """
+    An interface describing an algorithm which can invert CTIS observations
+    to yield a reconstruction of the observed scene.
+    """
+
+    @property
+    @abc.abstractmethod
+    def instrument(self) -> ctis.instruments.AbstractInstrument:
+        """
+        A model of a CTIS instrument which transforms the radiance of an observed
+        scene to photons measured by the sensors.
+        """
+
+    @abc.abstractmethod
+    def __call__(
+        self,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> InversionResult:
+        """
+        Reconstruct a scene using the observed images.
+
+        Parameters
+        ----------
+        images
+            The observed images used to calculate the reconstruction.
+            Must be evaluated on the same coordinates as
+            :attr:`~ctis.instruments.AbstractInstrument.coordinates_sensor`
+            attribute of :attr:`instrument`.
+        kwargs
+            Additional keyword arguments which can be used by subclass
+            implementations.
+        """
diff --git a/ctis/inverters/_inverters_test.py b/ctis/inverters/_inverters_test.py
new file mode 100644
index 0000000..b7232a4
--- /dev/null
+++ b/ctis/inverters/_inverters_test.py
@@ -0,0 +1,31 @@
+import abc
+import numpy as np
+import named_arrays as na
+import ctis
+
+
+class AbstractTestAbstractInverter(
+    abc.ABC,
+):
+
+    def test_instrument(self, a: ctis.inverters.AbstractInverter):
+        result = a.instrument
+        assert isinstance(result, ctis.instruments.AbstractInstrument)
+
+    def test__call__(
+        self,
+        a: ctis.inverters.AbstractInverter,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> ctis.inverters.AbstractInversionResult:
+        result = a(images, **kwargs)
+
+        assert isinstance(result, ctis.inverters.AbstractInversionResult)
+
+        assert result.solution.sum() > 0
+        assert isinstance(result.success, bool)
+        assert isinstance(result.message, str)
+        assert np.all(result.images == images)
+        assert result.inverter == a
+
+        return result
diff --git a/ctis/inverters/_iterative/__init__.py b/ctis/inverters/_iterative/__init__.py
new file mode 100644
index 0000000..e301a20
--- /dev/null
+++ b/ctis/inverters/_iterative/__init__.py
@@ -0,0 +1,8 @@
+from ._iterative import AbstractIterativeInverter, IterativeInversionResult
+from ._mart import MartInverter
+
+__all__ = [
+    "AbstractIterativeInverter",
+    "IterativeInversionResult",
+    "MartInverter",
+]
diff --git a/ctis/inverters/_iterative/_iterative.py b/ctis/inverters/_iterative/_iterative.py
new file mode 100644
index 0000000..9d3a0f5
--- /dev/null
+++ b/ctis/inverters/_iterative/_iterative.py
@@ -0,0 +1,150 @@
+from typing import ClassVar
+import dataclasses
+import named_arrays as na
+import ctis
+from .. import AbstractInverter, AbstractInversionResult
+
+__all__ = [
+    "AbstractIterativeInverter",
+    "IterativeInversionResult",
+]
+
+
+@dataclasses.dataclass
+class AbstractIterativeInverter(
+    AbstractInverter,
+):
+    """
+    An abstract inversion algorithm which reconstructs an observed scene
+    using iterative methods.
+
+    These methods will apply some operation repeatedly until a specified
+    convergence criteria is met.
+    """
+
+    axis_iteration: ClassVar[str] = "iteration"
+    """The logical axis associated with changing iteration index."""
+
+    num_iteration: int = dataclasses.field(default=100, kw_only=True)
+    """
+    The maximum number of iterations to perform.
+
+    If convergence is not reached before this number is exceeded,
+    a warning is raised and an unsuccessful result is returned. 
+    """
+
+    intermediate: bool = dataclasses.field(default=False, kw_only=True)
+    """
+    Whether to save intermediate solutions.
+
+    This is set to :obj:`False` during normal operation, but can be useful for
+    debugging or demonstration purposes.
+    """
+
+    def mean_chi_squared(
+        self,
+        images_observed: na.ScalarArray,
+        images_predicted: na.ScalarArray,
+    ) -> na.ScalarArray:
+        r"""
+        Evaluate :math:`\langle \chi^2 \rangle` for each observed/predicted
+        image pair.
+
+        Parameters
+        ----------
+        images_observed
+            The actual measured images.
+        images_predicted
+            The images predicted by the inversion.
+        """
+
+        uncertainty = self.instrument.uncertainty(images_predicted)
+
+        return ctis.inverters.merit.mean_chi_squared(
+            observed=images_observed,
+            expected=images_predicted,
+            uncertainty=uncertainty,
+            axis=self.instrument.axis_sensor_xy,
+        )
+
+    def correlation_residual(
+        self,
+        images_observed: na.ScalarArray,
+        images_predicted: na.ScalarArray,
+    ) -> na.ScalarArray:
+        """
+        Evaluate the correlation between the predicted images and the residual.
+
+        Parameters
+        ----------
+        images_observed
+            The actual measured images.
+        images_predicted
+            The images predicted by the inversion.
+        """
+        return ctis.inverters.merit.correlation_residual(
+            observed=images_observed,
+            expected=images_predicted,
+            axis=self.instrument.axis_sensor_xy,
+        )
+
+
+@dataclasses.dataclass
+class IterativeInversionResult(
+    AbstractInversionResult,
+):
+    """The results of an iterative inversion attempt."""
+
+    solutions: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray]
+    """
+    Intermediate solutions from each iteration.
+    
+    If :attr:`AbstractIterativeInverter.intermediate` is set to :obj:`True`,
+    this has up to :attr:`~AbstractIterativeInverter.num_iteration` elements
+    along the :attr:`~AbstractIterativeInverter.axis_iteration` logical axis.
+    Otherwise this has only one element along the
+    :attr:`~AbstractIterativeInverter.axis_iteration` axis.
+    """
+
+    success: bool = dataclasses.MISSING
+    """A boolean flag indicating whether the inversion was successful."""
+
+    images: na.FunctionArray[
+        na.SpectralPositionalVectorArray,
+        na.ScalarArray,
+    ] = dataclasses.MISSING
+    """The observed images on which the inversion was performed."""
+
+    inverter: "ctis.inverters.AbstractInverter" = dataclasses.MISSING
+    """The inversion algorithm instance that produced these results."""
+
+    message: str = dataclasses.MISSING
+    """Any message from the inversion routine concerning the results."""
+
+    num_iteration: int = dataclasses.MISSING
+    """The number of iterations performed by the inverter."""
+
+    mean_chi_squared: na.ScalarArray = dataclasses.MISSING
+    """The mean chi squared statistic for each iteration."""
+
+    correlation_residual: na.ScalarArray = dataclasses.MISSING
+    """
+    The correlation between the predicted images and the residuals 
+    for each iteration.
+    """
+
+    @property
+    def solution(
+        self,
+    ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray]:
+        axis_iteration = self.inverter.axis_iteration
+        return self.solutions[{axis_iteration: ~0}]
+
+    @property
+    def iteration(self) -> na.ScalarArray:
+        """The iteration value for each iteration."""
+        return na.arange(
+            start=0,
+            stop=self.num_iteration,
+            axis=self.inverter.axis_iteration,
+        )
diff --git a/ctis/inverters/_iterative/_iterative_test.py b/ctis/inverters/_iterative/_iterative_test.py
new file mode 100644
index 0000000..43fa9a6
--- /dev/null
+++ b/ctis/inverters/_iterative/_iterative_test.py
@@ -0,0 +1,34 @@
+import ctis
+import named_arrays as na
+from .._inverters_test import AbstractTestAbstractInverter
+
+
+class AbstractTestAbstractIterativeInverter(
+    AbstractTestAbstractInverter,
+):
+
+    def test__call__(
+        self,
+        a: ctis.inverters.AbstractIterativeInverter,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> ctis.inverters.IterativeInversionResult:
+
+        result = super().test__call__(
+            a=a,
+            images=images,
+            **kwargs,
+        )
+
+        axis_iteration = result.inverter.axis_iteration
+
+        assert result.iteration.size == result.num_iteration
+        assert result.mean_chi_squared.shape[axis_iteration] == result.num_iteration
+        assert result.correlation_residual.shape[axis_iteration] == result.num_iteration
+
+        return result
+
+    def test_num_iteration(self, a: ctis.inverters.AbstractIterativeInverter):
+        result = a.num_iteration
+        assert isinstance(result, int)
+        assert result > 0
diff --git a/ctis/inverters/_iterative/_mart/__init__.py b/ctis/inverters/_iterative/_mart/__init__.py
new file mode 100644
index 0000000..0b70d70
--- /dev/null
+++ b/ctis/inverters/_iterative/_mart/__init__.py
@@ -0,0 +1,5 @@
+from ._mart import MartInverter
+
+__all__ = [
+    "MartInverter",
+]
diff --git a/ctis/inverters/_iterative/_mart/_mart.py b/ctis/inverters/_iterative/_mart/_mart.py
new file mode 100644
index 0000000..2124913
--- /dev/null
+++ b/ctis/inverters/_iterative/_mart/_mart.py
@@ -0,0 +1,206 @@
+import warnings
+import dataclasses
+import numpy as np
+import named_arrays as na
+import ctis
+from .. import AbstractIterativeInverter, IterativeInversionResult
+
+__all__ = [
+    "MartInverter",
+]
+
+
+@dataclasses.dataclass
+class MartInverter(
+    AbstractIterativeInverter,
+):
+    """
+    An inversion routine based on the multiplicative algebraic reconstruction
+    technique (MART) :cite:t:`Gordon1970`.
+
+    For further information, see the discussion :doc:`../discussions/mart-discussion`.
+    """
+
+    instrument: ctis.instruments.AbstractInstrument = dataclasses.MISSING
+    """
+    A model of a CTIS instrument which transforms the radiance of an observed
+    scene to photons measured by the sensors.
+    """
+
+    gamma: None | float = None
+    r"""
+    Learning rate, :math:`\gamma`.
+    
+    At every iteration, the current correction, :math:`C`, is replaced by 
+    :math:`C^\gamma`.
+    
+    If :obj:`None`, :math:`\gamma = 2 / N`, where :math:`N` is the number of
+    channels.
+    """
+
+    threshold_convergence: float = 1e-3
+    r"""
+    The convergence threshold, :math:`T`, which halts the iteration.
+    
+    If :math:`\langle \chi_{i-1}^2 \rangle - \langle \chi_{i}^2 \rangle < T`,
+    then the algorithm is considered to be converged.
+    """
+
+    def __post_init__(self):
+
+        if self.gamma is None:
+            self.gamma = 2 / self.instrument.num_channel
+
+    def __call__(
+        self,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        guess: None | na.ScalarArray = None,
+        verbose: bool = False,
+    ) -> IterativeInversionResult:
+        """
+        Reconstruct a scene using the observed images.
+
+        Parameters
+        ----------
+        images
+            The observed images used to calculate the reconstruction.
+            Must be evaluated on the same position coordinates as
+            :attr:`~ctis.instruments.AbstractInstrument.coordinates_sensor`
+            attribute of :attr:`instrument`.
+        guess
+            The initial guess at the reconstructed scene.
+            Must be evaluated on the same coordinates as
+            :attr:`~ctis.instruments.AbstractInstrument.coordinates_scene`
+            attribute of :attr:`instrument`.
+        """
+
+        instrument = self.instrument
+
+        axis_channel = instrument.axis_channel
+
+        position_images = images.inputs.position
+        position_sensor = instrument.coordinates_sensor.position
+        if not np.all(position_images == position_sensor):
+            raise ValueError(
+                "`images.inputs.position` and `self.coordinates_sensor.position` "
+                "are not equal."
+            )
+        images_inputs = images.inputs
+        images = images.outputs
+
+        if guess is None:
+            scene = instrument.backproject(images).outputs
+            scene = scene.mean(axis_channel)
+            scene.ndarray[:] = scene.ndarray.mean()
+        else:
+            scene = guess.copy()
+
+        num_channel = instrument.num_channel
+
+        gamma = self.gamma
+
+        backprojected = instrument.backproject(images).outputs
+
+        backprojected = np.maximum(backprojected, 0)
+
+        intermediate = []
+
+        merit_old = np.inf
+
+        chi2 = []
+        correlation_residual = []
+
+        for i in range(self.num_iteration):
+
+            if self.intermediate:
+                intermediate.append(scene)
+
+            if verbose:  # pragma: nocover
+                print(f"{i=}")
+
+            images_new = instrument.image(scene, noise=False).outputs
+
+            chi2_ij = self.mean_chi_squared(images, images_new)
+            r_ij = self.correlation_residual(images, images_new)
+
+            chi2.append(chi2_ij)
+            correlation_residual.append(r_ij)
+
+            merit = chi2_ij.mean(axis_channel)
+
+            if verbose:  # pragma: nocover
+                print(f"merit: {merit}")
+
+            if merit > merit_old:  # pragma: nocover
+                message = "Failure: merit increasing."
+                success = False
+                num_iteration = i + 1
+                warnings.warn(message)
+                break
+
+            elif (merit_old - merit) < self.threshold_convergence:
+                message = f"Achieved merit less than {self.threshold_convergence}."
+                success = True
+                num_iteration = i + 1
+                break
+
+            backprojected_new = instrument.backproject(images_new).outputs
+
+            backprojected_new = np.maximum(backprojected_new, 0)
+
+            correction = backprojected / backprojected_new
+
+            correction = np.nan_to_num(
+                x=correction,
+                nan=1,
+                posinf=1,
+                neginf=1,
+            )
+
+            correction = correction**gamma
+
+            correction = np.prod(correction, axis=instrument.axis_channel)
+            correction = correction ** (1 / num_channel)
+
+            if self.intermediate:
+                scene = scene * correction
+            else:
+                scene *= correction
+
+            merit_old = merit
+
+        else:
+            message = f"Max number of iterations ({self.num_iteration}) exceeded."
+            warnings.warn(message)
+            success = False
+            num_iteration = self.num_iteration
+
+        if self.intermediate:
+            intermediate = na.stack(intermediate, axis=self.axis_iteration)
+            solutions = intermediate
+        else:
+            solutions = scene.add_axes(self.axis_iteration)
+
+        solutions = na.FunctionArray(
+            inputs=self.instrument.coordinates_scene,
+            outputs=solutions,
+        )
+
+        images = na.FunctionArray(
+            inputs=images_inputs,
+            outputs=images,
+        )
+
+        mean_chi_squared = na.stack(chi2, axis=self.axis_iteration)
+        correlation_residual = na.stack(correlation_residual, axis=self.axis_iteration)
+
+        return IterativeInversionResult(
+            solutions=solutions,
+            success=success,
+            images=images,
+            inverter=self,
+            message=message,
+            num_iteration=num_iteration,
+            mean_chi_squared=mean_chi_squared,
+            correlation_residual=correlation_residual,
+        )
diff --git a/ctis/inverters/_iterative/_mart/_mart_test.py b/ctis/inverters/_iterative/_mart/_mart_test.py
new file mode 100644
index 0000000..02d6d0e
--- /dev/null
+++ b/ctis/inverters/_iterative/_mart/_mart_test.py
@@ -0,0 +1,116 @@
+import matplotlib.pyplot as plt
+import pytest
+import astropy.units as u
+import named_arrays as na
+import ctis
+from .._iterative_test import AbstractTestAbstractIterativeInverter
+
+velocity = na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s
+
+wavelength_rest = 171 * u.AA
+
+AA = dict(unit=u.AA, equivalencies=u.doppler_optical(wavelength_rest))
+
+wavelength = velocity.to(**AA)
+
+position_scene = na.Cartesian2dVectorLinearSpace(
+    start=-10 * u.arcsec,
+    stop=10 * u.arcsec,
+    axis=na.Cartesian2dVectorArray("scene_x", "scene_y"),
+    num=na.Cartesian2dVectorArray(64, 64),
+)
+
+position_sensor = na.Cartesian2dVectorArray(
+    x=na.arange(0, 128 + 1, axis="sensor_x") * u.pix,
+    y=na.arange(0, 64 + 1, axis="sensor_y") * u.pix,
+)
+
+coordinates_scene = na.DopplerPositionalVectorArray.from_velocity(
+    velocity=velocity,
+    wavelength_rest=wavelength_rest,
+    position=position_scene,
+)
+coordinates_sensor = na.DopplerPositionalVectorArray.from_velocity(
+    velocity=velocity,
+    wavelength_rest=wavelength_rest,
+    position=position_sensor,
+)
+
+scene = ctis.scenes.gaussians(coordinates_scene)
+
+coordinates_scene.wavelength = wavelength
+coordinates_sensor.wavelength = wavelength
+
+angle = na.linspace(0, 360, num=4, axis="channel", endpoint=False) * u.deg
+
+instrument = ctis.instruments.IdealInstrument(
+    area_effective=1 * u.cm**2,
+    timedelta_exposure=20 * u.s,
+    plate_scale=0.4 * u.arcsec / u.pix,
+    dispersion=((10 * u.km / u.s).to(**AA) - wavelength_rest) / u.pix,
+    angle=angle,
+    wavelength_ref=wavelength_rest,
+    position_ref=na.Cartesian2dVectorArray(64, 32) * u.pix,
+    coordinates_scene=coordinates_scene,
+    coordinates_sensor=coordinates_sensor,
+    channel=angle,
+    axis_channel="channel",
+    axis_wavelength="wavelength",
+    axis_scene_xy=("scene_x", "scene_y"),
+    axis_sensor_xy=("sensor_x", "sensor_y"),
+)
+
+images = instrument.image(scene)
+
+inverter = ctis.inverters.MartInverter(
+    instrument=instrument,
+)
+
+
+@pytest.mark.parametrize(
+    argnames="a",
+    argvalues=[
+        ctis.inverters.MartInverter(
+            instrument=instrument,
+            threshold_convergence=1e-2,
+        ),
+        ctis.inverters.MartInverter(
+            instrument=instrument,
+            num_iteration=2,
+            threshold_convergence=1e-2,
+        ),
+        ctis.inverters.MartInverter(
+            instrument=instrument,
+            intermediate=True,
+            threshold_convergence=1e-2,
+        ),
+    ],
+)
+class TestMartInverter(
+    AbstractTestAbstractIterativeInverter,
+):
+    @pytest.mark.parametrize("images", [images])
+    @pytest.mark.parametrize(
+        argnames="guess",
+        argvalues=[
+            None,
+            na.ScalarArray.ones(scene.outputs.shape) * scene.outputs.unit,
+        ],
+    )
+    def test__call__(
+        self,
+        a: ctis.inverters.AbstractInverter,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        guess: na.ScalarArray,
+    ):
+        result = super().test__call__(
+            a=a,
+            images=images,
+            guess=guess,
+        )
+
+        fig, axs = result.plot_moments(scene)
+
+        assert isinstance(fig, plt.Figure)
+        for ax in axs:
+            assert isinstance(ax, plt.Axes)
diff --git a/ctis/inverters/_results.py b/ctis/inverters/_results.py
new file mode 100644
index 0000000..c5378ad
--- /dev/null
+++ b/ctis/inverters/_results.py
@@ -0,0 +1,280 @@
+import abc
+import dataclasses
+import numpy as np
+import matplotlib.pyplot as plt
+import astropy.units as u
+import astropy.visualization
+import named_arrays as na
+import ctis
+
+__all__ = [
+    "InversionResult",
+]
+
+
+@dataclasses.dataclass
+class AbstractInversionResult(
+    abc.ABC,
+):
+    """An interface describing the results of an inversion attempt."""
+
+    @property
+    @abc.abstractmethod
+    def solution(
+        self,
+    ) -> na.FunctionArray[
+        na.AbstractDopplerPositionalVectorArray,
+        na.ScalarArray,
+    ]:
+        """The reconstructed scene found by the inversion."""
+
+    @property
+    @abc.abstractmethod
+    def success(self) -> bool:
+        """Whether the inversion was successful."""
+
+    @property
+    @abc.abstractmethod
+    def images(
+        self,
+    ) -> na.FunctionArray[
+        na.AbstractDopplerPositionalVectorArray,
+        na.ScalarArray,
+    ]:
+        """The observed images used to calculate the inversion."""
+
+    @property
+    @abc.abstractmethod
+    def inverter(self) -> "ctis.inverters.AbstractInverter":
+        """The inversion algorithm instance that produced these results."""
+
+    @property
+    @abc.abstractmethod
+    def message(self) -> str:
+        """Any message from the inverter regarding these results."""
+
+    def plot_moments(
+        self,
+        truth: na.FunctionArray[
+            na.AbstractDopplerPositionalVectorArray,
+            na.ScalarArray,
+        ],
+        num_bins: int = 50,
+        range_radiance: None | tuple[u.Quantity, u.Quantity] = None,
+        range_median: None | tuple[u.Quantity, u.Quantity] = None,
+        range_iqr: None | tuple[u.Quantity, u.Quantity] = None,
+    ) -> tuple[plt.Figure, np.ndarray]:
+        recon = self.solution
+
+        axis_wavelength = self.inverter.instrument.axis_wavelength
+
+        wavelength_truth = truth.inputs.wavelength
+        wavelength_recon = recon.inputs.wavelength
+
+        dw_truth = wavelength_truth.volume_cell(axis_wavelength)
+        dw_recon = wavelength_recon.volume_cell(axis_wavelength)
+
+        radiance_truth = (truth.outputs * dw_truth).sum(axis_wavelength)
+        radiance_recon = (recon.outputs * dw_recon).sum(axis_wavelength)
+
+        median_truth = na.pdf.median(
+            x=truth.inputs.velocity,
+            f=truth.outputs,
+            axis="wavelength",
+        )
+        median_recon = na.pdf.median(
+            x=recon.inputs.velocity,
+            f=recon.outputs,
+            axis="wavelength",
+        )
+
+        iqr_truth = na.pdf.iqr(
+            x=truth.inputs.velocity,
+            f=truth.outputs,
+            axis="wavelength",
+        )
+        iqr_recon = na.pdf.iqr(
+            x=recon.inputs.velocity,
+            f=recon.outputs,
+            axis="wavelength",
+        )
+
+        bins = dict(true=num_bins, reconstructed=num_bins)
+
+        if range_radiance is None:
+            range_radiance = (None, None)
+
+        if range_median is None:
+            range_median = (None, None)
+
+        if range_iqr is None:
+            range_iqr = (None, None)
+
+        min_radiance, max_radiance = range_radiance
+        min_median, max_median = range_median
+        min_iqr, max_iqr = range_iqr
+
+        if min_radiance is None:
+            min_radiance = 0 * radiance_truth.unit
+        if max_radiance is None:
+            max_radiance = radiance_truth.max()
+
+        if min_median is None:
+            min_median = np.nanmin(median_truth)
+        if max_median is None:
+            max_median = np.nanmax(median_truth)
+
+        if min_iqr is None:
+            min_iqr = 0 * iqr_truth.unit
+        if max_iqr is None:
+            max_iqr = iqr_truth.max()
+
+        hist_radiance = na.histogram2d(
+            radiance_truth,
+            radiance_recon,
+            bins=bins,
+            min=min_radiance,
+            max=max_radiance,
+        )
+        hist_median = na.histogram2d(
+            median_truth,
+            median_recon,
+            bins=bins,
+            min=min_median,
+            max=max_median,
+        )
+        hist_iqr = na.histogram2d(
+            iqr_truth,
+            iqr_recon,
+            bins=bins,
+            min=min_iqr,
+            max=max_iqr,
+        )
+
+        hist_radiance = hist_radiance / hist_radiance.sum("reconstructed")
+        hist_median = hist_median / hist_median.sum("reconstructed")
+        hist_iqr = hist_iqr / hist_iqr.sum("reconstructed")
+
+        hist_radiance.outputs = np.nan_to_num(
+            x=hist_radiance.outputs,
+            posinf=0,
+            neginf=0,
+        )
+        hist_median.outputs = np.nan_to_num(hist_median.outputs)
+        hist_iqr.outputs = np.nan_to_num(hist_iqr.outputs)
+
+        with astropy.visualization.quantity_support():
+            fig, axs = plt.subplots(
+                constrained_layout=True,
+                figsize=(10, 4),
+                ncols=3,
+            )
+            ax_radiance, ax_median, ax_iqr = axs
+            img_radiance = na.plt.pcolormesh(
+                C=hist_radiance,
+                ax=ax_radiance,
+                vmax=np.nanpercentile(hist_radiance.outputs, 99.5),
+            )
+            img_median = na.plt.pcolormesh(
+                C=hist_median,
+                ax=ax_median,
+                vmax=np.nanpercentile(hist_median.outputs, 99.5),
+            )
+            img_iqr = na.plt.pcolormesh(
+                C=hist_iqr,
+                ax=ax_iqr,
+                vmax=np.nanpercentile(hist_iqr.outputs, 99.5),
+            )
+
+            pt_radiance = np.nanmean(radiance_truth).ndarray.value
+            pt_median = np.nanmean(median_truth).ndarray.value
+            pt_iqr = np.nanmean(iqr_truth).ndarray.value
+
+            ax_radiance.axline(
+                (pt_radiance, pt_radiance),
+                slope=1,
+                color="tab:red",
+                linestyle="dashed",
+            )
+            ax_median.axline(
+                (pt_median, pt_median),
+                slope=1,
+                color="tab:red",
+                linestyle="dashed",
+            )
+            ax_iqr.axline(
+                (pt_iqr, pt_iqr),
+                slope=1,
+                color="tab:red",
+                linestyle="dashed",
+            )
+            plt.colorbar(
+                img_radiance.ndarray.item(),
+                ax=ax_radiance,
+                location="top",
+                label="probability",
+            )
+            plt.colorbar(
+                img_median.ndarray.item(),
+                ax=ax_median,
+                location="top",
+                label="probability",
+            )
+            plt.colorbar(
+                img_iqr.ndarray.item(),
+                ax=ax_iqr,
+                location="top",
+                label="probability",
+            )
+            ax_radiance.set_xlabel(
+                f"true radiance ({ax_radiance.get_xlabel()})",
+            )
+            ax_radiance.set_ylabel(
+                f"reconstructed radiance ({ax_radiance.get_ylabel()})",
+            )
+            ax_median.set_xlabel(
+                f"true median ({ax_median.get_xlabel()})",
+            )
+            ax_median.set_ylabel(
+                f"reconstructed median ({ax_median.get_ylabel()})",
+            )
+            ax_iqr.set_xlabel(
+                f"true IQR ({ax_iqr.get_xlabel()})",
+            )
+            ax_iqr.set_ylabel(
+                f"reconstructed IQR ({ax_iqr.get_ylabel()})",
+            )
+
+            ax_radiance.set_aspect("equal")
+            ax_median.set_aspect("equal")
+            ax_iqr.set_aspect("equal")
+
+        return fig, axs
+
+
+@dataclasses.dataclass
+class InversionResult(
+    AbstractInversionResult,
+):
+    """
+    The results of an inversion attempt.
+    """
+
+    solution: na.FunctionArray[
+        na.AbstractDopplerPositionalVectorArray, na.ScalarArray
+    ] = dataclasses.MISSING
+    """The reconstructed scene found by the inversion."""
+
+    success: bool = dataclasses.MISSING
+    """A boolean flag indicating whether the inversion was successful."""
+
+    images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray] = (
+        dataclasses.MISSING
+    )
+    """The observed images on which the inversion was performed."""
+
+    inverter: "ctis.inverters.AbstractInverter" = dataclasses.MISSING
+    """The inversion algorithm instance that produced these results."""
+
+    message: str = dataclasses.MISSING
+    """Any message from the inversion routine concerning the results."""
diff --git a/ctis/inverters/merit/__init__.py b/ctis/inverters/merit/__init__.py
new file mode 100644
index 0000000..901daf4
--- /dev/null
+++ b/ctis/inverters/merit/__init__.py
@@ -0,0 +1,11 @@
+"""Functions used to evaluate the quality of CTIS inversions."""
+
+from ._merit import (
+    mean_chi_squared,
+    correlation_residual,
+)
+
+__all__ = [
+    "mean_chi_squared",
+    "correlation_residual",
+]
diff --git a/ctis/inverters/merit/_merit.py b/ctis/inverters/merit/_merit.py
new file mode 100644
index 0000000..5c73a7f
--- /dev/null
+++ b/ctis/inverters/merit/_merit.py
@@ -0,0 +1,68 @@
+from typing import Sequence
+import numpy as np
+import named_arrays as na
+
+__all__ = [
+    "mean_chi_squared",
+    "correlation_residual",
+]
+
+
+def mean_chi_squared(
+    observed: na.ScalarArray,
+    expected: na.ScalarArray,
+    uncertainty: na.ScalarArray,
+    axis: None | str | Sequence[str] = None,
+) -> na.ScalarArray:
+    r"""
+    Compute :math:`\langle \chi^2 \rangle = \biggl\langle \left( \frac{O - E}{\sigma} \right)^2 \biggr \rangle` ,
+    where :math:`O` is the observed value,
+    :math:`E` is the expected value,
+    and :math:`\sigma` denotes the standard deviation of the uncertainty.
+
+    Parameters
+    ----------
+    observed
+        The measured values.
+    expected
+        The values predicted by the model.
+    uncertainty
+        The uncertainty of the values predicted by the model.
+    axis
+        The logical axis or axes over which to average the result.
+    """
+    chisq = np.square((observed - expected) / uncertainty)
+
+    where = uncertainty != 0
+
+    return np.mean(
+        a=chisq,
+        axis=axis,
+        where=where,
+    )
+
+
+def correlation_residual(
+    observed: na.ScalarArray,
+    expected: na.ScalarArray,
+    axis: None | str | Sequence[str] = None,
+) -> na.ScalarArray:
+    """
+    Compute Pearson's correlation coefficient between the expected values
+    and the residual.
+
+    Parameters
+    ----------
+    observed
+        The measured values.
+    expected
+        The values predicted by the model.
+    axis
+        The logical axis or axes over which to average the result.
+    """
+
+    residual = observed - expected
+
+    r = na.stats.pearsonr(expected, residual, axis=axis)
+
+    return r
diff --git a/ctis/scenes/_gaussians.py b/ctis/scenes/_gaussians.py
index ddd910c..5eff2ac 100644
--- a/ctis/scenes/_gaussians.py
+++ b/ctis/scenes/_gaussians.py
@@ -14,10 +14,10 @@
 
 
 def _gaussian(
-    inputs: na.AbstractSpectralPositionalVectorArray,
+    inputs: na.AbstractDopplerPositionalVectorArray,
     amplitude: u.Quantity | na.AbstractScalar,
-    center: na.AbstractSpectralPositionalVectorArray,
-    width: na.AbstractSpectralPositionalVectorArray,
+    center: na.AbstractDopplerPositionalVectorArray,
+    width: na.AbstractDopplerPositionalVectorArray,
 ) -> na.AbstractScalar:
     """
     Compute a Gaussian kernel.
@@ -41,14 +41,30 @@ def _gaussian(
 
     inputs = inputs.cell_centers()
 
+    inputs = na.SpectralPositionalVectorArray(
+        wavelength=inputs.velocity,
+        position=inputs.position,
+    )
+
+    center = na.SpectralPositionalVectorArray(
+        wavelength=center.velocity,
+        position=center.position,
+    )
+
+    width = na.SpectralPositionalVectorArray(
+        wavelength=width.velocity,
+        position=width.position,
+    )
+
     arg = -np.square(((inputs - center) / width).length) / 2
+
     return amplitude * np.exp(arg)
 
 
 def gaussians(
-    inputs: SpectralPositionalVectorT,
-    width: na.AbstractSpectralPositionalVectorArray,
-) -> na.FunctionArray[SpectralPositionalVectorT, na.ScalarArray]:
+    inputs: na.DopplerPositionalVectorArray,
+    width: None | na.DopplerPositionalVectorArray = None,
+) -> na.FunctionArray[na.DopplerPositionalVectorArray, na.ScalarArray]:
     r"""
     A scene with eight randomly-positioned Gaussian kernels originally
     prepared by Amy R. Winebarger.
@@ -59,6 +75,9 @@ def gaussians(
         The grid of wavelengths and positions on which to evaluate the scene.
     width
         The standard deviation of the Gaussian kernels.
+        If :obj:`None` (the default), the width will be :math:`30\,\text{km/s}
+        in the wavelength direction and :math:`1\,\text{arcsec}` in the spatial
+        directions.
 
     Examples
     --------
@@ -79,8 +98,9 @@ def gaussians(
 
         # Define the grid of positions and velocities on which to evaluate the
         # test pattern
-        inputs = na.SpectralPositionalVectorArray(
-            wavelength=na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s,
+        inputs = na.DopplerPositionalVectorArray.from_velocity(
+            velocity=na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s,
+            wavelength_rest=171 * u.AA,
             position=na.Cartesian2dVectorLinearSpace(
                 start=-20 * platescale * u.pix,
                 stop=20 * platescale * u.pix,
@@ -89,18 +109,9 @@ def gaussians(
             ),
         )
 
-        # Define the standard deviations of the Gaussians in space and velocity
-        width = na.SpectralPositionalVectorArray(
-            wavelength=27 * u.km / u.s,
-            position=2.4 / 2.35 * u.arcsec,
-        )
-
         # Compute the scene of random Gaussians for the
         # given input grid and standard deviations
-        scene = ctis.scenes.gaussians(
-            inputs=inputs,
-            width=width,
-        )
+        scene = ctis.scenes.gaussians(inputs)
 
         # Plot the result
         with astropy.visualization.quantity_support():
@@ -174,7 +185,7 @@ def gaussians(
     center_y = np.array(center_y)
     center_v = np.array(center_v)
 
-    scale = 1 * u.erg / (u.cm**2 * u.sr * u.mAA * u.s)
+    scale = 1000 * u.erg / (u.cm**2 * u.sr * u.AA * u.s)
     amplitude = amplitude * scale
 
     axis = "event"
@@ -184,11 +195,19 @@ def gaussians(
     center_y = na.ScalarArray(center_y, axis) << u.arcsec
     center_v = na.ScalarArray(center_v, axis) << (u.km / u.s)
 
-    center = na.SpectralPositionalVectorArray(
-        wavelength=center_v,
+    center = na.DopplerPositionalVectorArray.from_velocity(
+        velocity=center_v,
+        wavelength_rest=inputs.wavelength_rest,
         position=na.Cartesian2dVectorArray(center_x, center_y),
     )
 
+    if width is None:
+        width = na.DopplerPositionalVectorArray.from_velocity(
+            velocity=30 * u.km / u.s,
+            wavelength_rest=inputs.wavelength_rest,
+            position=1 * u.arcsec,
+        )
+
     outputs = _gaussian(
         inputs=inputs,
         amplitude=amplitude,
diff --git a/ctis/scenes/_gaussians_test.py b/ctis/scenes/_gaussians_test.py
index 30efefd..5034778 100644
--- a/ctis/scenes/_gaussians_test.py
+++ b/ctis/scenes/_gaussians_test.py
@@ -8,24 +8,27 @@
 @pytest.mark.parametrize(
     argnames="inputs",
     argvalues=[
-        na.SpectralPositionalVectorArray(
-            wavelength=na.linspace(-500, 500, axis="wavelength", num=101) * u.km / u.s,
+        na.DopplerPositionalVectorArray.from_velocity(
+            velocity=na.linspace(-500, 500, axis="wavelength", num=101) * u.km / u.s,
+            wavelength_rest=171 * u.AA,
             position=na.Cartesian2dVectorLinearSpace(
                 start=-10 * u.arcsec,
                 stop=10 * u.arcsec,
                 axis=na.Cartesian2dVectorArray("x", "y"),
                 num=41,
             ),
-        )
+        ),
     ],
 )
 @pytest.mark.parametrize(
     argnames="width",
     argvalues=[
-        na.SpectralPositionalVectorArray(
-            wavelength=30 * u.km / u.s,
-            position=1 * u.arcsec,
-        )
+        None,
+        na.DopplerPositionalVectorArray.from_velocity(
+            velocity=30 * u.km / u.s,
+            wavelength_rest=304 * u.AA,
+            position=2 * u.arcsec,
+        ),
     ],
 )
 def test_gaussians(
diff --git a/docs/discussions/mart-discussion.rst b/docs/discussions/mart-discussion.rst
new file mode 100644
index 0000000..24e4000
--- /dev/null
+++ b/docs/discussions/mart-discussion.rst
@@ -0,0 +1,57 @@
+Multiplicative Algebraic Reconstruction Technique (MART)
+========================================================
+
+Algebraic reconstruction techniques (ARTs) are a classic approach to solving the computed
+tomography problem :cite:p:`Gordon1970`.
+There are two possible types of this technique: additive and multiplicative.
+For limited-angle tomography problems (such as reconstructing a scene using a CTIS),
+the multiplicative method is generally preferred due to its positivity-preserving
+properties.
+The multiplicative algebraic reconstruction technique (MART)
+has become the de-facto standard algorithm for reconstructing the
+solar transition region using
+the Multi-Order Solar EUV Spectrograph (MOSES) :cite:p:`Fox2010`
+and the EUV Snapshot Imaging Spectrograph (ESIS) :cite:p:`Parker2022`.
+
+In this package, our implementation of MART will generally follow the version
+described in :cite:t:`Parker2022`, with some slight adaptations to make it more
+work on genereal, curvilinear meshes.
+
+Vanilla MART
+------------
+
+The basic version of MART starts with an initial guess at the solution, :math:`\hat{u}_0`,
+which can be all ones, or some other informed choice.
+Given this boundary condition, we then loop through the following steps until
+the convergence criterion is reached:
+
+- Compute the images corresponding to the current guess, :math:`d_i = P \hat{u}_i`,
+  where :math:`P` is a projection operator representing the forward model of
+  a CTIS instrument, and :math:`i` is the current iteration index.
+- Compute the mean chi squared,
+  :math:`\langle \chi_i^2 \rangle = \biggl\langle \left( \frac{d_i - d}{\sigma_i} \right)^2 \biggr \rangle`,
+  where :math:`d` are the actual images measured by the CTIS, and :math:`\sigma_i`
+  is the uncertainty of the predicted images, :math:`d_i`.
+- Determine if the algorithm has converged by checking if
+  :math:`\langle \chi^2 \rangle` has stopped decreasing,
+  :math:`\langle \chi_{i-1}^2 \rangle - \langle \chi_{i}^2 \rangle < T`,
+  where :math:`T` is some threshold close to zero.
+- If convergence has not been reached, compute the correction factor for each channel,
+  :math:`C_i = \frac{P^* d}{P^* d_i}`,
+  where :math:`P^*` is a deprojection operator, similar to :math:`P^T`,
+  which spreads the intensity gathered by each CTIS channel evenly along
+  the projection direction.
+- Generate an effective correction factor for each channel,
+  :math:`C_i' = C_i^\gamma`, where :math:`0<\gamma<1` is the learning rate.
+- Find the total correction factor,
+  :math:`\overline{C}_i` by taking the geometric average of each channel's
+  correction factor.
+- Finally, generate a new guess by applying the correction factor to the current
+  guess, :math:`\hat{u}_{i+1} = \overline{C}_i \hat{u}_i`
+
+The main difference of this implementation from the one described in :cite:t:`Parker2022`
+is that there is no contrast-enhancement filtering yet.
+Another difference is that the correction factor is calculated in the coordinate system of the scene
+instead of the sensors.
+This is to allow us to conserve flux on both the forward and backward passes,
+potentially increasing the stability of the algorithm.
diff --git a/docs/index.rst b/docs/index.rst
index e6f55f6..064163c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -15,6 +15,7 @@ Some explanations of the theory behind inversion
 .. toctree::
     :maxdepth: 1
 
+    discussions/mart-discussion
     discussions/richardson-lucy-analogy/richardson-lucy-analogy
 
 Tutorials
@@ -26,6 +27,7 @@ Examples on how to use this package.
     :maxdepth: 1
 
     tutorials/ideal-instrument
+    tutorials/simple-mart
 
 API Reference
 =============
diff --git a/docs/refs.bib b/docs/refs.bib
index 2a23954..ed02d46 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -36,5 +36,46 @@ @ARTICLE{Bertero2005
        adsurl = {https://ui.adsabs.harvard.edu/abs/2005A&A...437..369B},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
-
+@article{Gordon1970,
+    title = {Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography},
+    journal = {Journal of Theoretical Biology},
+    volume = {29},
+    number = {3},
+    pages = {471-481},
+    year = {1970},
+    issn = {0022-5193},
+    doi = {https://doi.org/10.1016/0022-5193(70)90109-8},
+    url = {https://www.sciencedirect.com/science/article/pii/0022519370901098},
+    author = {Richard Gordon and Robert Bender and Gabor T. Herman},
+    abstract = {We give a new method for direct reconstruction of three-dimensional objects from a few electron micrographs taken at angles which need not exceed a range of 60 degrees. The method works for totally asymmetric objects, and requires little computer time or storage. It is also applicable to X-ray photography, and may greatly reduce the exposure compared to current methods of body-section radiography.}
+}
+@ARTICLE{Fox2010,
+       author = {{Fox}, J. Lewis and {Kankelborg}, Charles C. and {Thomas}, Roger J.},
+        title = "{A Transition Region Explosive Event Observed in He II with the MOSES Sounding Rocket}",
+      journal = {\apj},
+     keywords = {instrumentation: spectrographs, space vehicles: instruments, Sun: transition region, Sun: UV radiation, techniques: imaging spectroscopy},
+         year = 2010,
+        month = aug,
+       volume = {719},
+       number = {2},
+        pages = {1132-1143},
+          doi = {10.1088/0004-637X/719/2/1132},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2010ApJ...719.1132F},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+@ARTICLE{Parker2022,
+       author = {{Parker}, Jacob D. and {Smart}, Roy T. and {Kankelborg}, Charles and {Winebarger}, Amy and {Goldsworth}, Nelson},
+        title = "{First Flight of the EUV Snapshot Imaging Spectrograph (ESIS)}",
+      journal = {\apj},
+     keywords = {Spectroscopy, 1558},
+         year = 2022,
+        month = oct,
+       volume = {938},
+       number = {2},
+          eid = {116},
+        pages = {116},
+          doi = {10.3847/1538-4357/ac8eaa},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2022ApJ...938..116P},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
 
diff --git a/docs/tutorials/ideal-instrument.ipynb b/docs/tutorials/ideal-instrument.ipynb
index 708ae7d..1714d1b 100644
--- a/docs/tutorials/ideal-instrument.ipynb
+++ b/docs/tutorials/ideal-instrument.ipynb
@@ -279,7 +279,7 @@
     "tags": []
    },
    "source": [
-    "Define an ideal CTIS instrument, an instrument characterized by an effective area, plate scale, dispersion magnitude/angle, and the coordinates of the scene/sensor.\n",
+    "Define the angle of dispersion for each channel of our CTIS instrument.\n",
     "In this case, we'll define an instrument with 36 dispersion angles to help visualize the behavior of this instrument. "
    ]
   },
@@ -295,17 +295,49 @@
     "tags": []
    },
    "outputs": [],
+   "source": [
+    "angle = na.linspace(0, 360, num=36, axis=\"channel\", endpoint=False) * u.deg"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "18",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Define an ideal CTIS instrument, an instrument characterized by an effective area, plate scale, dispersion magnitude/angle, and the coordinates of the scene/sensor."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
    "source": [
     "instrument = ctis.instruments.IdealInstrument(\n",
     "    area_effective=(0.1 * u.mm**2).to(u.cm**2),\n",
     "    timedelta_exposure=0.001 * u.s,\n",
     "    plate_scale=.4 * u.arcsec / u.pix,\n",
     "    dispersion=((20 * u.km / u.s).to(**AA) - wavelength_rest) / u.pix,\n",
-    "    angle=na.linspace(0, 360, num=36, axis=\"channel\", endpoint=False) * u.deg,\n",
+    "    angle=angle,\n",
     "    wavelength_ref=wavelength_rest,\n",
     "    position_ref=32 * u.pix,\n",
     "    coordinates_scene=coordinates_scene,\n",
-    "    coordinates_sensor=coordinates_sensor,\n",
+    "    coordinates_sensor=coordinates_sensor,\\\n",
+    "    channel=\"dispersion angle = \" + angle.to_string_array(\"%03d\"),\n",
     "    axis_channel=\"channel\",\n",
     "    axis_wavelength=\"wavelength\",\n",
     "    axis_scene_xy=(\"scene_x\", \"scene_y\"),\n",
@@ -315,7 +347,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "18",
+   "id": "20",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -334,7 +366,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "19",
+   "id": "21",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -373,7 +405,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "20",
+   "id": "22",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -389,7 +421,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "21",
+   "id": "23",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -429,7 +461,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "22",
+   "id": "24",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -445,7 +477,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "23",
+   "id": "25",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -475,7 +507,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "24",
+   "id": "26",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -492,7 +524,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "25",
+   "id": "27",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -501,11 +533,13 @@
     "tags": []
    },
    "outputs": [],
-   "source": "image = instrument.image(scene, integrate=False)"
+   "source": [
+    "image = instrument.image(scene, integrate=False)"
+   ]
   },
   {
    "cell_type": "raw",
-   "id": "26",
+   "id": "28",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -522,7 +556,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "27",
+   "id": "29",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -539,9 +573,8 @@
     "        constrained_layout=True,\n",
     "    )\n",
     "    ax, cax = axs\n",
-    "    label = \"dispersion angle = \" + instrument.angle.to_string_array(\"%03d\")\n",
     "    ani, colorbar = na.plt.rgbmovie(\n",
-    "        label,\n",
+    "        instrument.channel,\n",
     "        image.inputs.wavelength,\n",
     "        image.inputs.position.x,\n",
     "        image.inputs.position.y,\n",
@@ -574,7 +607,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "28",
+   "id": "30",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -591,7 +624,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "29",
+   "id": "31",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -600,11 +633,13 @@
     "tags": []
    },
    "outputs": [],
-   "source": "image_sum = instrument.image(scene)"
+   "source": [
+    "image_sum = instrument.image(scene)"
+   ]
   },
   {
    "cell_type": "raw",
-   "id": "30",
+   "id": "32",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -620,7 +655,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "31",
+   "id": "33",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -629,11 +664,13 @@
     "tags": []
    },
    "outputs": [],
-   "source": "backprojected = instrument.backproject(image_sum)"
+   "source": [
+    "backprojected = instrument.backproject(image_sum)"
+   ]
   },
   {
    "cell_type": "raw",
-   "id": "32",
+   "id": "34",
    "metadata": {
     "editable": true,
     "raw_mimetype": "text/x-rst",
@@ -649,7 +686,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "33",
+   "id": "35",
    "metadata": {
     "editable": true,
     "raw_mimetype": "",
@@ -668,7 +705,7 @@
     "    )\n",
     "    ax, cax = axs\n",
     "    ani, colorbar = na.plt.rgbmovie(\n",
-    "        label,\n",
+    "        instrument.channel,\n",
     "        backprojected.inputs.wavelength,\n",
     "        backprojected.inputs.position.x,\n",
     "        backprojected.inputs.position.y,\n",
diff --git a/docs/tutorials/simple-mart.ipynb b/docs/tutorials/simple-mart.ipynb
new file mode 100644
index 0000000..51420d4
--- /dev/null
+++ b/docs/tutorials/simple-mart.ipynb
@@ -0,0 +1,986 @@
+{
+ "cells": [
+  {
+   "cell_type": "raw",
+   "id": "0",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Invert a Synthetic Scene with MART\n",
+    "==================================\n",
+    "In this tutorial, we'll use the Multiplicative Algebraic Reconstruction Technique (MART) :cite:p:`Gordon1970` to invert the :func:`~ctis.scenes.gaussians` sample scene developed by Amy R. Winebarger.\n",
+    "We will use a simple instrument with four projections to image this scene and then use the MART algorithm to reconstruct the original scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "1",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "import IPython.display\n",
+    "import matplotlib.pyplot as plt\n",
+    "import astropy.units as u\n",
+    "import astropy.visualization\n",
+    "import named_arrays as na\n",
+    "import ctis"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "2",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Start by defining a grid of Doppler velocities on which to reconstruct the scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "3",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "velocity = na.linspace(-500, 500, axis=\"wavelength\", num=21) * u.km / u.s"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "4",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Define the rest wavelength for converting between velocity and wavelength."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "5",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "wavelength_rest = 171 * u.AA"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "6",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Now define a grid of positions on which to reconstruct the scene,"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "position_scene = na.Cartesian2dVectorLinearSpace(\n",
+    "    start=-10 * u.arcsec,\n",
+    "    stop=10 * u.arcsec,\n",
+    "    axis=na.Cartesian2dVectorArray(\"scene_x\", \"scene_y\"),\n",
+    "    num=na.Cartesian2dVectorArray(64 + 1, 64 + 1),\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "8",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "and a grid of positions on the sensor representing the vertices of each pixel."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "9",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "position_sensor = na.Cartesian2dVectorArray(\n",
+    "    x=na.arange(0, 128 + 1, axis=\"sensor_x\") * u.pix,\n",
+    "    y=na.arange(0, 64 + 1, axis=\"sensor_y\") * u.pix,\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "10",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Combine the 1D velocity grid and the 2D position grid into a single 3D grid for both the scene and sensor coordinates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "11",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "coordinates_scene = na.DopplerPositionalVectorArray.from_velocity(\n",
+    "    velocity=velocity,\n",
+    "    wavelength_rest=wavelength_rest,\n",
+    "    position=position_scene,\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "12",
+   "metadata": {},
+   "source": [
+    "coordinates_sensor = na.DopplerPositionalVectorArray.from_velocity(\n",
+    "    velocity=velocity,\n",
+    "    wavelength_rest=wavelength_rest,\n",
+    "    position=position_sensor,\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "13",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Create a synthetic scene composed of spatial/spectral 3D Gaussians with various Doppler shifts."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "14",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "scene = ctis.scenes.gaussians(coordinates_scene)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "15",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Add a small background equal to 1 percent of the maximum value of the scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "16",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "scene = scene + scene.max() / 100"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "17",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Display the scene as a false-color image."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "18",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, axs = plt.subplots(\n",
+    "        ncols=2,\n",
+    "        gridspec_kw=dict(width_ratios=[.9,.1]),\n",
+    "        constrained_layout=True,\n",
+    "    )\n",
+    "    ax, cax = axs\n",
+    "    colorbar = na.plt.rgbmesh(\n",
+    "        C=scene,\n",
+    "        axis_wavelength=\"wavelength\",\n",
+    "        ax=ax,\n",
+    "        vmin=0,\n",
+    "        vmax=scene.outputs.max(),\n",
+    "    )\n",
+    "    na.plt.pcolormesh(\n",
+    "        C=colorbar,\n",
+    "        axis_rgb=\"wavelength\",\n",
+    "        ax=cax,\n",
+    "    )\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_xlabel(f\"scene $x$ ({ax.get_xlabel()})\")\n",
+    "    ax.set_ylabel(f\"scene $y$ ({ax.get_ylabel()})\")\n",
+    "    cax.xaxis.set_ticks_position(\"top\")\n",
+    "    cax.xaxis.set_label_position(\"top\")\n",
+    "    cax.yaxis.tick_right()\n",
+    "    cax.yaxis.set_label_position(\"right\")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "19",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Compute the average spectrum of the scene"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "20",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "spectrum = scene.outputs.mean((\"scene_x\", \"scene_y\"))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "21",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Plot the average spectrum of the scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "22",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(constrained_layout=True)\n",
+    "    ax2 = ax.twiny()\n",
+    "    na.plt.stairs(\n",
+    "        velocity,\n",
+    "        spectrum,\n",
+    "        ax=ax,\n",
+    "    )\n",
+    "    na.plt.stairs(\n",
+    "        scene.inputs.wavelength,\n",
+    "        spectrum,\n",
+    "        ax=ax2\n",
+    "    )\n",
+    "    ax.set_xlabel(f\"Doppler velocity ({ax.get_xlabel()})\")\n",
+    "    ax2.set_xlabel(f\"wavelength ({ax2.get_xlabel()})\")\n",
+    "    ax.set_ylabel(f\"average radiance ({ax.get_ylabel()})\")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "23",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Define the dispersion angles for our instrument.\n",
+    "In this case we'll define four channels, each each separated by :math:`90^\\circ` degrees."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "24",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "angle = na.linspace(0, 360, num=4, axis=\"channel\", endpoint=False) * u.deg + 5.64 * u.deg"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "25",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Define the magnitude of dispersion for our instrument in terms of Doppler velocity and then convert to wavelength units."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "26",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "dispersion = 10 * u.km / u.s\n",
+    "dispersion = dispersion.to(u.AA, equivalencies=u.doppler_optical(wavelength_rest))\n",
+    "dispersion = (dispersion - wavelength_rest) / u.pix\n",
+    "dispersion.to(u.mAA / u.pix)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "27",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Create an ideal CTIS using the dispersion magnitude and angles."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "28",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "instrument = ctis.instruments.IdealInstrument(\n",
+    "    area_effective=1 * u.cm ** 2,\n",
+    "    timedelta_exposure=20 * u.s,\n",
+    "    plate_scale=.4 * u.arcsec / u.pix,\n",
+    "    dispersion=dispersion,\n",
+    "    angle=angle,\n",
+    "    wavelength_ref=wavelength_rest,\n",
+    "    position_ref=na.Cartesian2dVectorArray(64, 32) * u.pix,\n",
+    "    coordinates_scene=coordinates_scene,\n",
+    "    coordinates_sensor=coordinates_sensor,\n",
+    "    channel=\"dispersion angle = \" + angle.to_string_array(\"%03d\"),\n",
+    "    axis_channel=\"channel\",\n",
+    "    axis_wavelength=\"wavelength\",\n",
+    "    axis_scene_xy=(\"scene_x\", \"scene_y\"),\n",
+    "    axis_sensor_xy=(\"sensor_x\", \"sensor_y\"),\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "29",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Apply the forward model of this instrument to the scene to calculate the observed images."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "30",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": "images = instrument.image(scene)",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "31",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Display the images as an animation, where each frame represents a different channel / dispersion direction."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "32",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(\n",
+    "        constrained_layout=True,\n",
+    "        figsize=(9.2, 4),\n",
+    "    )\n",
+    "    norm = plt.Normalize(\n",
+    "        vmin=0,\n",
+    "        vmax=images.outputs.value.ndarray.max(),\n",
+    "    )\n",
+    "    colorizer = plt.Colorizer(\n",
+    "        cmap=\"gray\",\n",
+    "        norm=norm,\n",
+    "    )\n",
+    "    ani = na.plt.pcolormovie(\n",
+    "        instrument.channel,\n",
+    "        images.inputs.position.x,\n",
+    "        images.inputs.position.y,\n",
+    "        C=images.outputs.value,\n",
+    "        axis_time=\"channel\",\n",
+    "        ax=ax,\n",
+    "        kwargs_pcolormesh=dict(\n",
+    "            colorizer=colorizer,\n",
+    "        ),\n",
+    "    )\n",
+    "    plt.colorbar(\n",
+    "        mappable=plt.cm.ScalarMappable(colorizer=colorizer),\n",
+    "        ax=ax,\n",
+    "        label=f\"signal ({images.outputs.unit:latex_inline})\",\n",
+    "    )\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_xlabel(f\"sensor $x$ ({images.inputs.position.x.unit})\")\n",
+    "    ax.set_ylabel(f\"sensor $y$ ({images.inputs.position.y.unit})\")\n",
+    "\n",
+    "result = ani.to_jshtml(fps=2)\n",
+    "result = IPython.display.HTML(result)\n",
+    "\n",
+    "plt.close(ani._fig)\n",
+    "\n",
+    "result"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "33",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Initialize the MART inversion algorithm with the instrument model. We'll also enable saving intermediate results so that we can visualize the behavior of the algorithm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "34",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "mart = ctis.inverters.MartInverter(\n",
+    "    instrument=instrument,\n",
+    "    intermediate=True,\n",
+    ")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "35",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Invert the images using our instance of MART and the initial guess."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "36",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": "inversion = mart(images)",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "37",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Display the results as a false-color movie, where each frame represents subsequent iterations of the MART algorithm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "38",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, axs = plt.subplots(\n",
+    "        ncols=3,\n",
+    "        gridspec_kw=dict(width_ratios=[.5, .5, .1]),\n",
+    "        constrained_layout=True,\n",
+    "        figsize=(10, 4.5),\n",
+    "    )\n",
+    "    ax1, ax2, cax = axs\n",
+    "    ax2.set_yticklabels([])\n",
+    "    na.plt.rgbmesh(\n",
+    "        C=scene,\n",
+    "        axis_wavelength=\"wavelength\",\n",
+    "        ax=ax1,\n",
+    "        vmin=0,\n",
+    "        vmax=scene.outputs.max(),\n",
+    "    )\n",
+    "    label = \"iteration = \" + inversion.iteration.to_string_array(\"%d\") + \"\\n\"\n",
+    "    name = r\"$\\langle \\chi^2 \\rangle$\"\n",
+    "    label = label + f\"{name} = \" + inversion.mean_chi_squared.mean(instrument.axis_channel).to_string_array()\n",
+    "    ani, colorbar = na.plt.rgbmovie(\n",
+    "        label,\n",
+    "        scene.inputs.wavelength,\n",
+    "        scene.inputs.position.x,\n",
+    "        scene.inputs.position.y,\n",
+    "        C=inversion.solutions.outputs,\n",
+    "        axis_time=inversion.inverter.axis_iteration,\n",
+    "        axis_wavelength=\"wavelength\",\n",
+    "        ax=ax2,\n",
+    "        vmin=0,\n",
+    "        vmax=scene.outputs.max(),\n",
+    "    )\n",
+    "    na.plt.pcolormesh(\n",
+    "        C=colorbar,\n",
+    "        axis_rgb=\"wavelength\",\n",
+    "        ax=cax,\n",
+    "    )\n",
+    "    ax1.set_title(\"original\")\n",
+    "    ax2.set_title(\"reconstructed\")\n",
+    "    unit_x = scene.inputs.position.x.unit\n",
+    "    unit_y = scene.inputs.position.y.unit\n",
+    "    ax1.set_xlabel(f\"scene $x$ ({unit_x:latex_inline})\")\n",
+    "    ax2.set_xlabel(f\"scene $x$ ({unit_x:latex_inline})\")\n",
+    "    ax1.set_ylabel(f\"scene $y$ ({unit_y:latex_inline})\")\n",
+    "    cax.xaxis.set_ticks_position(\"top\")\n",
+    "    cax.xaxis.set_label_position(\"top\")\n",
+    "    cax.yaxis.tick_right()\n",
+    "    cax.yaxis.set_label_position(\"right\")\n",
+    "\n",
+    "result = ani.to_jshtml(fps=20)\n",
+    "result = IPython.display.HTML(result)\n",
+    "\n",
+    "plt.close(ani._fig)\n",
+    "\n",
+    "result"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "39",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Isolate the solution array from the inversion result object."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "40",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "solution = inversion.solution"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "41",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Compute the average spectrum of the reconstructed scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "42",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "spectrum_inverted = solution.outputs.mean((\"scene_x\", \"scene_y\"))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "43",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Plot the average spectrum of the original scene vs. the average spectrum of the reconstructed scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "44",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(constrained_layout=True)\n",
+    "    na.plt.stairs(\n",
+    "        scene.inputs.wavelength,\n",
+    "        spectrum,\n",
+    "        ax=ax,\n",
+    "        label=\"original\",\n",
+    "    )\n",
+    "    na.plt.stairs(\n",
+    "        scene.inputs.wavelength,\n",
+    "        spectrum_inverted,\n",
+    "        ax=ax,\n",
+    "        label=\"reconstructed\",\n",
+    "    )\n",
+    "    ax.set_xlabel(f\"wavelength ({ax.get_xlabel()})\")\n",
+    "    ax2.set_xlabel(f\"wavelength ({ax2.get_xlabel()})\")\n",
+    "    ax.set_ylabel(f\"average radiance ({ax.get_ylabel()})\")\n",
+    "    ax.legend()"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "45",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Plot 2D histograms of the true vs. reconstructed value of the total radiance, median (Doppler shift), and interquartile range (Doppler width) for every pixel in the scene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "46",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "inversion.plot_moments(scene);"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "raw",
+   "id": "47",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "text/x-rst",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Plot :math:`\\langle \\chi^2 \\rangle` and the signal-correlated residual as a function of iteration."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "48",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "fig, ax = plt.subplots(\n",
+    "    nrows=2,\n",
+    "    sharex=True,\n",
+    "    constrained_layout=True,\n",
+    ")\n",
+    "na.plt.plot(\n",
+    "    inversion.iteration,\n",
+    "    inversion.mean_chi_squared,\n",
+    "    ax=ax[0],\n",
+    "    axis=inversion.inverter.axis_iteration,\n",
+    "    label=instrument.channel,\n",
+    ")\n",
+    "na.plt.plot(\n",
+    "    inversion.iteration,\n",
+    "    inversion.correlation_residual,\n",
+    "    ax=ax[1],\n",
+    "    axis=inversion.inverter.axis_iteration,\n",
+    "    label=instrument.channel,\n",
+    ")\n",
+    "ax[0].set_ylabel(r\"$\\langle \\chi^2 \\rangle$\")\n",
+    "ax[1].set_xlabel(\"iteration\")\n",
+    "ax[1].set_ylabel(\"signal-correlated residual\")\n",
+    "ax[0].set_yscale(\"log\")\n",
+    "ax[0].legend();"
+   ],
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
index cfbf262..851640b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,7 +17,7 @@ classifiers = [
 ]
 dependencies = [
     "astropy",
-    "named-arrays~=1.1",
+    "named-arrays~=1.4",
 ]
 dynamic = ["version"]
 
@@ -28,6 +28,7 @@ test = [
 doc = [
     "pytest",
     "matplotlib",
+    "interface-region-imaging-spectrograph",
     "graphviz",
     "sphinx-autodoc-typehints",
     "sphinxcontrib-bibtex",