diff --git a/docs/index.rst b/docs/index.rst
index 37ae8ad..bdb0e29 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -82,6 +82,7 @@ Jupyter notebook examples on how to use :mod:`optika`.
     :maxdepth: 1
 
     tutorials/prime_focus
+    tutorials/fzp_focus
 
 
 API Reference
diff --git a/docs/tutorials/fzp_focus.ipynb b/docs/tutorials/fzp_focus.ipynb
new file mode 100644
index 0000000..a87f868
--- /dev/null
+++ b/docs/tutorials/fzp_focus.ipynb
@@ -0,0 +1,762 @@
+{
+ "cells": [
+  {
+   "cell_type": "raw",
+   "id": "0",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Fresnel Zone Plate (FZP) Telescope Tutorial\n",
+    "============================================\n",
+    "\n",
+    "A `Fresnel Zone Plate (FZP) <https://en.wikipedia.org/wiki/Zone_plate>`_\n",
+    "is a diffractive optical element that focuses light by encoding the\n",
+    "interference pattern of two spherical waves on a flat surface.\n",
+    "In :mod:`optika`, an FZP is modeled as a transmissive surface with\n",
+    ":class:`~optika.rulings.HolographicRulingSpacing` rulings, where the\n",
+    "two recording-beam origins define the focal geometry.\n",
+    "\n",
+    "This tutorial builds a single-element FZP telescope that focuses\n",
+    "collimated 171 Å EUV light onto a detector and examines the resulting\n",
+    "spot diagrams."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import astropy.units as u\n",
+    "import astropy.visualization\n",
+    "import named_arrays as na\n",
+    "import optika"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "2",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "FZP\n",
+    "---\n",
+    "\n",
+    "We start by defining the aperture radius of the FZP"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "radius_aperture = 70 * u.mm"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "4",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "and the focal length"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "focal_length = 1000 * u.mm"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "6",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "the design wavelength is"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wavelength = 171 * u.AA"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "8",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Holographic recording geometry\n",
+    "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
+    "\n",
+    "An FZP is equivalent to a holographic grating whose rulings were\n",
+    "recorded by interfering two coherent beams.\n",
+    "Let :math:`f` denote the focal length.\n",
+    "To focus incoming collimated light (a plane wave from :math:`z = -\\infty`)\n",
+    "onto a focal point at :math:`z = +f`, we choose:\n",
+    "\n",
+    "- **Recording beam 1** (``x1``): originates very far away on the optical\n",
+    "  axis, representing a plane wave.  With ``is_diverging_1=True`` and\n",
+    "  :math:`x_1 \\to -\\infty`, the unit vectors from :math:`x_1` to every\n",
+    "  point on the surface become parallel, so the ruling pattern is\n",
+    "  identical to that produced by a true plane wave.  We place ``x1`` at\n",
+    "  1 AU (the Earth–Sun distance) so that the residual focal-length shift\n",
+    "  :math:`\\delta f = f^2 / d_{x_1} = (1000\\ \\mathrm{mm})^2 /\n",
+    "  (1.496 \\times 10^{14}\\ \\mathrm{mm}) \\approx 6.7\\ \\mathrm{pm}`\n",
+    "  is negligible compared to the diffraction-limited depth of focus\n",
+    "  :math:`\\sim\\lambda (f/D)^2 \\approx 32\\ \\mu\\mathrm{m}`.\n",
+    "- **Recording beam 2** (``x2``): converges *to* the desired focal point\n",
+    "  at :math:`z = +f`.  ``is_diverging_2=False`` means rays are directed\n",
+    "  toward :math:`x_2`, i.e. the second recording beam is a converging\n",
+    "  spherical wave.\n",
+    "\n",
+    "The :class:`~optika.rulings.HolographicRulingSpacing` class computes\n",
+    "the spatially-varying ruling vector :math:`\\mathbf{d}(\\mathbf{r})` from\n",
+    "these two source points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spacing = optika.rulings.HolographicRulingSpacing(\n",
+    "    x1=na.Cartesian3dVectorArray(0, 0, -1) * (1 * u.au).to(u.mm),\n",
+    "    x2=na.Cartesian3dVectorArray(\n",
+    "        x=0 * u.mm,\n",
+    "        y=0 * u.mm,\n",
+    "        z=focal_length,\n",
+    "    ),\n",
+    "    wavelength=wavelength,\n",
+    "    is_diverging_1=True,\n",
+    "    is_diverging_2=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "10",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "We wrap the spacing in an ideal :class:`~optika.rulings.Rulings` instance.\n",
+    ":class:`~optika.rulings.Rulings` assumes perfect diffraction efficiency\n",
+    "in the specified order, which is appropriate for an ideal FZP."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rulings_fzp = optika.rulings.Rulings(\n",
+    "    spacing=spacing,\n",
+    "    diffraction_order=1,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "12",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "We assemble the FZP as a flat, transmissive :class:`~optika.surfaces.Surface`.\n",
+    "No sag profile is needed (the surface is flat) and no mirror material is\n",
+    "assigned (the default :class:`~optika.materials.Vacuum` material makes the\n",
+    "surface transmissive, not reflective).\n",
+    "Setting ``is_pupil_stop=True`` marks the FZP aperture as the entrance pupil."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fzp = optika.surfaces.Surface(\n",
+    "    name=\"FZP\",\n",
+    "    aperture=optika.apertures.CircularAperture(radius_aperture),\n",
+    "    rulings=rulings_fzp,\n",
+    "    is_pupil_stop=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "14",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Sensor\n",
+    "------\n",
+    "\n",
+    "If the size of each pixel in the sensor is"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "width_pixel = 13 * u.um"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "16",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "the number of pixels along each axis is"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_pixel = na.Cartesian2dVectorArray(512, 512)"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "18",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "The sensor is placed at the focal plane, one focal length downstream\n",
+    "of the FZP.\n",
+    "Unlike a reflective system (such as the prime-focus telescope, where the\n",
+    "reflected beam travels in the :math:`-z` direction and the sensor must\n",
+    "be rotated 180°), here the beam continues in the :math:`+z` direction\n",
+    "after transmission through the FZP, so no rotation is needed.\n",
+    "Setting ``is_field_stop=True`` marks the sensor as the field stop."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sensor = optika.sensors.ImagingSensor(\n",
+    "    name=\"sensor\",\n",
+    "    width_pixel=width_pixel,\n",
+    "    axis_pixel=na.Cartesian2dVectorArray(\"detector_x\", \"detector_y\"),\n",
+    "    num_pixel=num_pixel,\n",
+    "    transformation=na.transformations.Cartesian3dTranslation(z=focal_length),\n",
+    "    is_field_stop=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20",
+   "metadata": {},
+   "source": [
+    "## Source\n",
+    "\n",
+    "We model the source as a solar disk — a circular patch of sky subtending\n",
+    "0.5° full diameter (0.25° half-angle). The aperture is specified as a\n",
+    "dimensionless cosine so that the surface is treated as being at infinity\n",
+    "(a directional source, not a physical aperture at a finite distance)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "radius_sun = 0.25 * u.deg  # 0.5 deg full angular diameter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source = optika.surfaces.Surface(\n",
+    "    name=\"solar disk\",\n",
+    "    aperture=optika.apertures.CircularAperture(\n",
+    "        radius=np.cos(radius_sun),\n",
+    "    ),\n",
+    "    transformation=na.transformations.Cartesian3dTranslation(\n",
+    "        z=-200 * u.mm,\n",
+    "    ),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "23",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Input rays\n",
+    "----------\n",
+    "\n",
+    "For a transmissive FZP (object at infinity), we specify the input rays\n",
+    "using **physical** coordinates rather than normalized ones.\n",
+    "The pupil grid gives physical positions on the FZP aperture in mm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "24",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Use an even number of samples so the grid never lands exactly at the\n",
+    "# FZP center (0, 0), where the ruling vector degenerates to zero.\n",
+    "pupil = na.Cartesian2dVectorLinearSpace(\n",
+    "    start=-radius_aperture,\n",
+    "    stop=radius_aperture,\n",
+    "    axis=na.Cartesian2dVectorArray(\"px\", \"py\"),\n",
+    "    num=10,\n",
+    "    centers=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "25",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "The field grid gives the angular directions of the incoming collimated\n",
+    "beams in degrees.\n",
+    "The half-field angle is set by the sensor size:\n",
+    ":math:`\\theta_{1/2} = \\arctan(N_\\mathrm{pix} \\, w_\\mathrm{pix} / 2 \\, f)`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "26",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "half_field = np.arctan(\n",
+    "    (num_pixel.x * width_pixel / (2 * focal_length)).to(u.dimensionless_unscaled).value\n",
+    ") * u.rad\n",
+    "half_field = half_field.to(u.deg)\n",
+    "half_field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "27",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "field = na.Cartesian2dVectorLinearSpace(\n",
+    "    start=-half_field,\n",
+    "    stop=half_field,\n",
+    "    axis=na.Cartesian2dVectorArray(\"fx\", \"fy\"),\n",
+    "    num=5,\n",
+    "    centers=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "28",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "We combine the field, pupil, and wavelength into an\n",
+    ":class:`~optika.vectors.ObjectVectorArray`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_input = optika.vectors.ObjectVectorArray(\n",
+    "    wavelength=wavelength,\n",
+    "    field=field,\n",
+    "    pupil=pupil,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "30",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Building the system\n",
+    "-------------------\n",
+    "\n",
+    ":class:`~optika.systems.SequentialSystem` assembles the optical surfaces\n",
+    "and the input ray grid.\n",
+    "No explicit ``object`` surface is provided, so the object is treated as\n",
+    "being at infinity (collimated input).\n",
+    "The FZP is the pupil stop and the sensor is the field stop."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "31",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "system = optika.systems.SequentialSystem(\n",
+    "    object=source,\n",
+    "    surfaces=[fzp],\n",
+    "    sensor=sensor,\n",
+    "    grid_input=grid_input,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "32",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "For a transmissive surface used as the pupil stop, the internal\n",
+    "normalization algorithm in\n",
+    ":class:`~optika.systems.SequentialSystem` does not correctly infer\n",
+    "the field-of-view extents from a backward raytrace.\n",
+    "We therefore compute the default ray function directly with physical\n",
+    "(non-normalized) coordinates and cache the result manually.\n",
+    "\n",
+    ".. note::\n",
+    "\n",
+    "    This workaround is required whenever a transmissive grating is the\n",
+    "    pupil stop.  Purely reflective systems (e.g. the prime-focus telescope)\n",
+    "    are not affected."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "33",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rays_default = system.rayfunction(\n",
+    "    field=field,\n",
+    "    pupil=pupil,\n",
+    "    normalized_field=False,\n",
+    "    normalized_pupil=False,\n",
+    ")\n",
+    "system.__dict__[\"rayfunction_default\"] = rays_default"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "34",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "We can plot the optical layout using\n",
+    ":meth:`~optika.systems.SequentialSystem.plot`.\n",
+    "To avoid the same normalization issue affecting the displayed rays,\n",
+    "we plot the surfaces only and overlay the ray paths from a separate\n",
+    ":meth:`~optika.systems.SequentialSystem.raytrace` call."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "35",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Use a coarser, on-axis pupil grid for the layout diagram\n",
+    "pupil_layout = na.Cartesian2dVectorLinearSpace(\n",
+    "    start=-radius_aperture,\n",
+    "    stop=radius_aperture,\n",
+    "    axis=na.Cartesian2dVectorArray(\"px\", \"py\"),\n",
+    "    num=5,\n",
+    "    centers=True,\n",
+    ")\n",
+    "\n",
+    "raytrace_layout = system.raytrace(\n",
+    "    field=na.Cartesian2dVectorArray(0 * u.deg, 0 * u.deg),\n",
+    "    pupil=pupil_layout,\n",
+    "    normalized_field=False,\n",
+    "    normalized_pupil=False,\n",
+    ")\n",
+    "\n",
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(constrained_layout=True)\n",
+    "    system.plot(\n",
+    "        ax=ax,\n",
+    "        components=(\"z\", \"y\"),\n",
+    "        color=\"black\",\n",
+    "        zorder=10,\n",
+    "        plot_rays=False,\n",
+    "    )\n",
+    "    na.plt.plot(\n",
+    "        raytrace_layout.outputs.position,\n",
+    "        ax=ax,\n",
+    "        axis=system.axis_surface,\n",
+    "        components=(\"z\", \"y\"),\n",
+    "        color=\"tab:blue\",\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "raw",
+   "id": "36",
+   "metadata": {
+    "raw_mimetype": "text/restructuredtext"
+   },
+   "source": [
+    "Spot Diagrams\n",
+    "-------------\n",
+    "\n",
+    "We can plot a spot diagram for each field angle using the\n",
+    ":meth:`~optika.systems.SequentialSystem.spot_diagram` method.\n",
+    "The on-axis field point (center of the grid) should show a tight\n",
+    "spot, while off-axis points will exhibit the coma and astigmatism\n",
+    "that are characteristic of a single on-axis zone plate."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = system.spot_diagram()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38",
+   "metadata": {},
+   "source": [
+    "## Diffraction-limited performance\n",
+    "\n",
+    "An ideal FZP focuses on-axis collimated light to a single geometric point — the only resolution limit is diffraction. The radius of the first dark ring of the Airy disk sets the diffraction limit:\n",
+    "\n",
+    "$$r_\\mathrm{Airy} = \\frac{1.22\\,\\lambda\\,f}{D}$$\n",
+    "\n",
+    "We verify this by tracing a dense on-axis pupil grid and comparing the geometric RMS spot radius to $r_\\mathrm{Airy}$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r_airy = (1.22 * wavelength * focal_length / (2 * radius_aperture)).to(u.um)\n",
+    "r_airy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "40",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Dense on-axis pupil grid — use even num to avoid (0, 0)\n",
+    "pupil_dense = na.Cartesian2dVectorLinearSpace(\n",
+    "    start=-radius_aperture,\n",
+    "    stop=radius_aperture,\n",
+    "    axis=na.Cartesian2dVectorArray(\"px_dense\", \"py_dense\"),\n",
+    "    num=50,\n",
+    "    centers=True,\n",
+    ")\n",
+    "\n",
+    "rays_onaxis = system.rayfunction(\n",
+    "    field=na.Cartesian2dVectorArray(0 * u.deg, 0 * u.deg),\n",
+    "    pupil=pupil_dense,\n",
+    "    normalized_field=False,\n",
+    "    normalized_pupil=False,\n",
+    ").outputs\n",
+    "\n",
+    "where = rays_onaxis.unvignetted.ndarray\n",
+    "pos_x = rays_onaxis.position.x.to(u.um).ndarray.value  # plain float array, in um\n",
+    "pos_y = rays_onaxis.position.y.to(u.um).ndarray.value\n",
+    "\n",
+    "centroid_x = np.mean(pos_x[where])\n",
+    "centroid_y = np.mean(pos_y[where])\n",
+    "dx = pos_x[where] - centroid_x\n",
+    "dy = pos_y[where] - centroid_y\n",
+    "\n",
+    "r_geometric = np.sqrt(np.mean(dx**2 + dy**2))  # um\n",
+    "\n",
+    "print(f\"Airy disk radius (first dark ring):  {r_airy:.4f}\")\n",
+    "print(f\"Geometric RMS spot radius (on-axis): {r_geometric:.4e} um\")\n",
+    "print(f\"Ratio (geometric / Airy):            {r_geometric / r_airy.value:.2e}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "41",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "theta = np.linspace(0, 2 * np.pi, 361)\n",
+    "x_airy = r_airy.value * np.cos(theta)\n",
+    "y_airy = r_airy.value * np.sin(theta)\n",
+    "\n",
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(constrained_layout=True)\n",
+    "    ax.scatter(dx, dy, s=4, alpha=0.4, color=\"tab:blue\",\n",
+    "               label=f\"ray positions (RMS = {r_geometric:.2e} um)\")\n",
+    "    ax.plot(x_airy, y_airy, color=\"tab:orange\", lw=2,\n",
+    "            label=f\"Airy disk  ($r$ = {r_airy:.4f})\")\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_xlabel(f\"$x$ ({u.um:latex_inline})\")\n",
+    "    ax.set_ylabel(f\"$y$ ({u.um:latex_inline})\")\n",
+    "    ax.legend()\n",
+    "    ax.set_title(\"On-axis geometric PSF vs. diffraction limit\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "42",
+   "metadata": {},
+   "source": [
+    "## Chromatic aberration\n",
+    "\n",
+    "Zone plates are notorious for strong chromatic aberration: focal length scales *inversely* with wavelength,\n",
+    "\n",
+    "$$f(\\lambda) = f_0 \\frac{\\lambda_0}{\\lambda}$$\n",
+    "\n",
+    "where $\\lambda_0 = 171$ Å and $f_0 = 1000$ mm are the design values.\n",
+    "A wavelength offset $\\Delta\\lambda$ produces a defocused disk of radius\n",
+    "\n",
+    "$$r_\\mathrm{blur} = R_\\mathrm{pupil}\\left|1 - \\frac{\\lambda}{\\lambda_0}\\right| \\approx R_\\mathrm{pupil} \\frac{|\\Delta\\lambda|}{\\lambda_0}$$\n",
+    "\n",
+    "at the design focal plane.\n",
+    "We use *asymmetric* offsets (−2 Å and +3 Å from 171 Å) so the three spots\n",
+    "have distinct radii (0, 0.82 mm, and 1.23 mm) and are visually separable."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "43",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Build the on-axis pupil grid from the 1-D coordinate arrays.\n",
+    "# Using meshgrid with indexing='ij' guarantees consistent (x, y) pairing\n",
+    "# and avoids axis-ordering ambiguity in the named_arrays output arrays.\n",
+    "px_1d = pupil_dense.x.to(u.mm).ndarray.value   # (50,)  px_dense\n",
+    "py_1d = pupil_dense.y.to(u.mm).ndarray.value   # (50,)  py_dense\n",
+    "px_2d, py_2d = np.meshgrid(px_1d, py_1d, indexing=\"ij\")  # (50, 50)\n",
+    "inside = np.sqrt(px_2d**2 + py_2d**2) <= radius_aperture.to(u.mm).value\n",
+    "\n",
+    "wl_design = wavelength.to(u.AA).value   # 171\n",
+    "f0        = focal_length.to(u.mm).value  # 1000\n",
+    "\n",
+    "# Asymmetric wavelengths so each spot has a distinct radius.\n",
+    "# Radii: 169 Å → 0.82 mm, 171 Å → 0 (design), 174 Å → 1.23 mm\n",
+    "wavelengths_plot = [169, 171, 174]\n",
+    "colors = [\"tab:blue\", \"tab:green\", \"tab:red\"]\n",
+    "labels = [\n",
+    "    f\"169 Å  ($f$ = {f0 * wl_design / 169:.0f} mm,  $r_\\\\mathrm{{blur}}$ = {70 * abs(1 - 169/wl_design):.2f} mm)\",\n",
+    "    \"171 Å  (design)\",\n",
+    "    f\"174 Å  ($f$ = {f0 * wl_design / 174:.0f} mm,  $r_\\\\mathrm{{blur}}$ = {70 * abs(1 - 174/wl_design):.2f} mm)\",\n",
+    "]\n",
+    "\n",
+    "# Compute max disk radius for axis limits\n",
+    "r_max = max(70 * abs(1 - wl / wl_design) for wl in wavelengths_plot)\n",
+    "\n",
+    "with astropy.visualization.quantity_support():\n",
+    "    fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 6))\n",
+    "    for wl_val, color, label in zip(wavelengths_plot, colors, labels):\n",
+    "        f_wl  = f0 * wl_design / wl_val     # focal length at this wavelength, mm\n",
+    "        scale = 1.0 - f0 / f_wl             # defocus magnification (0 at design λ)\n",
+    "        x_det = px_2d[inside] * scale\n",
+    "        y_det = py_2d[inside] * scale\n",
+    "        ax.scatter(x_det, y_det, s=2, alpha=0.5, color=color,\n",
+    "                   label=label, rasterized=True)\n",
+    "    margin = 1.15\n",
+    "    ax.set_xlim(-r_max * margin, r_max * margin)\n",
+    "    ax.set_ylim(-r_max * margin, r_max * margin)\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.legend(markerscale=4, fontsize=9)\n",
+    "    ax.set_xlabel(\"$x$ (mm)\")\n",
+    "    ax.set_ylabel(\"$y$ (mm)\")\n",
+    "    ax.set_title(\"On-axis spot at 171 Å focal plane: chromatic aberration\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}