diff --git a/CHANGELOG.md b/CHANGELOG.md index 89351e1..6a4c3e8 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/makewfs/sensors/base.py b/src/makewfs/sensors/base.py index 7a27e34..7c42ad6 100644 --- a/src/makewfs/sensors/base.py +++ b/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 diff --git a/src/makewfs/sensors/pyramid.py b/src/makewfs/sensors/pyramid.py index 05d951a..d20d2b7 100644 --- a/src/makewfs/sensors/pyramid.py +++ b/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, ) diff --git a/src/makewfs/sensors/shack_hartmann.py b/src/makewfs/sensors/shack_hartmann.py index b6ea4e0..e17b129 100644 --- a/src/makewfs/sensors/shack_hartmann.py +++ b/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( diff --git a/src/makewfs/wavefront.py b/src/makewfs/wavefront.py index 8dc6ae7..0132623 100644 --- a/src/makewfs/wavefront.py +++ b/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]: diff --git a/tests/test_provenance.py b/tests/test_provenance.py index 9514877..b8129b9 100644 --- a/tests/test_provenance.py +++ b/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 + )