Report phase input in OPD metres

CHANGELOG.md

@@ -4,6 +4,12 @@ All notable changes to `makewfs` are documented here.
 
 ## [Unreleased]
 
+- **Fixed: phase input was reported as metres.** With
+  `input.quantity = "phase"`, `OpticalResult.opd_m` and the
+  `wfs_input_opd_rms_m` frame metadata held the raw phase in radians, so a
+  13 nm wavefront was recorded as 0.083 m. Both now carry the input converted
+  to OPD metres at `input.reference_wavelength_m`. Images were never affected,
+  and OPD input is unchanged.
 - **Changed: generic optics helpers now come from
   [aocore](https://github.com/jacotay7/aocore).** `aocore>=0.1.2,<0.2` moves
   from the `dev` extra to a core dependency. Under CONVENTIONS.md §9, makewfs

src/makewfs/sensors/base.py

@@ -20,6 +20,7 @@ class OpticalResult:
     launched_rate_per_s: float
     captured_rate_per_s: Any
     opd_m: NDArray[np.float64]
+    """Input wavefront in OPD metres (phase input converted), without static OPD."""
     spectral_photon_rate: NDArray[np.float64] | None = None
     spectral_wavelengths_m: tuple[float, ...] | None = None
 

src/makewfs/sensors/pyramid.py

@@ -220,7 +220,7 @@ def render(self, wavefront: NDArray[np.float64]) -> OpticalResult:
             photon_rate,
             self.source_rate,
             captured,
-            wavefront,
+            self.wavefront.input_opd(wavefront),
             spectral_photon_rate,
             self._wavelengths,
         )

src/makewfs/sensors/shack_hartmann.py

@@ -460,7 +460,7 @@ def _render_samples(self, wavefronts: Sequence[NDArray[np.float64]]) -> OpticalR
         internal = (
             internal_samples[0] if sample_count == 1 else self.backend.stack(internal_samples)
         )
-        opd = (
+        opd = self.wavefront.input_opd(
             wavefronts[0]
             if sample_count == 1
             else self.backend.mean(

src/makewfs/wavefront.py

@@ -83,6 +83,23 @@ def __init__(
         )
         self._resample_coordinates: dict[tuple[int, int], Any] = {}
 
+    def input_opd(self, value: ArrayLike) -> NDArray[np.float64]:
+        """The input wavefront in OPD metres, before static OPD and regridding.
+
+        OPD input is returned as given. Phase input is converted at the
+        configured reference wavelength.
+        """
+        if self.config.input.quantity != "phase":
+            return cast(NDArray[np.float64], value)
+        assert self.config.input.reference_wavelength_m is not None
+        return cast(
+            NDArray[np.float64],
+            phase_to_opd(
+                self.backend.asarray(value, dtype=np.float64),
+                self.config.input.reference_wavelength_m,
+            ),
+        )
+
     def opd(
         self, value: ArrayLike, *, target_shape: tuple[int, int] | None = None
     ) -> NDArray[np.float64]:

tests/test_provenance.py

@@ -4,9 +4,12 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
 
 from makewfs import WavefrontSensor, load_config
 
+CONFIGS = Path(__file__).parents[1] / "examples" / "configs"
+
 
 def test_source_curve_digest_is_recorded(tmp_path: Path) -> None:
     path = tmp_path / "sed.txt"
@@ -48,3 +51,32 @@ def test_optical_blur_kernel_digest_is_recorded(tmp_path: Path) -> None:
         np.zeros(config.input.shape), seed=4
     )
     assert len(frame.metadata["wfs_shack_hartmann_optical_blur_kernel_sha256"]) == 16
+
+
+@pytest.mark.parametrize("name", ["shack_hartmann_minimal.toml", "pyramid_minimal.toml"])
+def test_phase_input_reports_its_opd_in_metres(name: str) -> None:
+    # Phase input is converted at the reference wavelength before anything
+    # downstream sees it, including the input-RMS metadata and OpticalResult.
+    opd_sensor = WavefrontSensor.from_toml(CONFIGS / name)
+    wavelength = 1.0e-6
+    phase_input = replace(
+        opd_sensor.config.input,
+        quantity="phase",
+        unit="rad",
+        reference_wavelength_m=wavelength,
+    )
+    phase_sensor = WavefrontSensor(replace(opd_sensor.config, input=phase_input))
+    shape = opd_sensor.config.input.shape
+    _, xx = np.indices(shape, dtype=np.float64)
+    opd = 30e-9 * (xx - (shape[1] - 1) / 2) / shape[1] + 10e-9
+    phase = 2.0 * np.pi * opd / wavelength
+
+    expected = float(np.sqrt(np.mean(opd**2)))
+    from_opd = opd_sensor.expose(opd, seed=1).metadata["wfs_input_opd_rms_m"]
+    from_phase = phase_sensor.expose(phase, seed=1).metadata["wfs_input_opd_rms_m"]
+    assert from_opd == pytest.approx(expected, rel=1e-12)
+    assert from_phase == pytest.approx(expected, rel=1e-12)
+    np.testing.assert_allclose(phase_sensor._render(phase).opd_m, opd, rtol=1e-12)
+    np.testing.assert_allclose(
+        phase_sensor.photon_rate(phase), opd_sensor.photon_rate(opd), rtol=1e-6, atol=0.0
+    )