Import generic optics helpers from aocore

AGENTS.md

@@ -62,6 +62,13 @@ pupil OPD/phase + static config
   because it is part of WFS image formation, but it never predicts LGS return
   flux or evolves the sodium layer.
 
+- Stack conventions and generic optics primitives belong to
+  [`aocore`](https://github.com/jacotay7/aocore), a core dependency. Its
+  `CONVENTIONS.md` is the contract, and §9 says to import the primitives rather
+  than re-implement them. Today makewfs takes `centered_coordinates`,
+  `ARCSEC_TO_RAD`/`RAD_TO_ARCSEC`, `opd_to_phase`/`phase_to_opd` and
+  `block_sum` from it. Change a convention in aocore, never here.
+
 Do not copy sibling physics for convenience. If their public API is insufficient,
 write a failing integration test/design note, use the conditional gates in
 `ROADMAP.md`, and make the smallest change in the owning repository.
@@ -141,6 +148,13 @@ Follow the target layout in `ROADMAP.md`:
 
 - `config.py` parses and validates; it does not propagate optics.
 - `wavefront.py`, `pupil.py`, and `sampling.py` hold shared numerical rules.
+  Some aocore look-alikes stay local on purpose; each says why in its module
+  docstring. `pupil.make_pupil` builds wedge spiders, square segment gaps,
+  rotation and rectangular grids on the selected backend, none of which
+  `aocore.Pupil` models. `backend.ArrayBackend` takes a dtype per array, explicit
+  FFT workers and `ndimage` helpers, which `aocore.Backend` does not, and its
+  centred FFTs keep the `fftshift` convention. `sampling.block_sum` wraps
+  `aocore.block_sum` but keeps its own factor-two fast path.
 - `sensors/` contains deterministic ideal optical engines and no camera noise.
   `_shack_hartmann_cuda.py` is a private first-use-JIT execution plan for exact
   compatible CUDA geometries; `shack_hartmann.py` remains the readable physics

CHANGELOG.md

@@ -4,6 +4,24 @@ All notable changes to `makewfs` are documented here.
 
 ## [Unreleased]
 
+- **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
+  now imports these primitives instead of keeping its own copies:
+  `centered_coordinates` for the pupil, field-stop and DFT detector grids;
+  `ARCSEC_TO_RAD`/`RAD_TO_ARCSEC` for source angles and the Shack-Hartmann
+  plate scale; `opd_to_phase`/`phase_to_opd` for phase input and the sensor
+  phasors; and `block_sum` for pixel integration. No public name changes.
+  `makewfs.sampling.block_sum` keeps its signature, error messages and
+  factor-two fast path, and delegates other factors to aocore. The behaviour
+  already followed the conventions. Float32 renders of the example
+  configurations are bit-for-bit unchanged, and float64 renders differ only by
+  floating-point rounding (at most about 6e-16 relative), from the changed
+  operation order in the angle and phase conversions and, for oversampling
+  factors above two, in pixel binning. The pupil masks, the
+  `ArrayBackend` and its `fftshift`-convention centred FFTs stay local
+  because aocore has no equivalent for them; the module docstrings say why.
+
 - **Tests: conformance with the AO stack conventions.** `tests/test_conformance.py`
   runs [aocore](https://github.com/jacotay7/aocore)'s checks against both
   sensors. Shack-Hartmann spots centre between pixels for a flat wavefront and

README.md

@@ -22,7 +22,9 @@ quadrature, and sodium-layer LGS geometry — everything from the pupil to the A
 The package owns the wavefront-sensor optics and nothing else. It deliberately
 reuses [`pyturb`](https://github.com/jacotay7/pyturb) for atmospheric OPD and
 [`getframes`](https://github.com/jacotay7/getframes) for detector response and
-noise; neither model is reimplemented here. It runs on NumPy by default and
+noise; neither model is reimplemented here. Shared conventions and generic optics
+primitives (coordinates, angle constants, OPD/phase conversion, pixel binning)
+come from [`aocore`](https://github.com/jacotay7/aocore). It runs on NumPy by default and
 switches to CUDA (via CuPy) with a single configuration field.
 
 ## Install

docs/units-and-coordinates.md

@@ -1,14 +1,24 @@
 # Units and coordinates
 
+`makewfs` follows the AO stack conventions defined by
+[`aocore`](https://github.com/jacotay7/aocore) (`CONVENTIONS.md`), and
+`tests/test_conformance.py` runs aocore's checks against both sensors. Pixel
+centres, arcsecond/radian constants, OPD/phase conversion and flux-conserving
+pixel binning come from `aocore` rather than local copies.
+
 - Input arrays are `(y, x)` NumPy arrays.
 - OPD is metres. Phase is radians and requires
   `input.reference_wavelength_m`.
 - `input.grid_extent_m` is the physical width/height represented by the input
   array; coordinates are centred on pixel centres.
-- `source.field_angle_arcsec` is `[x, y]` on-sky angle.
+- `source.field_angle_arcsec` is `[x, y]` on-sky angle, converted with
+  `aocore.ARCSEC_TO_RAD`.
 - `lgs_launch_position_m` is `[x, y]` in the entrance-pupil plane.
 - Optical output is photons/s/native detector pixel. Detector output is ADU.
 
-The centred FFT helpers use unitary normalization. Piston removal is a numerical
+The centred FFT helpers use unitary normalization. They keep the `fftshift`
+convention (zero frequency on pixel `n // 2`) internally; the detector images
+the sensors produce follow the stack's centring, with the optical axis at pixel
+`(n - 1) / 2`, which the conformance tests check. Piston removal is a numerical
 stabilization only; it does not alter intensity. Detector cropping reports lost
 flux rather than renormalizing the image.

pyproject.toml

@@ -25,6 +25,7 @@ classifiers = [
     "Typing :: Typed",
 ]
 dependencies = [
+    "aocore>=0.1.2,<0.2",
     "numpy>=1.23",
     "scipy>=1.10",
     "getframes>=2.2.0",
@@ -33,7 +34,6 @@ dependencies = [
 
 [project.optional-dependencies]
 dev = [
-    "aocore>=0.1.1",
     "build>=1.0",
     "mypy>=1.8",
     "pytest>=7.0",

src/makewfs/api.py

@@ -8,6 +8,7 @@
 from typing import Any, cast
 
 import numpy as np
+from aocore import RAD_TO_ARCSEC
 from numpy.typing import ArrayLike, NDArray
 
 from .backend import ArrayBackend, cpu_backend, cupy_backend
@@ -20,8 +21,6 @@
 from .sensors.shack_hartmann import ShackHartmannEngine
 from .wavefront import iter_phase_samples
 
-_ARCSEC_PER_RADIAN = 206264.80624709636
-
 
 class WavefrontSensor:
     """Configured wavefront-sensor facade.
@@ -125,7 +124,7 @@ def subaperture_plate_scale_arcsec(self) -> float:
         sampling = engine._spot_sampling(self.config.sensor.wavelength_m)
         subaperture_m = self.config.input.grid_extent_m / engine.n_lenslets
         lambda_over_d_rad = self.config.sensor.wavelength_m / subaperture_m
-        return float(lambda_over_d_rad / sampling * _ARCSEC_PER_RADIAN)
+        return float(lambda_over_d_rad / sampling * RAD_TO_ARCSEC)
 
     def subaperture_field_of_view_arcsec(self) -> float:
         """Return the on-sky field of view one subaperture's window spans.

src/makewfs/backend.py

@@ -5,6 +5,14 @@
 NumPy directly.  That boundary is deliberately private today; it gives a
 future CuPy implementation one place to provide array creation, reductions,
 FFT, and interpolation semantics without changing the optical equations.
+
+It is deliberately separate from ``aocore.Backend``. The sensors choose a dtype
+per array rather than one precision per backend, pass explicit FFT worker
+counts and ``overwrite_x``, and need SciPy/CuPy ``ndimage`` helpers. The
+centred FFTs here use the ``fftshift`` convention (zero frequency on pixel
+``n // 2``). The engines own the conversion to CONVENTIONS 1.3 centring: for
+example ``sampling.spot_intensity`` applies a half-sample phase ramp for even
+detectors, and ``tests/test_conformance.py`` checks both sensors.
 """
 
 from __future__ import annotations

src/makewfs/pupil.py

@@ -1,4 +1,12 @@
-"""Analytic and custom entrance-pupil masks."""
+"""Analytic and custom entrance-pupil masks.
+
+These stay local rather than using ``aocore.Pupil``: the configured geometry
+has no aocore equivalent. Spiders are angular wedges of ``width_fraction`` of
+pi, segment gaps lie on a square grid, the whole pattern rotates with
+``pupil_rotation_deg``, grids may be rectangular, and masks are built on the
+selected array backend in the working precision. Pixel centres follow
+CONVENTIONS 1.2 through ``wavefront._coordinates``.
+"""
 
 from __future__ import annotations
 

src/makewfs/sampling.py

@@ -8,6 +8,8 @@
 from typing import Any, cast
 
 import numpy as np
+from aocore import block_sum as _aocore_block_sum
+from aocore import centered_coordinates
 from numpy.typing import NDArray
 
 from .backend import ArrayBackend, centered_fft_intensity, cpu_backend
@@ -142,7 +144,13 @@ def crop_center(array: NDArray[Any], shape: tuple[int, int]) -> NDArray[Any]:
 def block_sum(
     array: NDArray[Any], factor: int, *, backend: ArrayBackend | None = None
 ) -> NDArray[Any]:
-    """Sum square pixel blocks while preserving total flux."""
+    """Sum square pixel blocks while preserving total flux.
+
+    A thin wrapper over ``aocore.block_sum``, which owns flux-conserving
+    binning (CONVENTIONS 9). It keeps this module's error messages and a
+    factor-two fast path. ``backend`` is accepted for compatibility; the
+    reduction runs on the array's own namespace.
+    """
     if factor < 1:
         raise ValueError("factor must be positive")
     if factor == 1:
@@ -161,10 +169,7 @@ def block_sum(
             + array[..., 1::2, 0::2]
             + array[..., 1::2, 1::2],
         )
-    resolved = backend or cpu_backend()
-    reshaped = array.reshape((*array.shape[:-2], height // factor, factor, width // factor, factor))
-    reduced_x = resolved.sum(reshaped, axis=-1)
-    return cast(NDArray[Any], resolved.sum(reduced_x, axis=-2))
+    return cast(NDArray[Any], _aocore_block_sum(array, factor))
 
 
 @dataclass
@@ -245,9 +250,8 @@ def build(
                 coordinate = backend.arange(nfft, dtype=np.float64)
                 plan.half_sample = backend.exp(-1j * math.pi * coordinate / nfft)
         else:
-            detector_coordinate = (
-                backend.arange(high_resolution_pixels, dtype=np.float64)
-                - (high_resolution_pixels - 1) / 2.0
+            detector_coordinate = backend.asarray(
+                centered_coordinates(high_resolution_pixels), dtype=np.float64
             ) / (sampling * oversampling)
             pupil_coordinate = backend.arange(samples_per_lenslet, dtype=np.float64)
             kernel = backend.exp(
@@ -259,12 +263,11 @@ def build(
             )
             plan.dft_kernel = backend.astype(kernel, plan.field_dtype)
         if field_stop_radius_lambda_over_d is not None:
-            coordinates = backend.arange(high_resolution_pixels, dtype=np.float64)
+            coordinates = backend.asarray(
+                centered_coordinates(high_resolution_pixels), dtype=np.float64
+            )
             y, x = backend.meshgrid(coordinates, coordinates, indexing="ij")
-            radius_lambda_over_d = backend.hypot(
-                x - (high_resolution_pixels - 1) / 2.0,
-                y - (high_resolution_pixels - 1) / 2.0,
-            ) / (oversampling * sampling)
+            radius_lambda_over_d = backend.hypot(x, y) / (oversampling * sampling)
             plan.field_stop_mask = radius_lambda_over_d <= field_stop_radius_lambda_over_d
         return plan
 

src/makewfs/sensors/pyramid.py

@@ -6,6 +6,7 @@
 from typing import Any, cast
 
 import numpy as np
+from aocore import opd_to_phase
 from numpy.typing import NDArray
 
 from ..backend import ArrayBackend, centered_fft2, centered_ifft2, complex_dtype, real_dtype
@@ -153,7 +154,7 @@ def _base_phasor(
         total_opd = internal + self._field_angle_opd[state_index]
         piston = total_opd[self._piston_index]
         relative_opd = total_opd - piston
-        phase = 2.0 * math.pi * relative_opd / state.wavelength_m
+        phase = opd_to_phase(relative_opd, state.wavelength_m)
         return cast(
             NDArray[Any],
             self.backend.asarray(self.backend.exp(1j * phase), dtype=self._complex_dtype),

src/makewfs/sensors/shack_hartmann.py

@@ -7,6 +7,7 @@
 from typing import Any, cast
 
 import numpy as np
+from aocore import opd_to_phase
 from numpy.typing import NDArray
 
 from ..backend import ArrayBackend, complex_dtype, real_dtype
@@ -379,7 +380,7 @@ def _field(
         # it makes the exact physical piston invariance survive finite precision
         # in exp(i*phase) and has no effect on the intensity.
         relative_opd = total_opd - piston
-        phase = 2.0 * math.pi * relative_opd / state.wavelength_m
+        phase = opd_to_phase(relative_opd, state.wavelength_m)
         return cast(
             NDArray[Any],
             self.backend.asarray(
@@ -416,7 +417,7 @@ def _fields(
             [self.source_states[index].wavelength_m for index in state_group] * samples,
             dtype=self._real_dtype,
         )[:, None, None]
-        phase = 2.0 * math.pi * (total_opd - piston) / wavelengths
+        phase = opd_to_phase(total_opd - piston, wavelengths)
         return cast(
             NDArray[Any],
             self.backend.asarray(

src/makewfs/source.py

@@ -7,6 +7,7 @@
 from pathlib import Path
 
 import numpy as np
+from aocore import ARCSEC_TO_RAD
 from numpy.typing import NDArray
 
 from .config import WFSConfig
@@ -33,8 +34,8 @@ def _normalised(values: tuple[float, ...], count: int) -> NDArray[np.float64]:
 def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]:
     """Return ``(x angle, y angle, weight)`` Gaussian angular quadrature."""
     source = config.source
-    base_x = math.radians(source.field_angle_arcsec[0] / 3600.0)
-    base_y = math.radians(source.field_angle_arcsec[1] / 3600.0)
+    base_x = source.field_angle_arcsec[0] * ARCSEC_TO_RAD
+    base_y = source.field_angle_arcsec[1] * ARCSEC_TO_RAD
     if source.angular_kernel_path is not None:
         data = np.atleast_2d(
             np.asarray(np.loadtxt(Path(source.angular_kernel_path)), dtype=np.float64)
@@ -48,8 +49,8 @@ def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]:
         weights = data[:, 2] / np.sum(data[:, 2])
         return [
             (
-                base_x + math.radians(float(x) / 3600.0),
-                base_y + math.radians(float(y) / 3600.0),
+                base_x + float(x) * ARCSEC_TO_RAD,
+                base_y + float(y) * ARCSEC_TO_RAD,
                 float(weight),
             )
             for (x, y), weight in zip(data[:, :2], weights)
@@ -60,7 +61,7 @@ def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]:
     # Three-point Gauss-Hermite is exact through fourth order for a normal
     # expectation. For higher configured orders use an evenly weighted grid;
     # this keeps the rule deterministic without requiring an optional package.
-    sigma = math.radians(source.angular_fwhm_arcsec / 3600.0) / 2.3548200450309493
+    sigma = source.angular_fwhm_arcsec * ARCSEC_TO_RAD / 2.3548200450309493
     if order == 3:
         nodes = np.array([-math.sqrt(3.0), 0.0, math.sqrt(3.0)])
         one_d = np.array([1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0])

src/makewfs/wavefront.py

@@ -7,6 +7,7 @@
 from typing import Any, cast
 
 import numpy as np
+from aocore import centered_coordinates, phase_to_opd
 from numpy.typing import ArrayLike, NDArray
 
 from .backend import ArrayBackend, cpu_backend
@@ -20,11 +21,16 @@ def _coordinates(
     backend: ArrayBackend | None = None,
     dtype: object = np.float64,
 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
-    """Return centered physical ``(x, y)`` coordinates for an array shape."""
+    """Return centered physical ``(x, y)`` coordinates for an array shape.
+
+    Pixel centres follow ``aocore.centered_coordinates`` (CONVENTIONS 1.2). The
+    unit-pitch offsets are exact in either precision; scaling to metres happens
+    on the backend in ``dtype``.
+    """
     resolved = backend or cpu_backend()
     height, width = shape
-    x = (resolved.arange(width, dtype=dtype) - (width - 1) / 2.0) * extent_m / width
-    y = (resolved.arange(height, dtype=dtype) - (height - 1) / 2.0) * extent_m / height
+    x = resolved.asarray(centered_coordinates(width), dtype=dtype) * extent_m / width
+    y = resolved.asarray(centered_coordinates(height), dtype=dtype) * extent_m / height
     xx, yy = resolved.meshgrid(x, y)
     return xx, yy
 
@@ -91,7 +97,7 @@ def opd(
         converted = self.backend.asarray(array, dtype=np.float64)
         if self.config.input.quantity == "phase":
             assert self.config.input.reference_wavelength_m is not None
-            converted = converted * self.config.input.reference_wavelength_m / (2.0 * np.pi)
+            converted = phase_to_opd(converted, self.config.input.reference_wavelength_m)
         if self.static_opd is not None:
             converted = converted + self.static_opd
         if not self.backend.scalar(self.backend.all(self.backend.isfinite(converted))):

tests/test_conformance.py

@@ -3,11 +3,10 @@
 from __future__ import annotations
 
 import numpy as np
-import pytest
+from aocore import ARCSEC_TO_RAD, conformance
 
 from makewfs import WavefrontSensor, load_config
 
-conformance = pytest.importorskip("aocore.conformance")
 SH_CONFIG = "examples/configs/shack_hartmann_minimal.toml"
 PYRAMID_CONFIG = "examples/configs/pyramid_minimal.toml"
 
@@ -19,8 +18,6 @@ def _sensor(path: str):
 
 
 def test_shack_hartmann_spots_centre_between_pixels_and_follow_tilt() -> None:
-    from aocore import ARCSEC_TO_RAD
-
     config, sensor, image = _sensor(SH_CONFIG)
     shape = config.input.shape
     conformance.check_image_centring(image, pupil_shape=shape)