From 9edd7b6757c5a179e5323a93cb94a3209d137241 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Wed, 10 Jun 2026 21:44:46 -0600 Subject: [PATCH 1/2] Added support for using the object surface as the field stop of a SequentialSystem Marking the object surface (with an angular, dimensionless aperture) as the field stop is the natural configuration for dispersive systems, where the sensor outline is not reachable at a single wavelength, but the stop root-finding problem was broken in three ways for this case: - The position-variable branch of `_calc_rayfunction_stops_only` called `sag()` with a 2-component vector, raising a TypeError before the solve even started. - The solved outputs were never broadcast over both stop axes, so the reductions in `field_min`/`field_max`/`pupil_min`/`pupil_max` raised an axis error. - The initial guess started every ray at the origin of the first stop with direction +z, which produces NaN residuals on the first iteration for surfaces far from that axis (e.g. the off-axis feed mirror of a Rowland-circle spectrograph) or on steep grazing-incidence flanks, and Newton's method cannot recover from NaN. The seed now aims each ray at its own target point on the last stop surface when no surface with optical power lies between the two stops (nearly exact), and otherwise at the center of the first powered surface. Adds a grazing-incidence spectrograph regression test whose source is the field stop, asserting that the recovered field equals the source's angular radius. Co-Authored-By: Claude Fable 5 --- optika/systems/_sequential.py | 99 +++++++++++++++++- optika/systems/_sequential_test.py | 162 +++++++++++++++++++++++++++++ 2 files changed, 257 insertions(+), 4 deletions(-) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index a0dbbbe4..d2f0f1cc 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -178,6 +178,26 @@ def pupil_stop(self) -> optika.surfaces.AbstractSurface: """ return self.surfaces_all[self.index_pupil_stop] + @staticmethod + def _anchor_surface( + subsystem: list[optika.surfaces.AbstractSurface], + ) -> optika.surfaces.AbstractSurface: + """ + The first surface after the start of the given subsystem with optical + power (a mirror, a curved sag, or rulings), used to aim the initial + guess of the stop root-finding problem. + """ + for surface in subsystem[1:]: + material = surface.material + if material is not None and material.is_mirror: + return surface + sag = surface.sag + if sag is not None and not isinstance(sag, optika.sags.NoSag): + return surface + if surface.rulings is not None: + return surface + return subsystem[~0] + @classmethod def _ray_error( cls, @@ -288,7 +308,7 @@ def _calc_rayfunction_stops_only( result.outputs.direction = na.Cartesian3dVectorArray(0, 0, 1) component_variable = "direction" - def zfunc(xy: na.AbstractCartesian3dVectorArray): + def zfunc(xy: na.AbstractCartesian2dVectorArray): return np.sqrt(1 - np.square(xy.length)) elif na.unit(grid_first).is_equivalent(u.dimensionless_unscaled): @@ -300,12 +320,66 @@ def zfunc(xy: na.AbstractCartesian3dVectorArray): result.outputs.position = na.Cartesian3dVectorArray() * u.mm component_variable = "position" - def zfunc(xy: na.AbstractCartesian3dVectorArray): - return surface_first.sag(xy) + def zfunc(xy: na.AbstractCartesian2dVectorArray): + position = na.Cartesian3dVectorArray( + x=xy.x, + y=xy.y, + z=0 * na.unit_normalized(xy.x), + ) + return surface_first.sag(position) else: raise ValueError(f"unrecognized input grid unit, {na.unit(grid_first)}") + # Seed the free ray component by aiming each ray at its own + # target point on the last stop surface when no surface with + # optical power lies between the two stops (the seed is then + # nearly exact), and otherwise at the center of the first powered + # surface, so that the root-finding starts inside its basin of + # convergence even for surfaces far from the axis of the first + # stop (e.g. an off-axis fold or feed mirror). + anchor = self._anchor_surface(subsystem) + if anchor is surface_last and na.unit(grid_last).is_equivalent(u.mm): + aim = na.Cartesian3dVectorArray( + x=grid_last.x, + y=grid_last.y, + z=0 * na.unit_normalized(grid_last.x), + ) + aim.z = surface_last.sag(aim) + if surface_last.transformation is not None: + aim = surface_last.transformation(aim) + else: + aim = na.Cartesian3dVectorArray() * u.mm + if anchor.transformation is not None: + aim = anchor.transformation(aim) + if surface_first.transformation is not None: + aim = surface_first.transformation.inverse(aim) + + if component_variable == "direction": + d = aim - result.outputs.position + d = d / d.length + # the sign of the seed direction is irrelevant to the + # root-finding problem (surface intercepts may have negative + # distance), but the z-component must be positive to be + # consistent with `zfunc` + flip = np.sign(d.z) + where = d.z != 0 + result.outputs.direction = na.Cartesian3dVectorArray( + x=np.where(where, flip * d.x, 0), + y=np.where(where, flip * d.y, 0), + z=np.where(where, flip * d.z, 1), + ) + else: + d = result.outputs.direction + t = aim.z / d.z + position_seed = na.Cartesian3dVectorArray( + x=aim.x - d.x * t, + y=aim.y - d.y * t, + z=0 * u.mm, + ) + position_seed.z = surface_first.sag(position_seed) + result.outputs.position = position_seed + if surface_first.transformation is not None: result.outputs = surface_first.transformation(result.outputs) @@ -347,7 +421,13 @@ def zfunc(xy: na.AbstractCartesian3dVectorArray): # measured in physical units, yielding a Jacobian made of noise. # Use a perturbation proportional to the scale of the problem # instead. - dx = 1e-6 + if component_variable == "direction": + dx = 1e-6 + else: + dx = 1e-6 * np.maximum( + scale, + 1 * na.unit_normalized(scale), + ) def jacobian(x, _function=function, _dx=dx): return na.jacobian(function=_function, x=x, dx=_dx) @@ -419,6 +499,17 @@ def _calc_rayfunction_stops( where = rays.direction @ obj.sag.normal(rays.position) > 0 result.outputs.direction[where] = -result.outputs.direction[where] + # If the first stop is the object surface, the solved variable is the + # position and the direction retains only the field-stop axis, so + # broadcast both components against each other to guarantee that + # reductions over both stop axes are well-defined downstream. + shape = na.shape_broadcasted( + result.outputs.position, + result.outputs.direction, + ) + result.outputs.position = result.outputs.position.broadcast_to(shape) + result.outputs.direction = result.outputs.direction.broadcast_to(shape) + if self.transformation is not None: result.outputs = self.transformation(result.outputs) diff --git a/optika/systems/_sequential_test.py b/optika/systems/_sequential_test.py index 8b129765..2143794a 100644 --- a/optika/systems/_sequential_test.py +++ b/optika/systems/_sequential_test.py @@ -268,3 +268,165 @@ def test_spot_diagram(self, a: optika.systems.AbstractSequentialSystem): ) class TestSequentialSystem(AbstractTestAbstractSequentialSystem): pass + + +def test__anchor_surface(): + first = optika.surfaces.Surface(name="first") + last = optika.surfaces.Surface(name="last") + mirror = optika.surfaces.Surface( + name="mirror", + material=optika.materials.Mirror(), + ) + curved = optika.surfaces.Surface( + name="curved", + sag=optika.sags.SphericalSag(radius=-100 * u.mm), + ) + grating = optika.surfaces.Surface( + name="grating", + rulings=optika.rulings.Rulings(spacing=1 * u.um, diffraction_order=1), + ) + flat = optika.surfaces.Surface(name="flat") + + anchor = optika.systems.SequentialSystem._anchor_surface + assert anchor([first, flat, mirror, last]) is mirror + assert anchor([first, curved, last]) is curved + assert anchor([first, grating, last]) is grating + assert anchor([first, flat, last]) is last + + +# small enough that the image of the field fits on the sensor +_radius_field_newtonian = 0.05 * u.deg + +_system_newtonian = optika.systems.SequentialSystem( + object=optika.surfaces.Surface( + name="source", + aperture=optika.apertures.CircularAperture( + radius=np.sin(_radius_field_newtonian), + ), + is_field_stop=True, + ), + surfaces=[ + optika.surfaces.Surface( + name="primary", + sag=optika.sags.SphericalSag(radius=-2000 * u.mm), + material=optika.materials.Mirror(), + aperture=optika.apertures.CircularAperture(radius=50 * u.mm), + transformation=na.transformations.Cartesian3dTranslation( + z=500 * u.mm, + ), + ), + optika.surfaces.Surface( + name="aperture", + aperture=optika.apertures.CircularAperture(radius=10 * u.mm), + transformation=na.transformations.Cartesian3dTranslation( + z=250 * u.mm, + ), + is_pupil_stop=True, + ), + ], + sensor=optika.sensors.ImagingSensor( + name="sensor", + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + timedelta_exposure=1 * u.s, + num_pixel=na.Cartesian2dVectorArray(128, 128), + transformation=na.transformations.Cartesian3dTranslation( + z=-500 * u.mm, + ), + ), + grid_input=_grid_input, +) + + +@pytest.mark.parametrize(argnames="a", argvalues=[_system_newtonian]) +class TestSequentialSystemNewtonian( + AbstractTestAbstractSequentialSystem, +): + """ + A Newtonian-style telescope where the pupil stop is downstream of the + primary mirror, so that the initial guess of the stop root-finding + problem must be aimed at the center of the primary instead of directly + at its own target on the pupil stop. + """ + + def test_field_max_matches_source_aperture( + self, + a: optika.systems.AbstractSequentialSystem, + ): + result = a.field_max + assert np.abs(result.x - _radius_field_newtonian) < 1e-6 * u.deg + assert np.abs(result.y - _radius_field_newtonian) < 1e-6 * u.deg + + +_radius_field_grazing = 0.25 * u.deg + +_system_grazing = optika.systems.SequentialSystem( + object=optika.surfaces.Surface( + name="source", + aperture=optika.apertures.CircularAperture( + radius=np.sin(_radius_field_grazing), + ), + is_field_stop=True, + ), + surfaces=[ + optika.surfaces.Surface( + name="paraboloid", + sag=optika.sags.ParabolicSag(focal_length=-2000 * u.mm), + material=optika.materials.Mirror(), + aperture=optika.apertures.CircularAperture(radius=260 * u.mm), + transformation=na.transformations.Cartesian3dTranslation( + z=2500 * u.mm, + ), + is_pupil_stop=True, + ), + optika.surfaces.Surface( + name="grating", + rulings=optika.rulings.Rulings( + spacing=10 * u.um, + diffraction_order=1, + ), + aperture=optika.apertures.RectangularAperture( + half_width=60 * u.mm, + ), + transformation=na.transformations.Cartesian3dTranslation( + z=1000 * u.mm, + ), + ), + ], + sensor=optika.sensors.ImagingSensor( + name="sensor", + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + # short exposure so that the Poisson lam stays representable for the + # large collecting area of the grazing primary + timedelta_exposure=1 * u.us, + num_pixel=na.Cartesian2dVectorArray(2048, 1024), + # offset by the deflection of the first diffraction order, + # (z_grating - z_sensor) * wavelength / spacing + transformation=na.transformations.Cartesian3dTranslation( + x=26 * u.mm, + z=480 * u.mm, + ), + ), + grid_input=_grid_input, +) + + +@pytest.mark.parametrize(argnames="a", argvalues=[_system_grazing]) +class TestSequentialSystemGrazingSpectrograph( + AbstractTestAbstractSequentialSystem, +): + """ + A grazing-incidence spectrograph with a transmission grating, where the + object surface (with an angular aperture) is the field stop. This guards + against regressions in the object-as-field-stop code path of the stop + root-finding problem. + """ + + def test_field_max_matches_source_aperture( + self, + a: optika.systems.AbstractSequentialSystem, + ): + result = a.field_max + assert np.abs(result.x - _radius_field_grazing) < 1e-6 * u.deg + assert np.abs(result.y - _radius_field_grazing) < 1e-6 * u.deg From 2a93b2d241c24de8d10dd3c267ff326a11858a2d Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Wed, 10 Jun 2026 21:48:23 -0600 Subject: [PATCH 2/2] Fixed the stop root-finding problem for transmissive stop surfaces MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The reversed trace in `_calc_rayfunction_stops` propagates the solved stop rays back to the object by re-applying every surface from the stop backward. That is only correct when the stop is a mirror: reflection is an involution, so applying the stop surface to its own forward output turns the rays around exactly. A transmissive stop that modifies the rays (e.g. a transmission grating or a Fresnel zone plate) is not an involution — its diffraction was applied a second time instead of undone, which corrupted the computed field and pupil extents by an order of magnitude. For non-involutory stops, the stop's own diffraction/refraction is now included inside the root-finding problem (the solved variable becomes the pre-stop ray), and the stop's material and rulings are stripped on the reversed leg so only its geometry participates. A true time-reversed trace through transmissive rulings would require negating the diffraction order (reflective rulings are time-reversal symmetric with the same order); this approach sidesteps reversed traces through transmissive elements entirely. Adds a Fresnel-zone-plate regression test asserting that the recovered field half-angle matches arctan(sensor half-width / focal length). Co-Authored-By: Claude Fable 5 --- optika/systems/_sequential.py | 63 ++++++++++++++++++++-- optika/systems/_sequential_test.py | 85 ++++++++++++++++++++++++++++++ 2 files changed, 144 insertions(+), 4 deletions(-) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index d2f0f1cc..269f6d41 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -178,6 +178,28 @@ def pupil_stop(self) -> optika.surfaces.AbstractSurface: """ return self.surfaces_all[self.index_pupil_stop] + @staticmethod + def _stop_is_involutory( + surface: optika.surfaces.AbstractSurface, + ) -> bool: + """ + Whether applying the given stop surface to its own forward output + exactly undoes its effect on the rays. This is true for mirrors + (reflection is an involution) and for surfaces that don't modify the + rays at all (e.g. a vacuum object surface), but false for surfaces + with rulings or refraction, whose effect would be applied a second + time instead of undone. + """ + material = surface.material + if material is not None and material.is_mirror: + return True + if surface.rulings is not None: + return False + if material is not None: + if not isinstance(material, optika.materials.Vacuum): + return False + return True + @staticmethod def _anchor_surface( subsystem: list[optika.surfaces.AbstractSurface], @@ -203,7 +225,8 @@ def _ray_error( cls, a: na.Cartesian2dVectorArray, rays: optika.rays.RayVectorArray, - subsystem: list[optika.surfaces.AbstractSurface], + propagators: list[optika.surfaces.AbstractSurface], + transformation_last: None | na.transformations.AbstractTransformation, grid_last: na.Cartesian2dVectorArray, component_variable: str, component_target: str, @@ -215,11 +238,10 @@ def _ray_error( rays_component_variable.z = zfunc(a) rays = optika.propagators.propagate_rays( - propagators=subsystem[1:], + propagators=propagators, rays=rays, ) - transformation_last = subsystem[~0].transformation if transformation_last is not None: rays = transformation_last.inverse(rays) @@ -390,6 +412,18 @@ def zfunc(xy: na.AbstractCartesian2dVectorArray): else: raise ValueError(f"unrecognized output grid unit, {na.unit(grid_last)}") + # A mirror stop reverses the rays, so the solved variable can be + # the post-reflection ray and the stop surface itself can be + # excluded from the root-finding propagators. A transmissive stop + # that modifies the rays (e.g. a transmission grating or Fresnel + # zone plate) does not reverse them, so its own diffraction or + # refraction must be included in the solve, and the solved + # variable is the pre-stop ray. + if self._stop_is_involutory(surface_first): + propagators = subsystem[1:] + else: + propagators = subsystem + # The residual of the root-finding problem has the same units as # the target grid, so the convergence tolerance must scale with # the size of the target aperture to be achievable in floating @@ -408,7 +442,8 @@ def zfunc(xy: na.AbstractCartesian2dVectorArray): function = functools.partial( self._ray_error, rays=result.outputs, - subsystem=subsystem, + propagators=propagators, + transformation_last=surface_last.transformation, grid_last=grid_last, component_variable=component_variable, component_target=component_target, @@ -477,6 +512,26 @@ def _calc_rayfunction_stops( subsystem = surfaces[index_stop::-1] + # This reversed trace works for a mirror stop because re-applying the + # stop surface reflects the rays back toward the object, and + # reflection is an involution (it exactly undoes itself on the second + # application). A transmissive stop (e.g. a transmission grating or a + # Fresnel zone plate) is not an involution: re-applying it would + # diffract/refract the rays a second time instead of undoing the + # first pass. (A true time-reversed trace through transmissive + # rulings would require negating the diffraction order, while + # reflective rulings are time-reversal symmetric with the same + # order.) So for a non-mirror stop, keep its geometry but strip its + # optical effect; `_calc_rayfunction_stops_only` solves for the + # pre-stop rays in this case. + surface_stop = subsystem[0] + if not self._stop_is_involutory(surface_stop): + subsystem[0] = dataclasses.replace( + surface_stop, + material=None, + rulings=None, + ) + rays_stop = self._calc_rayfunction_stops_only( wavelength_input=wavelength_input, axis_pupil_stop=axis_pupil_stop, diff --git a/optika/systems/_sequential_test.py b/optika/systems/_sequential_test.py index 2143794a..ea5607c8 100644 --- a/optika/systems/_sequential_test.py +++ b/optika/systems/_sequential_test.py @@ -294,6 +294,25 @@ def test__anchor_surface(): assert anchor([first, flat, last]) is last +def test__stop_is_involutory(): + involutory = optika.systems.SequentialSystem._stop_is_involutory + assert involutory( + optika.surfaces.Surface(material=optika.materials.Mirror()), + ) + assert not involutory( + optika.surfaces.Surface( + rulings=optika.rulings.Rulings(spacing=1 * u.um, diffraction_order=1), + ), + ) + assert not involutory( + optika.surfaces.Surface(material=optika.materials.MultilayerFilm()), + ) + assert involutory( + optika.surfaces.Surface(material=optika.materials.Vacuum()), + ) + assert involutory(optika.surfaces.Surface()) + + # small enough that the image of the field fits on the sensor _radius_field_newtonian = 0.05 * u.deg @@ -430,3 +449,69 @@ def test_field_max_matches_source_aperture( result = a.field_max assert np.abs(result.x - _radius_field_grazing) < 1e-6 * u.deg assert np.abs(result.y - _radius_field_grazing) < 1e-6 * u.deg + + +_focal_length_fzp = 100 * u.mm +_wavelength_fzp = 500 * u.nm + +_sensor_fzp = optika.sensors.ImagingSensor( + name="sensor", + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + timedelta_exposure=1 * u.s, + num_pixel=na.Cartesian2dVectorArray(64, 64), + transformation=na.transformations.Cartesian3dTranslation( + z=_focal_length_fzp, + ), + is_field_stop=True, +) + +_system_fzp = optika.systems.SequentialSystem( + object=optika.surfaces.Surface( + name="source", + aperture=optika.apertures.CircularAperture( + radius=np.sin(0.3 * u.deg), + ), + ), + surfaces=[ + optika.surfaces.Surface( + name="fzp", + rulings=optika.rulings.Rulings( + spacing=optika.rulings.HolographicRulingSpacing( + x1=na.Cartesian3dVectorArray(0, 0, -1) * u.AU, + x2=na.Cartesian3dVectorArray(0, 0, 1) * _focal_length_fzp, + wavelength=_wavelength_fzp, + ), + diffraction_order=1, + ), + aperture=optika.apertures.CircularAperture(radius=5 * u.mm), + is_pupil_stop=True, + ), + ], + sensor=_sensor_fzp, + grid_input=_grid_input, +) + + +@pytest.mark.parametrize(argnames="a", argvalues=[_system_fzp]) +class TestSequentialSystemFZP( + AbstractTestAbstractSequentialSystem, +): + """ + A telescope whose only optical element is a transmissive Fresnel zone + plate acting as the pupil stop. This guards against regressions in the + handling of non-involutory (transmissive) stop surfaces by the stop + root-finding problem, which previously re-applied the diffraction of the + stop surface and computed wildly incorrect field extents. + """ + + def test_field_max_matches_plate_scale( + self, + a: optika.systems.AbstractSequentialSystem, + ): + sensor = a.sensor + half_width = sensor.width_pixel * sensor.num_pixel / 2 + field_expected = np.arctan2(half_width, _focal_length_fzp) + result = a.field_max + assert np.abs(result.x - field_expected.x) < 1e-3 * u.deg + assert np.abs(result.y - field_expected.y) < 1e-3 * u.deg