diff --git a/.claude/skills/qprogram-docs/SKILL.md b/.claude/skills/qprogram-docs/SKILL.md
index 580b76c..dccb692 100644
--- a/.claude/skills/qprogram-docs/SKILL.md
+++ b/.claude/skills/qprogram-docs/SKILL.md
@@ -37,6 +37,17 @@ A change to the DSL is not finished when the code is. These move together:
    mention in `docs/reference/qp-format.md`.
 5. The docstrings of anything that appears in `docs/reference/api-qprogram.md`,
    since that page is generated from them.
+6. The figures, when a change touches a program an example page plots. They are
+   built by running those programs, so a stale one is a picture of code that no
+   longer exists:
+
+   ```bash
+   uv run --extra viz python .claude/skills/qprogram-docs/scripts/build_example_plots.py
+   ```
+
+   Pass figure names to rebuild only some. Each figure is written twice, once
+   per site theme, into `docs/assets/plots/`, and pages embed the pair with the
+   `#only-light` / `#only-dark` fragments.
 
 If you are asked to do only part of this, do that part and say plainly which
 steps are still outstanding.
diff --git a/.claude/skills/qprogram-docs/scripts/build_example_plots.py b/.claude/skills/qprogram-docs/scripts/build_example_plots.py
new file mode 100644
index 0000000..e795eef
--- /dev/null
+++ b/.claude/skills/qprogram-docs/scripts/build_example_plots.py
@@ -0,0 +1,761 @@
+# Copyright 2026 Qilimanjaro Quantum Tech
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Render the figures the example pages embed.
+
+Every figure comes from running that page's own program on the reference platform, so a plot can
+never drift from the code printed above it. Change a program on a page, re-run this, and the picture
+follows.
+
+Each figure is written twice, once per documentation theme, as ``<name>-light.png`` and
+``<name>-dark.png``. The pages embed them with the ``#only-light`` / ``#only-dark`` fragments the
+site's stylesheet keys on ``data-md-color-scheme``, so the reader gets the one built for the surface
+they are looking at rather than a light figure dimmed by CSS.
+
+Usage::
+
+    uv run --extra viz python .claude/skills/qprogram-docs/scripts/build_example_plots.py
+    uv run --extra viz python .claude/skills/qprogram-docs/scripts/build_example_plots.py rabi cpmg
+
+The colours are the data-visualisation palette's categorical slots, in slot order, with a separate
+set of steps chosen for the dark surface rather than the light ones re-used. The sequential ramp is
+one hue and runs away from the surface in both themes: light to dark on white, dark to light on
+slate, so that "more" is always the mark that stands furthest off the page.
+"""
+
+from __future__ import annotations
+
+import sys
+import time
+import warnings
+from pathlib import Path
+from typing import TYPE_CHECKING, Any
+
+import matplotlib as mpl
+
+mpl.use("Agg")
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import LinearSegmentedColormap
+
+import qprogram as qp
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+    from matplotlib.axes import Axes
+    from matplotlib.figure import Figure
+
+OUT = Path(__file__).resolve().parents[4] / "docs" / "assets" / "plots"
+
+# --------------------------------------------------------------------------------------------
+# Themes
+# --------------------------------------------------------------------------------------------
+#
+# The surfaces are the site's own: Material's ``default`` scheme is white, and ``slate`` is
+# hsl(225, 15%, 14%). Matching them exactly is what lets a figure sit on the page without a seam.
+# The series steps are the categorical slots 1-4; the dark column is the same four hues stepped for
+# the dark surface, validated against it rather than flipped.
+
+THEMES: dict[str, dict[str, Any]] = {
+    "light": {
+        "surface": "#ffffff",
+        "text": "#0b0b0b",
+        "muted": "#52514e",
+        "grid": "#e6e6e3",
+        "series": ("#2a78d6", "#eb6834", "#1baf7a", "#eda100"),
+        "ramp": ("#cde2fb", "#9ec5f4", "#5598e7", "#2a78d6", "#256abf", "#184f95", "#0d366b"),
+    },
+    "dark": {
+        "surface": "#1e2129",
+        "text": "#e2e4e9",
+        "muted": "#a2a7b3",
+        "grid": "#33373f",
+        "series": ("#3987e5", "#d95926", "#199e70", "#c98500"),
+        "ramp": ("#232733", "#1d3a63", "#1c5497", "#2a78d6", "#5598e7", "#9ec5f4", "#cde2fb"),
+    },
+}
+
+DPI = 160
+WIDE = (7.2, 4.0)
+
+
+def _style(ax: Axes, theme: dict[str, Any]) -> None:
+    """Push the frame back so the data reads first."""
+    ax.set_facecolor(theme["surface"])
+    ax.grid(visible=True, color=theme["grid"], linewidth=0.8, zorder=0)
+    ax.set_axisbelow(True)
+    for side in ("top", "right"):
+        ax.spines[side].set_visible(False)
+    for side in ("left", "bottom"):
+        ax.spines[side].set_color(theme["grid"])
+    ax.tick_params(colors=theme["muted"], labelsize=9, length=0)
+    ax.xaxis.label.set_color(theme["muted"])
+    ax.yaxis.label.set_color(theme["muted"])
+    ax.xaxis.label.set_fontsize(10)
+    ax.yaxis.label.set_fontsize(10)
+
+
+def _legend(ax: Axes, theme: dict[str, Any], **kwargs: Any) -> None:
+    """A legend with no box, so it does not compete with the marks."""
+    leg = ax.legend(frameon=False, fontsize=9, **kwargs)
+    for text in leg.get_texts():
+        text.set_color(theme["text"])
+
+
+def _save(fig: Figure, name: str, mode: str, theme: dict[str, Any]) -> None:
+    OUT.mkdir(parents=True, exist_ok=True)
+    path = OUT / f"{name}-{mode}.png"
+    fig.savefig(path, dpi=DPI, facecolor=theme["surface"], bbox_inches="tight", pad_inches=0.25)
+    plt.close(fig)
+    print(f"  wrote {path.relative_to(OUT.parents[2])}")
+
+
+FIGURES: dict[str, Callable[[], Callable[[str, dict[str, Any]], None]]] = {}
+
+
+def figure(name: str) -> Callable[[Callable[[], Any]], Callable[[], Any]]:
+    """Register a figure builder.
+
+    The decorated function runs the page's program once and returns a ``draw(mode, theme)``
+    closure, so the data is computed a single time and rendered once per theme.
+
+    Args:
+        name (str): The figure's name, used for the output files and on the command line.
+
+    Returns:
+        The decorator that registers the builder under ``name``.
+    """
+
+    def register(fn: Callable[[], Any]) -> Callable[[], Any]:
+        FIGURES[name] = fn
+        return fn
+
+    return register
+
+
+# --------------------------------------------------------------------------------------------
+# Shared calibration
+# --------------------------------------------------------------------------------------------
+
+SCHEMA = qp.BusSchema.transmon()
+Q = SCHEMA.q
+
+
+def _iq(amplitude: float, duration: int) -> qp.waveforms.IQPair:
+    return qp.waveforms.IQPair(qp.waveforms.Square(amplitude, duration), qp.waveforms.Square(0.0, duration))
+
+
+READOUT = _iq(1.0, 2000)
+WEIGHTS = qp.waveforms.IQPair(qp.waveforms.Square(1.0, 2000), qp.waveforms.Square(1.0, 2000))
+
+
+# --------------------------------------------------------------------------------------------
+# Rabi oscillation
+# --------------------------------------------------------------------------------------------
+
+
+@figure("rabi")
+def _rabi() -> Callable[[str, dict[str, Any]], None]:
+    program = qp.QProgram(label="rabi", schema=SCHEMA)
+    gain = program.variable("gain", label="Drive amplitude", units="V")
+    with program.average(shots=1000), program.sweep(gain, qp.Range(start=0.0, stop=1.0, step=0.01)):
+        program.set_gain(Q[0].drive, gain)
+        program.play(Q[0].drive, "pi_pulse")
+        program.sync()
+        m0 = program.measure(Q[0].readout, "readout", "weights")
+
+    library = {
+        "pi_pulse": qp.waveforms.IQDrag(amplitude=0.5, duration=40, sigma=8, beta=0.1),
+        "readout": READOUT,
+        "weights": WEIGHTS,
+    }
+    model = qp.MockMeasurementModel(
+        response=lambda bus, env: np.sin(np.pi * env["gain"]) ** 2 + 0j,  # ruff: ignore[unused-lambda-argument]
+        noise=0.02,
+    )
+    data = qp.simulate(program.with_waveforms(library), model=model).get(m0)
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        x = data.coords["gain"]
+        ax.plot(x, data.sel(IQ="I"), color=theme["series"][0], linewidth=1.8, label="I", zorder=3)
+        ax.plot(x, data.sel(IQ="Q"), color=theme["series"][1], linewidth=1.8, label="Q", zorder=3)
+        ax.set_xlabel("Drive amplitude (V)")
+        ax.set_ylabel("Readout response")
+        _legend(ax, theme, loc="upper left")
+        _save(fig, "rabi", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# Qubit spectroscopy
+# --------------------------------------------------------------------------------------------
+
+
+@figure("qubit-spectroscopy")
+def _qubit_spectroscopy() -> Callable[[str, dict[str, Any]], None]:
+    program = qp.QProgram(label="qubit_spectroscopy", schema=SCHEMA)
+    freq = program.variable("freq", label="Drive frequency", units="Hz")
+    with program.average(shots=1000), program.sweep(freq, qp.Linspace(4.6e9, 5.4e9, num=201)):
+        program.set_frequency(Q[0].drive, freq)
+        program.play(Q[0].drive, "saturation")
+        program.sync()
+        m0 = program.measure(Q[0].readout, "readout", "weights")
+
+    library = {"saturation": _iq(0.02, 20000), "readout": READOUT, "weights": WEIGHTS}
+
+    def lorentzian(bus: str, env: dict[str, float]) -> complex:  # ruff: ignore[unused-function-argument]
+        f0, hwhm = 5.0e9, 8e6
+        return 1.0 / (1.0 + ((env["freq"] - f0) / hwhm) ** 2) + 0j
+
+    model = qp.MockMeasurementModel(response=lorentzian, noise=0.01)
+    data = qp.simulate(program.with_waveforms(library), model=model).get(m0)
+    magnitude = np.hypot(data.sel(IQ="I"), data.sel(IQ="Q"))
+    ghz = data.coords["freq"] / 1e9
+    peak = float(ghz[int(np.argmax(magnitude.values))])
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        ax.plot(ghz, magnitude, color=theme["series"][0], linewidth=1.8, zorder=3)
+        ax.axvline(peak, color=theme["muted"], linewidth=1.0, linestyle=(0, (4, 3)), zorder=2)
+        ax.annotate(
+            f"f01 = {peak:.3f} GHz",
+            xy=(peak, float(magnitude.max())),
+            xytext=(8, -6),
+            textcoords="offset points",
+            color=theme["text"],
+            fontsize=9,
+        )
+        ax.set_xlabel("Drive frequency (GHz)")
+        ax.set_ylabel("Readout magnitude")
+        _save(fig, "qubit-spectroscopy", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# T1 and Ramsey
+# --------------------------------------------------------------------------------------------
+
+_COHERENCE_LIBRARY = {
+    "pi_pulse": qp.waveforms.IQDrag(amplitude=0.5, duration=40, sigma=8, beta=0.1),
+    "pi_half": qp.waveforms.IQDrag(amplitude=0.25, duration=40, sigma=8, beta=0.1),
+    "readout": READOUT,
+    "weights": WEIGHTS,
+}
+
+
+@figure("t1")
+def _t1() -> Callable[[str, dict[str, Any]], None]:
+    t1 = qp.QProgram(label="t1", schema=SCHEMA)
+    delay = t1.variable("delay", label="Delay", units="ns")
+    with t1.average(shots=1000), t1.sweep(delay, qp.Linspace(0.0, 40_000.0, num=41)):
+        t1.play(Q[0].drive, "pi_pulse")
+        t1.wait(Q[0].drive, delay)
+        t1.sync()
+        m0 = t1.measure(
+            Q[0].readout,
+            "readout",
+            "weights",
+            fields=(qp.MeasurementField.IQ, qp.MeasurementField.STATE),
+        )
+
+    model = qp.MockMeasurementModel(
+        p_excited=lambda bus, env: float(np.exp(-env["delay"] / 12_000.0)),  # ruff: ignore[unused-lambda-argument]
+        seed=3,
+    )
+    population = qp.simulate(t1.with_waveforms(_COHERENCE_LIBRARY), model=model).get(
+        m0, field=qp.MeasurementField.STATE
+    )
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        micro = population.coords["delay"] / 1000.0
+        ax.plot(
+            micro,
+            population,
+            color=theme["series"][0],
+            linewidth=1.8,
+            marker="o",
+            markersize=3.5,
+            zorder=3,
+        )
+        ax.set_xlabel("Delay (μs)")
+        ax.set_ylabel("Excited-state population")
+        ax.set_ylim(-0.05, 1.05)
+        _save(fig, "t1", mode, theme)
+
+    return draw
+
+
+@figure("ramsey")
+def _ramsey() -> Callable[[str, dict[str, Any]], None]:
+    ramsey = qp.QProgram(label="ramsey", schema=SCHEMA)
+    delay = ramsey.variable("delay", label="Free evolution", units="ns")
+    with ramsey.average(shots=1000), ramsey.sweep(delay, qp.Linspace(0.0, 3000.0, num=151)):
+        ramsey.reset_phase(Q[0].drive)
+        ramsey.play(Q[0].drive, "pi_half")
+        ramsey.wait(Q[0].drive, delay)
+        ramsey.set_phase(Q[0].drive, 2 * np.pi * 2e-3 * delay)
+        ramsey.play(Q[0].drive, "pi_half")
+        ramsey.sync()
+        m1 = ramsey.measure(
+            Q[0].readout,
+            "readout",
+            "weights",
+            fields=(qp.MeasurementField.IQ, qp.MeasurementField.STATE),
+        )
+
+    def fringe(bus: str, env: dict[str, float]) -> float:  # ruff: ignore[unused-function-argument]
+        t = env["delay"]
+        return 0.5 * (1 - np.cos(2 * np.pi * 2e-3 * t) * np.exp(-t / 1500.0))
+
+    model = qp.MockMeasurementModel(p_excited=fringe, seed=1)
+    population = qp.simulate(ramsey.with_waveforms(_COHERENCE_LIBRARY), model=model).get(
+        m1, field=qp.MeasurementField.STATE
+    )
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        t = population.coords["delay"]
+        ax.plot(t, population, color=theme["series"][0], linewidth=1.8, zorder=3)
+        envelope = 0.5 * (1 + np.exp(-t / 1500.0))
+        ax.plot(t, envelope, color=theme["muted"], linewidth=1.0, linestyle=(0, (4, 3)), zorder=2)
+        ax.annotate(
+            "decay envelope",
+            xy=(float(t[110]), float(envelope[110])),
+            xytext=(0, 10),
+            textcoords="offset points",
+            color=theme["muted"],
+            fontsize=9,
+        )
+        ax.set_xlabel("Free evolution (ns)")
+        ax.set_ylabel("Excited-state population")
+        ax.set_ylim(-0.05, 1.05)
+        _save(fig, "ramsey", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# CZ chevron
+# --------------------------------------------------------------------------------------------
+
+
+@figure("cz-chevron")
+def _cz_chevron() -> Callable[[str, dict[str, Any]], None]:
+    schema = qp.BusSchema.flux_tunable_transmon()
+    q = schema.q
+    program = qp.QProgram(label="cz_chevron", schema=schema)
+    amp = program.variable("amp", label="Flux amplitude", units="V")
+    dur = program.variable("dur", label="Flux duration", units="ns")
+
+    with program.average(shots=1000), program.sweep(amp, qp.Range(0.0, 1.0, 0.01)):  # ruff: ignore[multiple-with-statements]
+        with program.sweep(dur, qp.Range(10, 210, 2)):
+            program.play(q[1].drive, "pi")
+            program.sync()
+            program.play(
+                q[0].flux,
+                qp.waveforms.FlatTop(amplitude=amp, duration=dur, smooth_duration=5),
+            )
+            program.sync()
+            m0 = program.measure(q[0].readout, "readout", "weights")
+            program.measure(q[1].readout, "readout", "weights")
+
+    library = {
+        "pi": qp.waveforms.IQDrag(amplitude=0.5, duration=40, sigma=8, beta=0.1),
+        "readout": READOUT,
+        "weights": WEIGHTS,
+    }
+
+    def chevron(bus: str, env: dict[str, float]) -> complex:  # ruff: ignore[unused-function-argument]
+        coupling = 0.01
+        detuning = 0.05 * (env["amp"] - 0.5)
+        rate = np.hypot(coupling, detuning)
+        return (coupling / rate) ** 2 * np.sin(np.pi * rate * env["dur"]) ** 2 + 0j
+
+    # The documented size: a 101 x 101 grid at 1000 shots with two measurements is 20 million model
+    # samples and around eight minutes on one core. It is run as written rather than coarsened,
+    # which is the point of building the figure from the page's own program.
+    result = qp.simulate(program.with_waveforms(library), model=qp.MockMeasurementModel(response=chevron))
+    data0 = result.get(m0)
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=(7.2, 4.4))
+        _style(ax, theme)
+        ax.grid(visible=False)
+        cmap = LinearSegmentedColormap.from_list("qp-sequential", theme["ramp"])
+        mesh = ax.pcolormesh(
+            data0.coords["dur"],
+            data0.coords["amp"],
+            data0.sel(IQ="I"),
+            cmap=cmap,
+            shading="nearest",
+            rasterized=True,
+        )
+        bar = fig.colorbar(mesh, ax=ax, pad=0.02)
+        bar.set_label("Population transferred", color=theme["muted"], fontsize=10)
+        bar.ax.tick_params(colors=theme["muted"], labelsize=9, length=0)
+        bar.outline.set_visible(True)
+        bar.outline.set_edgecolor(theme["grid"])
+        bar.outline.set_linewidth(0.8)
+        ax.set_xlabel("Flux duration (ns)")
+        ax.set_ylabel("Flux amplitude (V)")
+        _save(fig, "cz-chevron", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# Active reset
+# --------------------------------------------------------------------------------------------
+
+
+@figure("active-reset")
+def _active_reset() -> Callable[[str, dict[str, Any]], None]:
+    program = qp.QProgram(label="active_reset", schema=SCHEMA)
+    amp = program.variable("amp", label="Drive amplitude", units="V")
+    with program.average(shots=1000), program.sweep(amp, qp.Linspace(0.0, 1.0, num=21)):
+        check = program.measure(Q[0].readout, "readout", "weights", fields=(qp.MeasurementField.STATE,))
+        with program.if_(check.state == 1):
+            program.play(Q[0].drive, "pi")
+        with program.else_():
+            program.wait(Q[0].drive, 40)
+        program.sync()
+        program.set_gain(Q[0].drive, amp)
+        program.play(Q[0].drive, "pi")
+        program.sync()
+        m0 = program.measure(
+            Q[0].readout,
+            "readout",
+            "weights",
+            fields=(qp.MeasurementField.IQ, qp.MeasurementField.STATE),
+        )
+
+    library = {
+        "pi": qp.waveforms.IQDrag(amplitude=0.5, duration=40, sigma=8, beta=0.1),
+        "readout": READOUT,
+        "weights": WEIGHTS,
+    }
+    model = qp.MockMeasurementModel(
+        response=lambda bus, env: np.sin(np.pi * env["amp"]) ** 2 + 0j,  # ruff: ignore[unused-lambda-argument]
+        p_excited=lambda bus, env: 0.1,  # ruff: ignore[unused-lambda-argument]
+    )
+    result = qp.simulate(program.with_waveforms(library), model=model)
+    rabi = result.get(m0)
+    heralds = result.get(check, field=qp.MeasurementField.STATE)
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, (top, bottom) = plt.subplots(2, 1, figsize=(7.2, 4.8), sharex=True, gridspec_kw={"height_ratios": [2, 1]})
+        for ax in (top, bottom):
+            _style(ax, theme)
+        x = rabi.coords["amp"]
+        top.plot(x, rabi.sel(IQ="I"), color=theme["series"][0], linewidth=1.8, marker="o", markersize=3.5, zorder=3)
+        top.set_ylabel("Readout response")
+        bottom.plot(x, heralds, color=theme["series"][1], linewidth=1.8, marker="o", markersize=3.5, zorder=3)
+        bottom.set_ylim(0.0, 0.2)
+        bottom.set_ylabel("Herald rate")
+        bottom.set_xlabel("Drive amplitude (V)")
+        for ax, title in ((top, "Rabi sweep behind the reset"), (bottom, "Shots that needed the pi pulse")):
+            ax.set_title(title, color=theme["text"], fontsize=10, loc="left", pad=6)
+        fig.align_ylabels((top, bottom))
+        _save(fig, "active-reset", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# Resonator spectroscopy
+# --------------------------------------------------------------------------------------------
+
+
+@figure("resonator-spectroscopy")
+def _resonator_spectroscopy() -> Callable[[str, dict[str, Any]], None]:
+    program = qp.QProgram(label="resonator_spectroscopy", schema=SCHEMA)
+    lo = program.variable("lo", label="Readout LO", units="Hz")
+    with program.average(shots=1000), program.sweep(lo, qp.Linspace(7.0e9, 7.4e9, num=101)):
+        program.set_parameter(Q[0].readout, "lo_frequency", lo)
+        m0 = program.measure(Q[0].readout, "readout", "weights")
+
+    def notch(bus: str, env: dict[str, float]) -> complex:  # ruff: ignore[unused-function-argument]
+        f = env["q0/readout.lo_frequency"]
+        return 1.0 - 1.0 / (1.0 + ((f - 7.2e9) / 4e6) ** 2) + 0j
+
+    platform = qp.ReferencePlatform(schema=SCHEMA, model=qp.MockMeasurementModel(response=notch, noise=0.005))
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", qp.ExecutionWarning)
+        data = platform.execute(program.with_waveforms({"readout": READOUT, "weights": WEIGHTS})).get(m0)
+
+    ghz = data.coords["lo"] / 1e9
+    magnitude = np.hypot(data.sel(IQ="I"), data.sel(IQ="Q"))
+    dip = float(ghz[int(np.argmin(magnitude.values))])
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        ax.plot(ghz, magnitude, color=theme["series"][0], linewidth=1.8, zorder=3)
+        ax.axvline(dip, color=theme["muted"], linewidth=1.0, linestyle=(0, (4, 3)), zorder=2)
+        ax.annotate(
+            f"{dip:.3f} GHz",
+            xy=(dip, float(magnitude.min())),
+            xytext=(8, 14),
+            textcoords="offset points",
+            color=theme["text"],
+            fontsize=9,
+        )
+        ax.set_xlabel("Readout LO frequency (GHz)")
+        ax.set_ylabel("Transmitted magnitude")
+        _save(fig, "resonator-spectroscopy", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# CPMG on two qubits
+# --------------------------------------------------------------------------------------------
+
+
+@figure("cpmg")
+def _cpmg() -> Callable[[str, dict[str, Any]], None]:
+    @qp.fragment
+    def cpmg(f, drive, readout, tau):  # ruff: ignore[missing-type-function-argument, missing-return-type-private-function]
+        """Two pi/2 pulses around a train of four refocusing pi pulses."""
+        f.play(drive, "pi_half")
+        for _ in range(4):
+            f.wait(drive, tau / 2)
+            f.play(drive, "pi")
+            f.wait(drive, tau / 2)
+        f.play(drive, "pi_half")
+        f.sync([drive, readout])
+        f.measure(readout, "readout", "weights", fields=(qp.MeasurementField.STATE,))
+
+    program = qp.QProgram(label="cpmg", schema=SCHEMA)
+    tau = program.variable("tau", label="Pulse spacing", units="ns")
+    with program.average(shots=1000), program.sweep(tau, qp.Linspace(20.0, 2000.0, num=100)):
+        program.call(cpmg, Q[0].drive, Q[0].readout, tau)
+        program.call(cpmg, Q[1].drive, Q[1].readout, tau)
+
+    library = {
+        "pi_half": qp.waveforms.IQDrag(amplitude=0.25, duration=40, sigma=8, beta=0.1),
+        "pi": qp.waveforms.IQDrag(amplitude=0.5, duration=40, sigma=8, beta=0.1),
+        "readout": READOUT,
+        "weights": WEIGHTS,
+    }
+
+    def coherence(bus: str, env: dict[str, float]) -> float:
+        t2 = 3000.0 if bus.startswith("q0") else 1200.0
+        return 0.5 * (1.0 - np.exp(-(4 * env["tau"]) / t2))
+
+    result = qp.simulate(program.with_waveforms(library), model=qp.MockMeasurementModel(p_excited=coherence, seed=0))
+    q0 = result.get("m0", field=qp.MeasurementField.STATE)
+    q1 = result.get("m0_2", field=qp.MeasurementField.STATE)
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        x = q0.coords["tau"]
+        for series, data, label in ((0, q0, "q0"), (1, q1, "q1")):
+            ax.plot(x, data, color=theme["series"][series], linewidth=1.8, label=label, zorder=3)
+        ax.set_xlabel("Pulse spacing (ns)")
+        ax.set_ylabel("Excited-state population")
+        # The two curves converge at the long-spacing end, so a direct label there would sit on
+        # top of its neighbour; with two series the legend carries identity on its own.
+        _legend(ax, theme, loc="lower right")
+        _save(fig, "cpmg", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# Multiplexed readout
+# --------------------------------------------------------------------------------------------
+
+
+@figure("multiplexed-readout")
+def _multiplexed_readout() -> Callable[[str, dict[str, Any]], None]:
+    qubits = (0, 1, 2, 3)
+    program = qp.QProgram(label="multiplexed_rabi", schema=SCHEMA)
+    amp = program.variable("amp", label="Drive amplitude", units="V")
+    with program.average(shots=1000), program.sweep(amp, qp.Linspace(0.0, 1.0, num=21)):
+        for i in qubits:
+            program.set_gain(Q[i].drive, amp)
+            program.play(Q[i].drive, "pi_pulse")
+        program.sync()
+        for i in qubits:
+            program.measure(
+                Q[i].readout,
+                "readout",
+                "weights",
+                fields=(qp.MeasurementField.IQ, qp.MeasurementField.STATE),
+            )
+
+    library = qp.WaveformLibrary()
+    library.set("readout", _iq(0.9, 1000), element="q", idx=0, kind="readout")
+    library.set("readout", _iq(0.7, 3000), element="q", idx=2, kind="readout")
+    library.set("readout", _iq(0.5, 2000), element="q", kind="readout")
+    library.set("weights", WEIGHTS)
+    library.set("pi_pulse", qp.waveforms.IQDrag(0.5, 40, 8, 0.1))
+
+    period = {"q0/readout": 1.0, "q1/readout": 2.0, "q2/readout": 3.0, "q3/readout": 4.0}
+
+    def rabi(bus: str, env: dict[str, float]) -> float:
+        return float(np.sin(np.pi * env["amp"] / period[bus]) ** 2)
+
+    model = qp.MockMeasurementModel(response=lambda bus, env: rabi(bus, env) + 0j, p_excited=rabi, noise=0.02, seed=0)
+    result = qp.simulate(program.with_waveforms(library), model=model)
+    series = [result.get(i, field=qp.MeasurementField.STATE) for i in qubits]
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=WIDE)
+        _style(ax, theme)
+        x = series[0].coords["amp"]
+        for i, data in enumerate(series):
+            label = f"q{i}"
+            ax.plot(x, data, color=theme["series"][i], linewidth=1.8, label=label, zorder=3)
+            ax.annotate(
+                label,
+                xy=(float(x[-1]), float(data[-1])),
+                xytext=(6, 0),
+                textcoords="offset points",
+                color=theme["text"],
+                fontsize=9,
+                va="center",
+            )
+        ax.set_xlabel("Drive amplitude (V)")
+        ax.set_ylabel("Excited-state population")
+        ax.set_xlim(float(x.min()), float(x.max()) * 1.07)
+        # q0 peaks at the top of the axes, so the legend goes above them rather than inside.
+        _legend(ax, theme, loc="lower left", bbox_to_anchor=(0.0, 1.01), ncol=4)
+        _save(fig, "multiplexed-readout", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+# Single-shot readout
+# --------------------------------------------------------------------------------------------
+
+
+@figure("single-shot-readout")
+def _single_shot_readout() -> Callable[[str, dict[str, Any]], None]:
+    program = qp.QProgram(label="single_shot_readout", schema=SCHEMA)
+    prepared = program.variable("prepared", label="Prepared state")
+    shot = program.variable("shot", label="Shot index")
+    with program.sweep(prepared, qp.Values([0, 1])), program.sweep(shot, qp.Range(0, 1999, 1)):
+        program.set_gain(Q[0].drive, prepared)
+        program.play(Q[0].drive, "pi_pulse")
+        program.sync()
+        shots = program.measure(
+            Q[0].readout,
+            "readout",
+            "weights",
+            fields=(qp.MeasurementField.IQ, qp.MeasurementField.STATE),
+            name="shots",
+        )
+
+    class ReadoutModel:
+        """Two gaussian blobs in the IQ plane, classified by a threshold on I."""
+
+        def __init__(self, separation: float = 4.0, sigma: float = 1.0, seed: int = 0) -> None:
+            self.separation = separation
+            self.sigma = sigma
+            self._rng = np.random.default_rng(seed)
+
+        def sample(self, bus: str, env: dict[str, float]) -> qp.MeasurementSample:  # ruff: ignore[unused-method-argument]
+            center = self.separation if env["prepared"] else 0.0
+            i = center + self._rng.normal(0.0, self.sigma)
+            qv = self._rng.normal(0.0, self.sigma)
+            return qp.MeasurementSample(i=i, q=qv, state=int(i > self.separation / 2))
+
+    library = {
+        "pi_pulse": qp.waveforms.IQDrag(0.5, 40, 8, 0.1),
+        "readout": READOUT,
+        "weights": WEIGHTS,
+    }
+    result = qp.simulate(program.with_waveforms(library), model=ReadoutModel(seed=0))
+    iq = result.get(shots)
+    state = result.get("shots", field=qp.MeasurementField.STATE).values
+    fidelity = 1 - (state[0].mean() + (1 - state[1].mean())) / 2
+
+    def draw(mode: str, theme: dict[str, Any]) -> None:
+        fig, ax = plt.subplots(figsize=(5.6, 5.0))
+        _style(ax, theme)
+        for prep, series, label in ((0, 0, "prepared |0>"), (1, 1, "prepared |1>")):
+            ax.scatter(
+                iq[prep].sel(IQ="I"),
+                iq[prep].sel(IQ="Q"),
+                s=5,
+                alpha=0.45,
+                linewidths=0,
+                color=theme["series"][series],
+                label=label,
+                zorder=3,
+            )
+        ax.axvline(2.0, color=theme["muted"], linewidth=1.0, linestyle=(0, (4, 3)), zorder=4)
+        ax.annotate(
+            "threshold",
+            xy=(2.0, 1.0),
+            xycoords=("data", "axes fraction"),
+            xytext=(6, -12),
+            textcoords="offset points",
+            color=theme["muted"],
+            fontsize=9,
+            va="top",
+        )
+        ax.set_xlabel("I (arb.)")
+        ax.set_ylabel("Q (arb.)")
+        ax.set_aspect("equal", adjustable="box")
+        ax.set_title(f"Assignment fidelity {fidelity:.3f}", color=theme["text"], fontsize=10, loc="left", pad=6)
+        _legend(ax, theme, loc="upper left", markerscale=2.4)
+        _save(fig, "single-shot-readout", mode, theme)
+
+    return draw
+
+
+# --------------------------------------------------------------------------------------------
+
+
+def main(argv: list[str]) -> int:
+    """Build the requested figures, or every one of them when given no names.
+
+    Args:
+        argv (list[str]): Figure names to build. Empty means all of them.
+
+    Returns:
+        A process exit status: ``0`` on success, ``1`` when a name is not a known figure.
+    """
+    wanted = argv or sorted(FIGURES)
+    unknown = [name for name in wanted if name not in FIGURES]
+    if unknown:
+        print(f"unknown figure(s): {', '.join(unknown)}", file=sys.stderr)
+        print(f"available: {', '.join(sorted(FIGURES))}", file=sys.stderr)
+        return 1
+    for name in wanted:
+        started = time.perf_counter()
+        print(f"{name}: running the program...")
+        draw = FIGURES[name]()
+        for mode, theme in THEMES.items():
+            draw(mode, theme)
+        print(f"  {time.perf_counter() - started:.1f}s")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv[1:]))
diff --git a/changelog/29.changed.md b/changelog/29.changed.md
new file mode 100644
index 0000000..fd2a043
--- /dev/null
+++ b/changelog/29.changed.md
@@ -0,0 +1 @@
+The example pages now carry the figures their programs produce. Each one is rendered by `.claude/skills/qprogram-docs/scripts/build_example_plots.py`, which builds the page's own program, runs it on the reference platform, and writes the result to `docs/assets/plots/`, so a plot cannot drift from the code printed above it. Every figure is written once per site theme and the page picks the one built for the surface being read on, so the dark variant is its own render on the dark surface rather than a light figure behind a filter.
diff --git a/docs/assets/plots/active-reset-dark.png b/docs/assets/plots/active-reset-dark.png
new file mode 100644
index 0000000..640bd39
Binary files /dev/null and b/docs/assets/plots/active-reset-dark.png differ
diff --git a/docs/assets/plots/active-reset-light.png b/docs/assets/plots/active-reset-light.png
new file mode 100644
index 0000000..8ce809e
Binary files /dev/null and b/docs/assets/plots/active-reset-light.png differ
diff --git a/docs/assets/plots/cpmg-dark.png b/docs/assets/plots/cpmg-dark.png
new file mode 100644
index 0000000..3f0360e
Binary files /dev/null and b/docs/assets/plots/cpmg-dark.png differ
diff --git a/docs/assets/plots/cpmg-light.png b/docs/assets/plots/cpmg-light.png
new file mode 100644
index 0000000..e6bf6ce
Binary files /dev/null and b/docs/assets/plots/cpmg-light.png differ
diff --git a/docs/assets/plots/cz-chevron-dark.png b/docs/assets/plots/cz-chevron-dark.png
new file mode 100644
index 0000000..ebb4c04
Binary files /dev/null and b/docs/assets/plots/cz-chevron-dark.png differ
diff --git a/docs/assets/plots/cz-chevron-light.png b/docs/assets/plots/cz-chevron-light.png
new file mode 100644
index 0000000..fb7770d
Binary files /dev/null and b/docs/assets/plots/cz-chevron-light.png differ
diff --git a/docs/assets/plots/multiplexed-readout-dark.png b/docs/assets/plots/multiplexed-readout-dark.png
new file mode 100644
index 0000000..de5c570
Binary files /dev/null and b/docs/assets/plots/multiplexed-readout-dark.png differ
diff --git a/docs/assets/plots/multiplexed-readout-light.png b/docs/assets/plots/multiplexed-readout-light.png
new file mode 100644
index 0000000..8039e13
Binary files /dev/null and b/docs/assets/plots/multiplexed-readout-light.png differ
diff --git a/docs/assets/plots/qubit-spectroscopy-dark.png b/docs/assets/plots/qubit-spectroscopy-dark.png
new file mode 100644
index 0000000..427ef78
Binary files /dev/null and b/docs/assets/plots/qubit-spectroscopy-dark.png differ
diff --git a/docs/assets/plots/qubit-spectroscopy-light.png b/docs/assets/plots/qubit-spectroscopy-light.png
new file mode 100644
index 0000000..0ce91fa
Binary files /dev/null and b/docs/assets/plots/qubit-spectroscopy-light.png differ
diff --git a/docs/assets/plots/rabi-dark.png b/docs/assets/plots/rabi-dark.png
new file mode 100644
index 0000000..4ef317a
Binary files /dev/null and b/docs/assets/plots/rabi-dark.png differ
diff --git a/docs/assets/plots/rabi-light.png b/docs/assets/plots/rabi-light.png
new file mode 100644
index 0000000..c62690a
Binary files /dev/null and b/docs/assets/plots/rabi-light.png differ
diff --git a/docs/assets/plots/ramsey-dark.png b/docs/assets/plots/ramsey-dark.png
new file mode 100644
index 0000000..7a502f3
Binary files /dev/null and b/docs/assets/plots/ramsey-dark.png differ
diff --git a/docs/assets/plots/ramsey-light.png b/docs/assets/plots/ramsey-light.png
new file mode 100644
index 0000000..ea9630c
Binary files /dev/null and b/docs/assets/plots/ramsey-light.png differ
diff --git a/docs/assets/plots/resonator-spectroscopy-dark.png b/docs/assets/plots/resonator-spectroscopy-dark.png
new file mode 100644
index 0000000..a8b0d1a
Binary files /dev/null and b/docs/assets/plots/resonator-spectroscopy-dark.png differ
diff --git a/docs/assets/plots/resonator-spectroscopy-light.png b/docs/assets/plots/resonator-spectroscopy-light.png
new file mode 100644
index 0000000..0e53175
Binary files /dev/null and b/docs/assets/plots/resonator-spectroscopy-light.png differ
diff --git a/docs/assets/plots/single-shot-readout-dark.png b/docs/assets/plots/single-shot-readout-dark.png
new file mode 100644
index 0000000..f19d12b
Binary files /dev/null and b/docs/assets/plots/single-shot-readout-dark.png differ
diff --git a/docs/assets/plots/single-shot-readout-light.png b/docs/assets/plots/single-shot-readout-light.png
new file mode 100644
index 0000000..f9ec631
Binary files /dev/null and b/docs/assets/plots/single-shot-readout-light.png differ
diff --git a/docs/assets/plots/t1-dark.png b/docs/assets/plots/t1-dark.png
new file mode 100644
index 0000000..7ac2f7b
Binary files /dev/null and b/docs/assets/plots/t1-dark.png differ
diff --git a/docs/assets/plots/t1-light.png b/docs/assets/plots/t1-light.png
new file mode 100644
index 0000000..f0296a3
Binary files /dev/null and b/docs/assets/plots/t1-light.png differ
diff --git a/docs/examples/active-reset.md b/docs/examples/active-reset.md
index 6259477..5a98198 100644
--- a/docs/examples/active-reset.md
+++ b/docs/examples/active-reset.md
@@ -213,7 +213,7 @@ library = {
 result = qp.simulate(
     program.with_waveforms(library),
     model=qp.MockMeasurementModel(
-        response=lambda bus, env: np.sin(np.pi * env["amp"] / 2) ** 2 + 0j,
+        response=lambda bus, env: np.sin(np.pi * env["amp"]) ** 2 + 0j,
         p_excited=lambda bus, env: 0.1,
     ),
 )
@@ -232,6 +232,13 @@ worth doing: a herald rate that climbs over a run is the readout heating the
 qubit or the previous shot's pulse leaking, and neither shows up in the Rabi
 curve until it has already distorted it.
 
+Both records plotted against the same sweep, the Rabi curve above and the
+herald rate below, is the shape worth watching: the reset should hold flat
+while the experiment underneath it moves.
+
+![Two stacked panels sharing a drive-amplitude axis: a Rabi curve peaking at 0.5 V and falling back to 0, and a herald rate scattered around 0.1.](../assets/plots/active-reset-light.png#only-light)
+![Two stacked panels sharing a drive-amplitude axis: a Rabi curve peaking at 0.5 V and falling back to 0, and a herald rate scattered around 0.1.](../assets/plots/active-reset-dark.png#only-dark)
+
 Asking `check` for a field it never requested raises rather than substituting
 one:
 
diff --git a/docs/examples/cpmg-fragments.md b/docs/examples/cpmg-fragments.md
index fc7ff1a..04e9e55 100644
--- a/docs/examples/cpmg-fragments.md
+++ b/docs/examples/cpmg-fragments.md
@@ -207,6 +207,11 @@ both receive the bus string, so one model can describe a whole chip. Without
 that branch the two qubits would return the same numbers and the second
 measurement would be teaching nothing.
 
+One fragment, two call sites, two curves that separate because the chip does:
+
+![Excited-state population against pulse spacing for two qubits, both rising toward 0.5, with q1 losing coherence faster than q0.](../assets/plots/cpmg-light.png#only-light)
+![Excited-state population against pulse spacing for two qubits, both rising toward 0.5, with q1 losing coherence faster than q0.](../assets/plots/cpmg-dark.png#only-dark)
+
 The total free evolution is four times the spacing, which is why the model
 reads `4 * env["tau"]` rather than `env["tau"]`. The reference executor times
 nothing, so the `wait` operations contribute no duration of their own and that
diff --git a/docs/examples/cz-chevron.md b/docs/examples/cz-chevron.md
index 3e173b4..3b2c082 100644
--- a/docs/examples/cz-chevron.md
+++ b/docs/examples/cz-chevron.md
@@ -199,6 +199,9 @@ plt.xlabel("Flux duration (ns)")
 plt.ylabel("Flux amplitude (V)")
 ```
 
+![Heatmap of transferred population against flux duration and amplitude, with interference fringes converging to a chevron tip at 0.5 V.](../assets/plots/cz-chevron-light.png#only-light)
+![Heatmap of transferred population against flux duration and amplitude, with interference fringes converging to a chevron tip at 0.5 V.](../assets/plots/cz-chevron-dark.png#only-dark)
+
 `pcolormesh` takes the x axis first, so the inner sweep goes first and the
 outer one second, the opposite of the dimension order in `data0.dims`.
 matplotlib is not a runtime dependency; it comes with the `viz` extra,
diff --git a/docs/examples/index.md b/docs/examples/index.md
index b5cf4e0..6dc093f 100644
--- a/docs/examples/index.md
+++ b/docs/examples/index.md
@@ -7,6 +7,13 @@ the reference platform, the pure-Python interpreter that ships with the
 package. The pages are ordered so that each one adds a few pieces to the ones
 before it, and each names what those are in its opening paragraphs.
 
+The figures come from running the programs. Each one is produced by
+`.claude/skills/qprogram-docs/scripts/build_example_plots.py`, which builds the
+page's own program, executes it on the reference platform, and writes the
+result to `docs/assets/plots/`, so a plot cannot drift from the code printed
+above it. Every figure is written twice, once per site theme, and the page
+picks the one built for the surface you are reading on.
+
 | Example | The features it exercises |
 |---|---|
 | [Rabi oscillation](rabi.md) | One sweep inside an averaging block, one pulse, one measurement, and a result with a single sweep dimension plus the trailing `IQ` axis. |
diff --git a/docs/examples/multiplexed-readout.md b/docs/examples/multiplexed-readout.md
index 1f65cc6..05213c1 100644
--- a/docs/examples/multiplexed-readout.md
+++ b/docs/examples/multiplexed-readout.md
@@ -257,6 +257,12 @@ model = qp.MockMeasurementModel(
 Without the branch on `bus` all four records come back identical, which looks
 like a working multiplexed readout and is not one.
 
+Four records off one shot, each qubit reaching its pi amplitude somewhere
+different, which is the calibration this measurement exists to produce:
+
+![Excited-state population against drive amplitude for four qubits, each with a different Rabi period, labelled q0 through q3.](../assets/plots/multiplexed-readout-light.png#only-light)
+![Excited-state population against drive amplitude for four qubits, each with a different Rabi period, labelled q0 through q3.](../assets/plots/multiplexed-readout-dark.png#only-dark)
+
 ## Adapting it
 
 To move a program to different qubits, `rebind` maps element indices and
diff --git a/docs/examples/qubit-spectroscopy.md b/docs/examples/qubit-spectroscopy.md
index 565c590..82882d3 100644
--- a/docs/examples/qubit-spectroscopy.md
+++ b/docs/examples/qubit-spectroscopy.md
@@ -166,13 +166,19 @@ plt.plot(data.coords["freq"] / 1e9, magnitude)
 plt.xlabel("Drive frequency (GHz)")
 plt.ylabel("Readout magnitude")
 
-f01 = float(data.coords["freq"][int(np.argmax(magnitude))])  # 5.0e9
+f01 = float(data.coords["freq"][int(np.argmax(magnitude.values))])  # 5.0e9
 ```
 
+![Readout magnitude against drive frequency, flat except for a sharp peak at 5.000 GHz marked as f01.](../assets/plots/qubit-spectroscopy-light.png#only-light)
+![Readout magnitude against drive frequency, flat except for a sharp peak at 5.000 GHz marked as f01.](../assets/plots/qubit-spectroscopy-dark.png#only-dark)
+
 `np.hypot` over two `sel` results returns a `DataArray` with dims `("freq",)`,
 so the coordinate survives the arithmetic and the peak can be read back as a
-frequency. matplotlib is not a runtime dependency; it comes with the `viz`
-extra, installed with `pip install "qprogram[viz]"`.
+frequency. `np.argmax` wants the underlying array rather than the `DataArray`,
+which is what `.values` is for; handing it the labelled array raises
+`ValueError: dimensions ('freq',) must have the same length as the number of
+data dimensions, ndim=0`. matplotlib is not a runtime dependency; it comes with
+the `viz` extra, installed with `pip install "qprogram[viz]"`.
 
 ## Adapting it
 
diff --git a/docs/examples/rabi.md b/docs/examples/rabi.md
index d815e83..f9451be 100644
--- a/docs/examples/rabi.md
+++ b/docs/examples/rabi.md
@@ -170,7 +170,7 @@ of exactly `0.0`. Shape the response to see a curve:
 import numpy as np
 
 model = qp.MockMeasurementModel(
-    response=lambda bus, env: np.sin(np.pi * env["gain"] / 2) ** 2 + 0j,
+    response=lambda bus, env: np.sin(np.pi * env["gain"]) ** 2 + 0j,
     noise=0.02,
 )
 result = qp.simulate(resolved, model=model)
@@ -230,6 +230,9 @@ plt.xlabel("Drive amplitude (V)")
 plt.legend()
 ```
 
+![Readout response against drive amplitude. I rises to a maximum of 1 at 0.5 V and falls back to 0 by 1.0 V, while Q stays flat at 0.](../assets/plots/rabi-light.png#only-light)
+![Readout response against drive amplitude. I rises to a maximum of 1 at 0.5 V and falls back to 0 by 1.0 V, while Q stays flat at 0.](../assets/plots/rabi-dark.png#only-dark)
+
 The axis label is written out here rather than read from the variable. The
 `label` and `units` given to `program.variable` travel with the program into
 `.qp` and are there for tooling and for whoever reads the file; the executor
diff --git a/docs/examples/resonator-spectroscopy.md b/docs/examples/resonator-spectroscopy.md
index 31123b3..94de032 100644
--- a/docs/examples/resonator-spectroscopy.md
+++ b/docs/examples/resonator-spectroscopy.md
@@ -206,6 +206,12 @@ result.get(m0).shape  # (101, 2)
 platform.parameters  # {"q0/readout.lo_frequency": 7400000000.0}
 ```
 
+![Transmitted magnitude against readout LO frequency, flat near 1 except for a sharp dip to 0 at 7.200 GHz.](../assets/plots/resonator-spectroscopy-light.png#only-light)
+![Transmitted magnitude against readout LO frequency, flat near 1 except for a sharp dip to 0 at 7.200 GHz.](../assets/plots/resonator-spectroscopy-dark.png#only-dark)
+
+The dip is the resonator, and its centre is the frequency the readout pulse
+wants to be at.
+
 A swept parameter reaches the model differently from a swept variable. The
 earlier pages read `env["freq"]` or `env["delay"]`, the id of the loop
 variable; a parameter write puts its value in the same `env` under
diff --git a/docs/examples/single-shot-readout.md b/docs/examples/single-shot-readout.md
index 2e5ed97..19db9f2 100644
--- a/docs/examples/single-shot-readout.md
+++ b/docs/examples/single-shot-readout.md
@@ -198,6 +198,9 @@ plt.axvline(2.0, color="k", lw=0.5)  # the classifier threshold
 plt.legend()
 ```
 
+![Scatter of single shots in the IQ plane: two well-separated gaussian blobs for the ground and excited preparations, split by a threshold at I = 2.](../assets/plots/single-shot-readout-light.png#only-light)
+![Scatter of single shots in the IQ plane: two well-separated gaussian blobs for the ground and excited preparations, split by a threshold at I = 2.](../assets/plots/single-shot-readout-dark.png#only-dark)
+
 Four thousand shots run in well under a tenth of a second, so this is the
 cheapest program in the section despite having the most records. The
 combination to be careful with is not the shot count but the shot count times
diff --git a/docs/examples/t1-and-ramsey.md b/docs/examples/t1-and-ramsey.md
index bc49b39..4eb5733 100644
--- a/docs/examples/t1-and-ramsey.md
+++ b/docs/examples/t1-and-ramsey.md
@@ -199,6 +199,12 @@ from it per shot. Both receive the same `(bus, env)` pair, so a decay written
 against `env["delay"]` is all it takes to give the curve a shape. Without a
 `p_excited` argument every shot classifies as 0.
 
+The curve is the exponential the model was given, sampled a thousand shots per
+point, and the scatter around it is the Bernoulli noise of that count:
+
+![Excited-state population against delay, decaying exponentially from 1 toward 0 over 40 microseconds.](../assets/plots/t1-light.png#only-light)
+![Excited-state population against delay, decaying exponentially from 1 toward 0 over 40 microseconds.](../assets/plots/t1-dark.png#only-dark)
+
 The `STATE` array has no trailing `"IQ"` dimension, because a classified
 outcome is one number per shot rather than a pair. `result.get(m0)` on the same
 handle still returns the IQ field with dims `("delay", "IQ")` and shape
@@ -232,6 +238,13 @@ The 20 ns spacing that `Linspace(0.0, 3000.0, num=151)` resolves to samples the
 count: the delay axis has to resolve the artificial detuning, not just reach
 far enough to see the envelope.
 
+The fringe and the envelope are the two things the measurement separates. The
+oscillation is the 2 MHz the program put there; the decay it sits inside is the
+one the qubit contributed:
+
+![Excited-state population against free evolution time, oscillating at 2 MHz inside a decaying envelope drawn as a dashed line.](../assets/plots/ramsey-light.png#only-light)
+![Excited-state population against free evolution time, oscillating at 2 MHz inside a decaying envelope drawn as a dashed line.](../assets/plots/ramsey-dark.png#only-dark)
+
 Nothing in either run knows about relaxation or precession. `wait` evaluates
 its expression and returns, `set_phase` and `reset_phase` do the same, and the
 curves come entirely from the two model callbacks reading `env["delay"]`. See
diff --git a/docs/getting-started.md b/docs/getting-started.md
index fb8052a..357e7ca 100644
--- a/docs/getting-started.md
+++ b/docs/getting-started.md
@@ -223,7 +223,7 @@ import numpy as np
 import qprogram as qp
 
 model = qp.MockMeasurementModel(
-    response=lambda bus, env: np.sin(np.pi * env["gain"] / 2) ** 2 + 0j,
+    response=lambda bus, env: np.sin(np.pi * env["gain"]) ** 2 + 0j,
     noise=0.02,
     seed=7,
 )
diff --git a/docs/guide/execution.md b/docs/guide/execution.md
index c6d03da..d29dda3 100644
--- a/docs/guide/execution.md
+++ b/docs/guide/execution.md
@@ -34,7 +34,7 @@ import numpy as np
 import qprogram as qp
 
 model = qp.MockMeasurementModel(
-    response=lambda bus, env: np.sin(np.pi * env["g"] / 2) ** 2 + 0j,
+    response=lambda bus, env: np.sin(np.pi * env["g"]) ** 2 + 0j,
     noise=0.02,
     seed=7,
 )