he11lib/docs/api.md

# he11lib API Reference

`he11lib` reconstructs the Laguerre-Gauss (LG) modal content ("mode purity")
of a free-space-propagating gyrotron RF beam from a set of thermal (flux)
images taken at different distances from the output window.

See `examples/full_pipeline_example.py` for a runnable end-to-end
demonstration, and
`docs/superpowers/specs/2026-07-02-gyrotron-mode-purity-design.md` for the
full design rationale.

Every class/function below is exported from the top-level `he11lib` package
(e.g. `from he11lib import BeamReconstructor`), except where noted.

## Quick start

```python
from he11lib import (
    BeamReconstructor,
    CameraModel,
    CameraModelTolerance,
    MeasurementPlane,
)

# planes: a list of >=3 MeasurementPlane objects built from your own
# flux arrays (see MeasurementPlane below).

# Nominal camera pose/intrinsics from calibration; every field here is
# refined jointly with the mode fit because its tolerance is nonzero.
camera = CameraModel(
    focal_length_px=2000.0,
    position=(0.0, 0.0, -2.0),
    orientation_deg=(0.0, 0.0, 0.0),
)
camera_tolerance = CameraModelTolerance(
    focal_length_px=20.0,
    position=(0.01, 0.01, 0.05),
    orientation_deg=(2.0, 2.0, 2.0),
)

reconstructor = BeamReconstructor(
    w0=5e-3, z0=0.5, wavelength=1.76e-3,
    camera=camera, camera_tolerance=camera_tolerance,
)
result = reconstructor.reconstruct(planes)

for mode, (power_fraction, phase_rad) in result.purity.items():
    print(mode, power_fraction, phase_rad)
```

## `data` — `MeasurementPlane`, `ReconstructionResult`

### `MeasurementPlane(flux, z, z_tolerance=0.0, label=None)`

One measurement: a 2D flux array plus its acquisition metadata.

- `flux` — 2D `np.ndarray` of flux values. Dead-pixel correction, background
  subtraction, and saturation clipping are assumed already handled upstream.
- `z` — nominal distance from the output window, in meters. Must be `> 0`.
- `z_tolerance` — `+/-` bound, in meters, around the nominal `z` within
  which the true distance is jointly refined by `ModalFitter`. Must be
  `>= 0`; `0` (the default) means `z` is trusted exactly and held fixed.
- `label` — optional human-readable identifier.

Per-plane camera geometry (`pixel_scale`/`viewing_angle_deg`) no longer
lives on `MeasurementPlane` — camera pose/intrinsics are a single shared
`CameraModel` for the whole reconstruction (see `geometry` below).

### `validate_planes(planes)`

Raises `ValueError` if there are fewer than 3 planes, planes have
mismatched flux shapes, or `z` values are not all distinct. Called
internally by `ModalFitter.fit`/`fit_auto`, `PhaseRetriever.retrieve`, and
`BeamReconstructor.reconstruct` — you generally don't need to call it
yourself. Not exported from the top-level package; import via
`from he11lib.data import validate_planes` if needed.

### `ReconstructionResult`

Output of a full reconstruction (returned by `ModalFitter.fit`/`fit_auto`
and `BeamReconstructor.reconstruct`):

- `purity: dict[(p, l), (power_fraction, phase_rad)]`
- `reconstructed_field: np.ndarray` — reconstructed complex field.
- `centers: list[(x, y)]` — fitted beam transverse center per plane, meters.
- `pointing_angle_horizontal_deg`, `pointing_angle_vertical_deg: float` —
  fitted shared beam pointing (tilt) angles, independent horizontal and
  vertical.
- `geometry: dict[str, float]` — geometry parameters used or fitted: the 9
  `CameraModel` field names from `he11lib.geometry.CAMERA_FIELD_NAMES`
  (`focal_length_px`, `position_x`, `position_y`, `position_z`, `yaw_deg`,
  `pitch_deg`, `roll_deg`, `principal_point_x`, `principal_point_y`), plus
  `z_{i}` per plane index `i` (that plane's fitted/held distance).
- `residuals: list[np.ndarray]` — per-plane (measured − modeled) flux maps.
  Empty when `used_phase_retrieval` is `True`.
- `coefficient_uncertainty: dict[(p, l), float]` — 1-sigma uncertainty on
  each mode's fitted power fraction. `NaN` per mode when
  `used_phase_retrieval` is `True`.
- `used_phase_retrieval: bool` — whether the phase-retrieval fallback (not
  the modal fit) produced this result.

## `modes` — `LGBasis`

`LGBasis(w0, z0, wavelength)` — the LG mode basis referenced to a known
waist radius `w0` (m), waist location `z0` (m), and radiation `wavelength`
(m).

- `beam_radius(z)` — `w(z)`.
- `inverse_radius_of_curvature(z)` — `1/R(z)` (well-defined, `0`, at the
  waist).
- `gouy_phase(z, p, l)` — Gouy phase of mode `(p, l)` at `z`.
- `field(x, y, z, p, l)` — complex `LG_{p,l}` field sampled on the `(x, y)`
  grid at distance `z`.
- `field_superposition(x, y, z, coefficients)` — complex field for
  `coefficients: dict[(p, l), complex]`.
- `project(complex_field, x, y, dx, z, modes)` — projects `complex_field`
  onto each `(p, l)` in `modes`, returning `dict[(p, l), complex]`
  coefficients (Riemann-sum inner product; `dx` is the grid spacing).

## `geometry` — `CameraModel`, `CameraModelTolerance`, `GeometryCalibration`

### `CameraModel(focal_length_px, position, orientation_deg, principal_point=(0.0, 0.0))`

A nominal pinhole camera pose/intrinsics shared across every plane in one
reconstruction. Always a point estimate — never trusted as exact by
itself; trust is expressed via the paired `CameraModelTolerance`.

- `focal_length_px` — focal length in pixel units.
- `position` — `(x, y, z)` camera position in the beam-axis world frame,
  meters; `z=0` is the output window.
- `orientation_deg` — `(yaw, pitch, roll)`, degrees. All-zero means the
  boresight is normal to every `z=const` target plane with no in-plane
  rotation.
- `principal_point` — `(px, px)` offset from the frame center.

### `CameraModelTolerance(focal_length_px, position, orientation_deg, principal_point=(0.0, 0.0))`

Per-field `+/-` refinement bound, same shape as `CameraModel`. Every field
must be `>= 0` (raises `ValueError` otherwise). A field's tolerance of `0`
holds that `CameraModel` field fixed at its nominal value during fitting;
`> 0` lets `ModalFitter` refine it within `[nominal - tolerance, nominal +
tolerance]`.

### `GeometryCalibration(camera)`

Wraps a `CameraModel` and resolves pixel <-> physical coordinate mappings
via true pinhole projection (not a uniform affine/cosine approximation).

- `pixel_coordinates(x, y, z) -> (row, col)` — forward-projects physical
  `(x, y)` at depth `z` to pixel coordinates. Raises `ValueError` if the
  point is behind the camera (`Z_cam <= 0`).
- `physical_coordinates(image_shape, z) -> (x, y)` — inverse-projects every
  pixel in a frame of `image_shape` to physical `(x, y)` on the `z=const`
  plane, via ray-plane intersection (this is what produces genuine
  keystoning — non-uniform spacing across the frame — for tilted/off-axis
  poses). Raises `ValueError` if the plane is edge-on to or behind the
  camera.
- `effective_pixel_scale(image_shape, z) -> float` — a single isotropic
  meters/pixel figure (finite-difference approximation at the frame
  center), for callers like `DiffusionDeconvolver` that assume one
  isotropic pixel-space kernel.

### `CAMERA_FIELD_NAMES`, `camera_to_values`, `tolerance_to_values`, `camera_from_values`

Module-level helpers used internally by `ModalFitter` to flatten/unflatten
`CameraModel`/`CameraModelTolerance` into the optimizer's parameter vector.
Not usually needed by application code, but exported for advanced use
(e.g. inspecting `CAMERA_FIELD_NAMES` to interpret `ReconstructionResult.geometry` keys).

## `noise` — `NoiseEstimator`

`NoiseEstimator()` — automatic per-image noise estimation (no
user-supplied noise parameter needed).

- `estimate_std(image)` — fast Laplacian-based (Immerkær 1996) noise
  standard-deviation estimate.
- `weights(image)` — per-pixel weights (`1/sigma**2`) for noise-weighted
  least squares.

## `deconvolution` — `DiffusionDeconvolver`

`DiffusionDeconvolver(thermal_diffusivity, dwell_time)` — optional
correction for lateral thermal-diffusion blur in the absorbing target
(`thermal_diffusivity` in m²/s, `dwell_time` in s). Disabled unless you
pass a `deconvolver` to `BeamReconstructor`.

- `blur_sigma_m()` — Gaussian blur standard deviation, in meters.
- `blur(image, pixel_scale)` — forward blur (for synthetic testing).
- `deconvolve(image, pixel_scale, noise_to_signal_ratio=1e-3)` — regularized
  (Wiener) removal of the blur.

Note: the blur/deconvolution kernel is isotropic in pixel space. A tilted
or off-axis `CameraModel` produces a pixel scale that varies across the
frame and between `x`/`y` (keystoning), so `deconvolve` uses
`GeometryCalibration.effective_pixel_scale` — a single isotropic
approximation evaluated at the frame center. This is exact only for an
on-axis, untilted camera; at oblique poses it is an accepted
approximation (see `CLAUDE.md`).

## `synthetic` — `SyntheticBeamGenerator`

`SyntheticBeamGenerator(basis, camera)` — forward model used to validate
the pipeline against known ground truth, and to evaluate experimental
design. `camera` is the ground-truth `CameraModel` (position/orientation/
intrinsics) used to render each plane via true perspective projection.

- `generate(coefficients, z_list, image_shape, *, center=(0.0, 0.0), pointing_angle_horizontal_deg=0.0, pointing_angle_vertical_deg=0.0, z_tolerance=0.0, nominal_z_offsets=None, noise_std=0.0, seed=None) -> list[MeasurementPlane]`
  — returns one `MeasurementPlane` per (true) `z` in `z_list`. The beam's
  transverse center drifts linearly with `z` according to the two
  independent pointing angles, starting from `center` at the basis's
  `z0`. `nominal_z_offsets`, if given, maps a true `z` to an offset
  applied to that plane's *nominal* `z` — letting a reconstruction be
  tested against a deliberately-offset nominal input while the plane's
  flux is still rendered at the true `z`. Every resulting plane shares
  `z_tolerance`.

## `fitting` — `ModalFitter`, `generate_mode_shells`

### `generate_mode_shells(max_order)`

Groups candidate `LG_{p,l}` modes into shells of increasing order
`2p + |l|`, up to and including `max_order`. Returns
`list[list[(p, l)]]`, one list of modes per order.

### `ModalFitter(basis, noise_estimator=None)`

Core reconstruction path: a joint nonlinear least-squares fit of complex LG
coefficients, beam center/pointing, and any nonzero-tolerance camera/`z`
geometry.

- `fit(planes, modes, camera, camera_tolerance, initial_coefficients=None, initial_center=(0.0, 0.0), initial_pointing_deg=(0.0, 0.0)) -> ReconstructionResult`
  — fits exactly the given candidate `modes`. Every `CameraModel` field
  with a nonzero `camera_tolerance` entry, and every plane whose
  `z_tolerance` is nonzero, is refined within `[nominal - tolerance,
  nominal + tolerance]`; zero-tolerance fields are held fixed at their
  nominal value.
- `fit_auto(planes, camera, camera_tolerance, max_order=4, bic_improvement_threshold=10.0) -> ReconstructionResult`
  — starts from `LG_00` and grows the candidate mode set shell-by-shell
  (via `generate_mode_shells`), stopping once BIC no longer improves by
  more than `bic_improvement_threshold`, capped at `max_order`. Emits a
  `UserWarning` (does not raise) if the cap is reached while the fit is
  still improving, or if the number of free camera+`z` parameters is large
  relative to the number of planes (see `CLAUDE.md`'s degeneracy pitfall).

## `phase_retrieval` — `PhaseRetriever`, `propagate_angular_spectrum`

Fallback reconstruction path for when the modal fit's residual stays high,
or when the mode content isn't well described by a small finite mode set.

### `propagate_angular_spectrum(field, dx, dz, wavelength)`

Free-space-propagates a complex `field` (pixel spacing `dx`) by distance
`dz` via the (paraxial) angular-spectrum method — the same propagation
model implicitly assumed by `LGBasis`'s closed-form paraxial modes.

### `PhaseRetriever(wavelength)`

- `retrieve(planes, camera, max_iterations=200) -> PhaseRetrievalResult`
  — multi-plane Gerchberg-Saxton phase retrieval: propagates a trial
  complex field back and forth between planes, enforcing the measured
  amplitude (`sqrt(flux)`) at each plane, without assuming a finite mode
  basis. All planes are propagated on one common physical grid, derived
  from `camera` at the smallest-`z` plane's depth.

### `PhaseRetrievalResult`

`field, x, y, z, center, residual` — the recovered complex field (at the
smallest-`z` plane) on its `(x, y)` grid, the estimated beam center
(intensity centroid), and the final RMS amplitude-mismatch residual.
Project `field` onto `LGBasis` (via `LGBasis.project`) to get a purity
table, as `BeamReconstructor` does internally for its fallback path.

## `reconstruct` — `BeamReconstructor`

`BeamReconstructor(w0, z0, wavelength, camera, camera_tolerance, max_order=4, noise_estimator=None, deconvolver=None, force_phase_retrieval=False, phase_retrieval_residual_threshold=None)`

High-level orchestrator wiring together the full pipeline: optional
diffusion deblurring → `ModalFitter.fit_auto` → optional
`PhaseRetriever` fallback. `camera`/`camera_tolerance` are the nominal
shared `CameraModel` and its per-field refinement bounds for this
reconstruction.

- `reconstruct(planes) -> ReconstructionResult`
  1. Validates `planes` (see `validate_planes`).
  2. If `deconvolver` is set, deblurs each plane using
     `GeometryCalibration(camera).effective_pixel_scale(plane.flux.shape, plane.z)`.
  3. Runs `ModalFitter(basis, noise_estimator).fit_auto(planes, camera, camera_tolerance, max_order)`.
  4. Runs the `PhaseRetriever` fallback instead, projecting its recovered
     field onto all modes up to `max_order`, if `force_phase_retrieval` is
     `True`, or if `phase_retrieval_residual_threshold` is set and the
     modal fit's noise-weighted RMS residual exceeds it. In that case
     `result.residuals` is empty, `coefficient_uncertainty` is `NaN` per
     mode, `geometry` is empty, and both pointing-angle fields are `NaN`
     (phase retrieval doesn't fit geometry/pointing or produce a fit
     covariance).

## `plotting` — diagnostic visualizations

Each function returns a `matplotlib.figure.Figure` for the caller to
display (`fig.show()`) or save (`fig.savefig(...)`); none of them call
`plt.show()` themselves.

- `plot_mode_purity(result)` — bar chart of power fraction per mode.
- `plot_center_trace(planes, result)` — fitted beam center `(x, y)` vs. `z`.
- `plot_residuals(planes, result)` — per-plane residual maps. Raises
  `ValueError` if `result.residuals` is empty (e.g. after the
  phase-retrieval fallback).