he11lib/docs/superpowers/specs/2026-07-02-gyrotron-mode-purity-design.md

# he11lib — Gyrotron Beam Mode Purity Reconstruction

**Date:** 2026-07-02
**Status:** Approved for planning

## Purpose

A general-purpose, reusable Python library for reconstructing the Laguerre-Gauss
(LG) modal content ("mode purity") of a free-space-propagating gyrotron RF beam,
from a set of thermal (flux) images captured at different distances from the
gyrotron output window. The library must account for real-world measurement
issues: unknown/partial camera geometry, sensor noise, target thermal-diffusion
blur, and unknown beam pointing/centering — not just an idealized intensity fit.

This is a from-scratch library (empty project directory), intended to be used by
others beyond the initial author.

## Scope

- Beam regime: **free-space quasi-optical propagation** (post mode-converter /
  output window), not waveguide-confined hybrid modes.
- Mode basis: **Laguerre-Gauss modes** `LG_{p,l}`, referenced to a **known**
  waist size `w0` and waist location `z0` supplied by the user (design values
  from the gyrotron/mode-converter specs) — these are inputs, not fit
  parameters.
- Input data: **pre-processed NumPy flux arrays** (already extracted from raw
  NVF radiometric camera files by the user's own pipeline). Dead-pixel
  correction, background/ambient subtraction, and saturation-clipping detection
  are already handled upstream and out of scope for this library. Because the
  data is flux (not temperature), there is no temperature-to-power nonlinearity
  to correct.
- Number of measurement planes: **3 to 10**, at distances spanning roughly
  0–100 cm from the output window (e.g. 30/40/50/60 cm), known to a nominal
  precision (set by a translation stage or tape measure).
- Camera geometry (pixel-to-length scale, viewing angle relative to beam axis):
  may be **known from a prior calibration** (supplied as input) **or unknown**
  (estimated jointly with everything else).
- Beam **transverse center and pointing angle** are always unknown and must be
  estimated from the images.
- Residual sensor noise (NETD-type, after upstream preprocessing) is
  **auto-estimated per image**; no user-supplied noise parameter required for
  the default path.
- Target thermal-diffusion blur correction (deconvolution) is **optional**,
  parametrized by target thermal diffusivity + exposure/dwell time.
- Deliverables include full API docs and an end-to-end example
  script/notebook demonstrating the complete pipeline.

## Architecture

Modular/composable design: each concern is an independent, testable component
with a clear interface, wired together by a high-level orchestrator. This
supports both simple end-to-end use and power-user access to individual
stages (e.g., using `LGBasis` standalone, or swapping in a custom noise
model).

### Components

- **`data.py` — `MeasurementPlane`, `ReconstructionResult`**
  `MeasurementPlane`: container for one measurement — flux array, nominal `z`
  distance, optional known pixel scale (mm/px) and viewing angle (deg),
  optional metadata (timestamp, label). `ReconstructionResult`: container for
  all pipeline outputs (see below).

- **`geometry.py` — `GeometryCalibration`**
  Applies projective/perspective correction to compensate for oblique camera
  viewing angle, and converts pixel coordinates to physical length units. If
  pixel scale and/or viewing angle are supplied on a `MeasurementPlane`, uses
  them directly; if not supplied, exposes them as free parameters to be
  solved jointly by the `ModalFitter`.

- **`noise.py` — `NoiseEstimator`**
  Estimates residual per-image noise standard deviation automatically (e.g.
  from low-signal/background regions or high-frequency residual content).
  Produces per-pixel or per-image weights used by the `ModalFitter`'s
  noise-weighted least squares.

- **`deconvolution.py` — `DiffusionDeconvolver`**
  Optional step. Given target material thermal diffusivity and
  exposure/dwell time, builds a blur kernel modeling lateral heat spreading
  in the absorbing target and deconvolves it from each plane before fitting.
  Disabled by default.

- **`modes.py` — `LGBasis`**
  Given reference `w0`, `z0`, generates complex LG mode fields `LG_{p,l}(x,
  y, z)` at arbitrary transverse coordinates and axial distance `z`,
  including correct Gouy phase and beam-radius evolution. Supports
  evaluating a finite candidate set of `(p, l)` indices and projecting an
  arbitrary complex field onto the basis (used both by the modal fit and by
  the phase-retrieval fallback).

- **`fitting.py` — `ModalFitter`**
  Core reconstruction path. Parameters: complex LG coefficients
  (amplitude + phase) for each candidate mode, beam transverse center `(x,
  y)` per plane, a shared beam pointing angle (tilt), and any uncalibrated
  geometry parameters (pixel scale / viewing angle) not supplied on the
  `MeasurementPlane`s. Objective: noise-weighted nonlinear least-squares
  residual between modeled `|Σ c_j · LG_j(x, y, z)|²` and measured flux,
  summed over all planes.

  **Automatic mode-set growth**: starts from `LG_00`, incrementally adds
  candidate modes in shells of increasing order, and stops growing once an
  information-criterion (e.g. BIC) / residual-reduction test shows no
  meaningful improvement, subject to a configurable maximum order cap (for
  tractability and as a safety bound). Warns (does not error) if the cap is
  hit while still improving, or if final residuals remain large.

- **`phase_retrieval.py` — `PhaseRetriever`**
  Optional/fallback path. Gerchberg-Saxton-style iterative multi-plane phase
  retrieval: propagates a trial complex field back and forth between
  measurement planes (via free-space propagation, e.g. angular spectrum),
  enforcing the measured amplitude (sqrt of flux) at each plane, without
  assuming a finite mode basis. Used when explicitly requested, or
  automatically as a fallback when `ModalFitter`'s residual stays high after
  mode-set growth completes. The recovered field is then projected onto
  `LGBasis` to produce a purity table.

- **`synthetic.py` — `SyntheticBeamGenerator`**
  Forward model: given known mode coefficients, `w0`/`z0`, geometry
  (center, tilt, viewing angle, pixel scale), noise level, and optional
  diffusion blur, generates synthetic multi-plane flux images. Used for
  library validation (recover known ground truth) and for users to test
  experimental design (e.g., "would these 4 distances separate my modes?").

- **`reconstruct.py` — `BeamReconstructor`**
  High-level orchestrator: given a list of `MeasurementPlane`s and
  configuration (known `w0`/`z0`, optional deconvolution params, optional
  mode-set cap, optional forced phase-retrieval mode), runs geometry
  correction → noise estimation → optional deconvolution → modal fit (with
  automatic growth) → optional phase-retrieval fallback → produces a
  `ReconstructionResult`. Each stage remains independently accessible for
  power users.

- **`plotting.py`**
  Diagnostic visualizations: measured vs. reconstructed intensity per plane,
  residual maps, mode purity bar chart, beam center/pointing trace across
  planes.

### Data flow

1. Build a list of `MeasurementPlane` objects (flux array + nominal z +
   optional known geometry per plane).
2. `GeometryCalibration` applies projective correction using supplied
   pixel-scale/viewing-angle, or flags them as unknowns.
3. `NoiseEstimator` computes per-plane noise weights.
4. `DiffusionDeconvolver` optionally deblurs each plane.
5. `ModalFitter` runs the joint noise-weighted nonlinear least-squares fit
   over LG coefficients + center + pointing + any uncalibrated geometry
   params, growing the mode set automatically.
6. If requested, or if fit residual remains high, `PhaseRetriever` runs as a
   fallback and its result is projected onto `LGBasis`.
7. `BeamReconstructor` assembles a `ReconstructionResult` containing: mode
   purity table (power fraction + phase per mode), reconstructed complex
   field, fitted beam center/pointing per plane, geometry parameters used or
   fitted, per-plane residual maps, and coefficient uncertainties (from the
   fit's covariance).

## Testing strategy

`SyntheticBeamGenerator` is the backbone of validation: generate synthetic
multi-plane data from known ground-truth mode content, geometry, and noise,
run it through the full pipeline (and through individual components in
isolation), and assert recovered parameters match ground truth within
tolerance. Individual components also get targeted unit tests against known
analytic cases (e.g. `LGBasis` orthogonality and known Gouy phase values,
single-pure-mode recovery, geometry correction on synthetic projective
distortions).

## Error handling

Validate only at boundaries: reject malformed `MeasurementPlane` inputs
(mismatched array shapes, non-positive `z`, fewer than 3 planes). Warn
(rather than raise) when automatic mode-set growth hits its configured cap
while still improving, or when final fit residuals remain large — the caller
gets the result plus a diagnostic flag, not a crash.

## Dependencies

NumPy, SciPy (`optimize` for the nonlinear least-squares fit, `ndimage` for
projective transforms and deconvolution), Matplotlib for diagnostic
plotting. No GPU requirement.

## Package layout

```
he11lib/
  data.py            # MeasurementPlane, ReconstructionResult
  geometry.py        # GeometryCalibration
  noise.py           # NoiseEstimator
  deconvolution.py   # DiffusionDeconvolver
  modes.py           # LGBasis
  fitting.py         # ModalFitter (+ automatic mode-set growth)
  phase_retrieval.py # PhaseRetriever
  synthetic.py       # SyntheticBeamGenerator
  reconstruct.py      # BeamReconstructor (orchestrator)
  plotting.py         # diagnostic visualizations
docs/
  ...                 # API docs
examples/
  full_pipeline_example.py  # end-to-end demo of the whole pipeline
tests/
  ...
```

## Deliverables

- Full implementation of all components above.
- API documentation for the public interface of every module.
- An end-to-end example (script or notebook) exercising the complete
  pipeline on synthetic data, from `SyntheticBeamGenerator` through
  `BeamReconstructor` to `plotting.py` diagnostics.
- Test suite covering synthetic ground-truth recovery and per-component unit
  tests.