diff --git a/docs/index.rst b/docs/index.rst index 49caa701..fb5fbb32 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -86,6 +86,17 @@ Jupyter notebook examples on how to use :mod:`optika`. tutorials/prime_focus +Guides +====== + +In-depth explanations of how :mod:`optika` works. + +.. toctree:: + :maxdepth: 1 + + stop_finding + + API Reference ============= diff --git a/docs/stop_finding.rst b/docs/stop_finding.rst new file mode 100644 index 00000000..8450744e --- /dev/null +++ b/docs/stop_finding.rst @@ -0,0 +1,315 @@ +Stop Finding +============ + +One of the design goals of :mod:`optika` is to automatically compute the +field of view and entrance pupil of an optical system (see +:attr:`~optika.systems.AbstractSequentialSystem.field_min`, +:attr:`~optika.systems.AbstractSequentialSystem.field_max`, +:attr:`~optika.systems.AbstractSequentialSystem.pupil_min`, and +:attr:`~optika.systems.AbstractSequentialSystem.pupil_max`). +Unlike Zemax, the user never has to specify the extent of the field or pupil. +This page describes the strategy :mod:`optika` uses to discover that extent and +to sample the field and pupil for imaging, since the strategy spans several +methods and the reasoning behind it is easy to lose. + + +Overview +-------- + +The trajectory of a ray through a sequential system is fixed by two surfaces: + +* the **field stop**, which limits the region of the object that is imaged, and +* the **pupil stop** (or aperture stop), which limits the bundle of rays + accepted from each point of the object. + +A surface is marked as a stop by setting ``is_field_stop=True`` or +``is_pupil_stop=True`` on it. Every ray that survives the system passes inside +both stops, so a ray that grazes the *border* of one stop while passing through +a chosen point of the other traces out the boundary of the accepted light. This +is what :attr:`~optika.systems.AbstractSequentialSystem.rayfunction_stops` +returns: a :class:`~optika.rays.RayFunctionArray`, defined on the first surface +of the system, whose rays are constructed to strike prescribed points on both +the field stop and the pupil stop. The field of view, the entrance pupil, and +the sampling of rays used for image simulation are all derived from it. + +The difficulty is that a ray is launched from the *first* surface, but the +constraints live on the *stop* surfaces further downstream. There is no +closed-form expression for the launch coordinate that lands a ray on a given +point of a given stop, so the launch coordinate is found by root-finding. + + +Two principles +-------------- + +Everything below follows from two physical requirements. They are worth stating +up front because they dictate *where* each coordinate is measured, and getting +the measurement plane wrong produces answers that look plausible but are subtly +biased. + +**The field is anchored on the object plane.** + A point-spread function is the image of a single point of the object. For + the simulated PSF to be correct, every ray in a given field bundle must + share one object direction (for a distant object) or one object position + (for a nearby object). If the field were instead anchored on an internal + surface, rays with the same field label but different pupil labels would + correspond to slightly different object points, smearing the PSF. So the + field coordinate is always resolved on the object plane. + +**The pupil is measured at the entrance pupil.** + The entrance pupil is the image of the pupil stop in object space, i.e. the + plane on which the incoming wavefront is uniform for a uniform distant + source. Measuring and sampling the pupil *there*, rather than on the + angular object plane or on the pupil stop itself, matters for two reasons: + + * **Robustness.** For a system with a tiny entrance aperture (a feed optic + much smaller than the beam, as in FURST), the entrance-pupil extent is + essentially independent of field, so rays aimed through it always land on + the aperture. Measured instead as an angle on the object plane, the same + extent swings with field and a single global box makes most field angles + miss the aperture. + + * **Radiometry.** Vignetting is the fraction of the accepted bundle that + survives to the detector, an area integral over the pupil. Sampling + uniformly on the entrance pupil is equal-area sampling in the plane where + the wavefront is uniform, so the vignetting fraction is an unweighted mean + of the surviving samples. Sampling uniformly on the pupil stop instead + would, for a system with pupil distortion (such as the grating in ESIS, + which is the pupil stop), place unequal areas of the entrance pupil under + each sample and bias the result. + +The consequence is that the field is resolved as an object-space coordinate and +the pupil as an entrance-pupil coordinate, and the machinery below exists to +connect those object-space coordinates to the physical stop surfaces. + + +Input coordinates +----------------- + +A stop ray is labelled by three input coordinates, gathered together in an +:class:`~optika.vectors.ObjectVectorArray`: + +``wavelength`` + The vacuum wavelength of the ray. It is carried through the solve + unchanged; it only matters because dispersive surfaces (gratings) bend + different wavelengths differently. + +``field`` + A two-dimensional coordinate locating the ray on the field stop. + +``pupil`` + A two-dimensional coordinate locating the ray on the pupil stop. + +Both ``field`` and ``pupil`` may be given in either **normalized** or +**physical** units, and the units alone tell :mod:`optika` how to interpret +them: + +.. list-table:: + :header-rows: 1 + :widths: 25 25 50 + + * - Units + - Meaning + - Interpretation + * - dimensionless + - normalized + - A value in :math:`[-1, 1]`, denormalized against the bounding box of + the relevant stop's aperture. + * - length (e.g. ``mm``) + - physical position + - A position on the relevant surface. + * - angle (e.g. ``deg``) + - physical direction + - A direction, converted to direction cosines with + :func:`~optika.direction`. + +This convention (dimensionless means normalized, length means position, angle +means direction) is used consistently throughout the stop-finding code; see +``_coordinates_are_normalized`` and +:attr:`~optika.systems.AbstractSequentialSystem.object_is_at_infinity`. + + +The two-point ray solve +----------------------- + +The core primitive is the private method ``_shoot_rays``. Given a +``subsystem`` (a contiguous slice of the system's surfaces), a wavelength, a +grid ``grid_first`` on the first surface, and a grid ``grid_last`` on the last +surface, it finds the launch ray that connects them. Both grids are given in +*physical* units; denormalizing a normalized stop coordinate against its +aperture is the caller's job, so at this boundary a dimensionless grid is an +unambiguous direction cosine rather than a normalized coordinate. + +A ray leaving the first surface has two degrees of freedom that are *not* +pinned by ``grid_first``. If ``grid_first`` is a **position**, the free degrees +of freedom are the launch **direction**; if ``grid_first`` is a **direction**, +they are the launch **position**. ``_shoot_rays`` reads which case applies from +the units of ``grid_first`` (length means position, angle or dimensionless +direction cosine means direction), builds a +:class:`~optika.rays.RayVectorArray` with the fixed coordinate filled in, and +solves for the free coordinate so that the ray lands on ``grid_last`` at the +last surface. The residual whose root is sought is the miss distance at the +last surface, + +.. math:: + + \vec{r}(\vec{a}) = \vec{g}_\text{trial}(\vec{a}) - \vec{g}_\text{last}, + +where :math:`\vec{a}` is the trial value of the free coordinate and +:math:`\vec{g}_\text{trial}` is where the resulting ray actually crosses the +target coordinate of the last surface (its position if ``grid_last`` is a +position, its direction if ``grid_last`` is a direction). This residual is +evaluated by ``_ray_error`` and driven to zero with +:func:`named_arrays.optimize.root_newton`. + +Because ``grid_last`` may itself be either a position or a direction, the same +routine works in either direction along the system. In particular, running it +with the subsystem reversed and an angular ``grid_last`` back-traces rays to a +target *direction* on an object at infinity, which is how the field extent on +the object plane is recovered. + +Two details make the solve robust across systems of wildly different physical +scale: + +* **Seeding.** The initial guess aims each ray from the first surface toward a + sensible target: directly at its point on the last surface when no surface + with optical power lies in between (the guess is then nearly exact), and + otherwise at the center of the first powered surface (a mirror, a curved sag, + or a ruled surface, found by ``_anchor_surface``). This keeps the guess + inside the basin of convergence even for strongly off-axis feed or fold + mirrors. + +* **Scaling.** Both the convergence tolerance and the finite-difference step + used to estimate the Jacobian are scaled by the size of the target aperture, + so that a millimeter-scale spectrograph and a meter-scale telescope are + solved to the same *relative* precision. The default absolute step of + :func:`named_arrays.jacobian` is otherwise below the floating-point noise + floor of the raytrace and yields a Jacobian made of noise. + +``_shoot_rays`` returns the launch rays *at the first surface* (in global +coordinates), with both the given and the solved coordinate filled in, so that +propagating them through the subsystem reproduces ``grid_last``. + + +The strategy +------------ + +The two principles rule out computing a single global field box and a single +global pupil box, because the entrance-pupil extent depends on field. Instead +the field of view and the per-field entrance pupil are calibrated in stages, +and the expensive root-finding is confined to a coarse grid. + +**1. Field extent on the object plane.** + A small number of rays connecting the field-stop border to the pupil-stop + border are traced back to the object plane, and their extent gives the + field of view. This is the object-plane field domain that every later stage + samples within. + +**2. Per-field entrance-pupil extent.** + The object surface is treated as the field stop, and a *coarse* grid of + field directions, sampled evenly across the field of view from stage 1, is + connected to every point on the *wire* (border) of the pupil stop with + ``_shoot_rays``. Only the wire is needed: it is the boundary of the pupil + stop, so its image on the entrance pupil is the boundary of the entrance + pupil, and the minimum and maximum of those positions give a per-field + entrance-pupil bounding box. This is the stage that captures the + field dependence of the pupil (the pupil distortion), and it is the only + stage that pays for dense root-finding, kept affordable by using a coarse + field grid and only the one-dimensional pupil-stop wire. + +**3. Interpolation.** + The per-field bounding box from stage 2 is a smooth function of field (it is + the pupil distortion), so it is fit on the coarse field grid and + interpolated onto the dense field grid, in the same spirit as + :class:`~optika.distortion.PolynomialDistortionModel` and the vignetting + models. The coarse grid must be fine enough to resolve the distortion, not + the scene. + +**4. Dense forward trace.** + For image simulation, each dense ray is fully specified in object space: its + direction is an object-plane field angle, and its position is a normalized + pupil coordinate mapped onto the interpolated per-field entrance-pupil box. + Both coordinates are known without any root-finding, so the dense pass is a + pure *forward* trace through the system. This is what keeps a + :math:`1000 \times 1000` field affordable: the Newton solves live entirely + in the coarse calibration of stage 2, never in the dense grid. + +Because a uniform grid on the entrance-pupil box is equal-area in the plane +where the wavefront is uniform, the vignetted fraction of that grid (with the +``where`` keyword marking the surviving rays) is the vignetting directly, with +no Jacobian weight. The axis-aligned bounding box slightly over-covers a +rotated or astigmatic entrance pupil, but the excess samples fall outside the +aperture and are removed by vignetting, so the result is correct if marginally +less sample-efficient. + + +The object as the field stop +---------------------------- + +Stage 2 relies on being able to mark the **object surface itself** as the field +stop, with an angular (dimensionless, sine-of-half-angle) aperture. The +``field`` coordinate is then a direction on the object plane rather than a +position on an internal surface, which is exactly the object-plane anchoring the +first principle requires. It is also the only workable option in two common +cases where an internal field stop has no solution: + +* **Spectrographs**, where the field stop is the detector. A single wavelength + illuminates only part of the detector, so connecting the *border* of the + detector to the pupil stop has no solution at that wavelength. + +* **Systems with a tiny entrance aperture**, where a single global pupil box, + back-projected to the object, misses the feed optic for most field angles. + +When ``_shoot_rays`` cannot connect an internal field stop at a single +wavelength, the raised error points the user toward this option. Whether the +object is treated as being at infinity is inferred from the units of its +aperture: a length aperture is a finite object, a dimensionless aperture is an +object at infinity (see +:attr:`~optika.systems.AbstractSequentialSystem.object_is_at_infinity`). + + +Field of view and entrance pupil +-------------------------------- + +The field of view and entrance pupil are read off from the stop rayfunction by +reducing over *both* the field and pupil axes: + +* :attr:`~optika.systems.AbstractSequentialSystem.field_min` / + :attr:`~optika.systems.AbstractSequentialSystem.field_max` give the corners + of the field of view, expressed as angles (via :func:`~optika.angles`) when + the object is at infinity and as positions otherwise. + +* :attr:`~optika.systems.AbstractSequentialSystem.pupil_min` / + :attr:`~optika.systems.AbstractSequentialSystem.pupil_max` give the corners + of the entrance pupil, expressed as positions when the object is at infinity + and as angles otherwise. + +The two swap roles with object distance because field and pupil are conjugate: +for a distant object the field is naturally angular and the pupil is a physical +aperture, while for a nearby object the field is a physical extent and the +pupil subtends an angle. + + +Conventions and assumptions +--------------------------- + +A few conventions are load-bearing and worth stating explicitly: + +* **Dimensionless means normalized at the system input, and a direction cosine + inside the solve.** A ``field`` or ``pupil`` grid supplied to the system is + interpreted as normalized when it is dimensionless (see the table above). + Denormalization happens before ``_shoot_rays`` is called, so once inside the + solve a dimensionless grid is unambiguously a physical direction cosine. + +* **The pupil is measured in an object-space plane.** The entrance-pupil extent + and the sampling grid live on the entrance pupil, not on the pupil stop. + Measuring on the pupil stop would reintroduce the distortion bias the second + principle exists to avoid. + +* **Grids are in the local frame** of the surface they belong to. Aperture + denormalization produces coordinates in the surface's own frame, and + ``_shoot_rays`` returns launch rays in the global frame after applying the + first surface's transformation. + +* **Angles are direction cosines.** :func:`~optika.direction` and + :func:`~optika.angles` convert between a pair of azimuth/elevation angles and + a three-dimensional direction cosine, and are inverses of each other. diff --git a/optika/_util.py b/optika/_util.py index e2dde757..2066bf30 100644 --- a/optika/_util.py +++ b/optika/_util.py @@ -74,21 +74,29 @@ def direction( def angles( - direction: na.AbstractCartesian3dVectorArray, + direction: na.AbstractCartesian2dVectorArray | na.AbstractCartesian3dVectorArray, ) -> na.Cartesian2dVectorArray: - """ - Convert a 3D vector of direction cosines to a 2D vector of azimuth and + r""" + Convert a vector of direction cosines to a 2D vector of azimuth and elevation angles. Parameters ---------- direction A vector of direction cosines. + If a 2D vector is given, the :math:`z` component is reconstructed from + the unit-length constraint, :math:`d_z = \sqrt{1 - d_x^2 - d_y^2}`. See Also -------- :func:`direction` : Inverse of this function """ + if not isinstance(direction, na.AbstractCartesian3dVectorArray): + direction = na.Cartesian3dVectorArray( + x=direction.x, + y=direction.y, + z=np.sqrt(1 - np.square(direction.length)), + ) if na.unit(direction) is None: direction = direction << u.dimensionless_unscaled return na.Cartesian2dVectorArray( diff --git a/optika/_util_test.py b/optika/_util_test.py index 1abb9bc9..7480da5d 100644 --- a/optika/_util_test.py +++ b/optika/_util_test.py @@ -36,6 +36,11 @@ def test_direction(angles: na.AbstractCartesian2dVectorArray): argnames="direction", argvalues=[ na.Cartesian3dVectorArray(1, 2, 5).normalized, + na.Cartesian3dVectorArray( + x=na.linspace(-0.4, 0.4, axis="x", num=5), + y=na.linspace(-0.4, 0.4, axis="y", num=5), + z=1, + ).normalized, ], ) def test_angles(direction: na.AbstractCartesian3dVectorArray): @@ -46,3 +51,9 @@ def test_angles(direction: na.AbstractCartesian3dVectorArray): assert isinstance(result, na.AbstractCartesian2dVectorArray) assert np.allclose(direction, optika.direction(result)) + + # a 2D vector of direction cosines reconstructs the same (positive) z from + # the unit-length constraint, so it yields the same angles as the 3D vector + result_from_cosines = optika.angles(direction.xy) + assert isinstance(result_from_cosines, na.AbstractCartesian2dVectorArray) + assert np.allclose(result_from_cosines, result) diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index 71f36074..0e6429fd 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -129,6 +129,30 @@ def wire(self, num: None | int = None) -> na.AbstractCartesian3dVectorArray: wire. """ + def denormalize( + self, + position: na.AbstractCartesian2dVectorArray, + ) -> na.Cartesian2dVectorArray: + """ + Convert a point in normalized (-1, 1) coordinates on the plane of the + aperture to physical coordinates. + + Parameters + ---------- + position + Points in normalized coordinates on the plane of the aperture. + """ + lower = self.bound_lower.xy + upper = self.bound_upper.xy + + ptp = upper - lower + + x = (position + 1) / 2 + + result = ptp * x + lower + + return result + def plot( self, ax: None | matplotlib.axes.Axes | na.ScalarArray[npt.NDArray] = None, diff --git a/optika/apertures/_apertures_test.py b/optika/apertures/_apertures_test.py index bdf01db7..02b2a5c2 100644 --- a/optika/apertures/_apertures_test.py +++ b/optika/apertures/_apertures_test.py @@ -131,6 +131,29 @@ def test_wire(self, a: optika.apertures.AbstractAperture): assert "wire" in wire.shape assert wire.shape["wire"] == a.samples_wire + @pytest.mark.parametrize( + argnames="position", + argvalues=[ + na.Cartesian2dVectorArray(0, 0), + na.Cartesian2dVectorArray( + x=na.linspace(-1, 1, axis="x", num=5), + y=na.linspace(-1, 1, axis="y", num=5), + ), + ], + ) + def test_denormalize( + self, + a: optika.apertures.AbstractAperture, + position: na.AbstractCartesian2dVectorArray, + ): + result = a.denormalize(position) + + assert isinstance(result, na.AbstractCartesian2dVectorArray) + + unit = na.unit_normalized(a.wire()) + + assert na.unit_normalized(result).is_equivalent(unit) + class TestPlot( test_mixins.AbstractTestPlottable.TestPlot, ): diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 8e2eed3f..72d502e2 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1,4 +1,5 @@ from __future__ import annotations + from typing import Sequence, Callable, Any, ClassVar import abc import dataclasses @@ -101,11 +102,24 @@ def surfaces_all(self) -> list[optika.surfaces.AbstractSurface]: else: result = [] - result += list(self.surfaces) + surfaces = list(self.surfaces) sensor = self.sensor if sensor is not None: - result += [sensor] + surfaces += [sensor] + + transformation = self.transformation + + result += [ + dataclasses.replace( + surface, + transformation=na.transformations.compose( + transformation, + surface.transformation, + ), + ) + for surface in surfaces + ] if not any(s.is_field_stop for s in result): result[0] = dataclasses.replace(result[0], is_field_stop=True) @@ -376,6 +390,10 @@ def _ray_error( rays_component_variable.y = a.y rays_component_variable.z = zfunc(a) + surface_first = subsystem[0] + if surface_first.transformation is not None: + rays = surface_first.transformation(rays) + rays = optika.propagators.propagate_rays( propagators=subsystem[1:], rays=rays, @@ -393,6 +411,230 @@ def _ray_error( result = grid_last_trial - grid_last return result + def _shoot_rays( + self, + subsystem: Sequence[optika.surfaces.AbstractSurface], + wavelength: u.Quantity | na.AbstractScalar, + grid_first: na.AbstractCartesian2dVectorArray, + grid_last: na.AbstractCartesian2dVectorArray, + ) -> optika.rays.RayVectorArray: + """ + Find the position or direction needed for rays launched from `grid_first` + on the first surface of `subsystem` to land at `grid_last` on the + last surface of `subsystem`. + + Both grids must be given in physical units; denormalizing normalized + stop coordinates against a surface's aperture is the responsibility of + the caller. + + Parameters + ---------- + subsystem + A sequence of surfaces to trace rays through. + wavelength + The vacuum wavelength of the rays. + grid_first + The physical position or direction of the rays launched from the + first surface of `subsystem`. + If in length units, `grid_first` represents positions on the first + surface, and this function solves for the starting direction of the + rays. + If in angular units, or dimensionless direction cosines, + `grid_first` represents directions on the first surface, and this + function solves for the starting position of the rays. + grid_last + The target physical position or direction of the rays incident on + the last surface of `subsystem`. + If in length units, `grid_last` represents target positions on the + last surface. + If in angular units, or dimensionless direction cosines, `grid_last` + represents target directions on the last surface. + + Returns + ------- + The rays launched from the first surface of `subsystem`, with both the + given `grid_first` coordinate and the solved-for complementary + coordinate filled in, so that propagating them through `subsystem` + lands them at `grid_last` on the last surface. + These rays are in global coordinates. + """ + surface_first = subsystem[0] + surface_last = subsystem[~0] + + # Fold angular grids into direction cosines, so that below the solve + # only needs to distinguish positions (length units) from directions + # (dimensionless direction cosines). + if na.unit_normalized(grid_first).is_equivalent(u.deg): + grid_first = optika.direction(grid_first).xy + if na.unit_normalized(grid_last).is_equivalent(u.deg): + grid_last = optika.direction(grid_last).xy + + rays = optika.rays.RayVectorArray(wavelength=wavelength) + + if na.unit_normalized(grid_first).is_equivalent(u.mm): + # `grid_first` is a position on the first surface, so solve for the + # starting direction of the rays. + grid_first = na.Cartesian3dVectorArray( + x=grid_first.x, + y=grid_first.y, + z=0 * u.mm, + ) + rays.position = grid_first.replace(z=surface_first.sag(grid_first)) + rays.direction = na.Cartesian3dVectorArray(0, 0, 1) + component_variable = "direction" + + def zfunc(xy: na.AbstractCartesian2dVectorArray): + return np.sqrt(1 - np.square(xy.length)) + + else: + # `grid_first` is a direction cosine on the first surface, so solve + # for the starting position of the rays. + rays.direction = na.Cartesian3dVectorArray( + x=grid_first.x, + y=grid_first.y, + z=np.sqrt(1 - np.square(grid_first.length)), + ) + rays.position = na.Cartesian3dVectorArray() * u.mm + component_variable = "position" + + 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) + + # Seed the free ray component by aiming each ray at its own target + # point on the last surface when no surface with optical power lies + # between the two ends (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 surface (e.g. an off-axis fold or feed mirror). + anchor = self._anchor_surface(subsystem) + if anchor is surface_last and na.unit_normalized(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 - rays.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 + rays.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 = rays.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) + rays.position = position_seed + + # if surface_first.transformation is not None: + # rays = surface_first.transformation(rays) + + if na.unit_normalized(grid_last).is_equivalent(u.mm): + component_target = "position" + else: + component_target = "direction" + + # 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 point for systems of + # any physical scale. + scale = np.maximum( + grid_last.x.ptp(), + grid_last.y.ptp(), + ) + max_abs_error = 1e-9 * np.maximum( + scale, + 1 * na.unit_normalized(scale), + ) + + variables = getattr(rays, component_variable) + + # The internal seed only carries the axes it was built from (the "rough" + # seed, aimed at a powered surface, has the field axes but no target-grid + # axis), so broadcast it to the full solution shape. Otherwise the solve + # is under-determined (one unknown per seed axis, but a residual over the + # combined grid), which cannot converge. + guess = na.Cartesian2dVectorArray(x=variables.x, y=variables.y) + guess = guess.broadcast_to(na.shape_broadcasted(guess, grid_last)) + + function = functools.partial( + self._ray_error, + rays=rays, + subsystem=subsystem, + grid_last=grid_last, + component_variable=component_variable, + component_target=component_target, + zfunc=zfunc, + ) + + # The default perturbation used by `na.jacobian` is an absolute 1e-10 + # in the units of the variable, which is below the floating-point noise + # floor of the raytrace for variables measured in physical units, + # yielding a Jacobian made of noise. Use a perturbation proportional to + # the scale of the problem instead. + 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) + + try: + root = na.optimize.root_newton( + function=function, + guess=guess, + jacobian=jacobian, + max_abs_error=max_abs_error, + ) + except ValueError as e: # pragma: nocover + raise ValueError( + f"Could not solve for the rays connecting surfaces " + f"{surface_first.name!r} and {surface_last.name!r}. " + f"If the target surface is only partially reachable at a " + f"single wavelength (e.g. a spectrograph sensor), consider " + f"marking the object surface, with an angular (dimensionless, " + f"sine of the half-angle) aperture, as the field stop instead." + ) from e + + variables.x = root.x + variables.y = root.y + variables.z = zfunc(root) + + if surface_first.transformation is not None: + rays = surface_first.transformation(rays) + + return rays + def _calc_rayfunction_stops_only( self, wavelength_input: na.ScalarLike, @@ -458,168 +700,13 @@ def _calc_rayfunction_stops_only( result.inputs.field = grid_first result.inputs.pupil = grid_last - if na.unit(grid_first).is_equivalent(u.mm): - grid_first = na.Cartesian3dVectorArray( - x=grid_first.x, - y=grid_first.y, - z=0 * u.mm, - ) - result.outputs.position = grid_first.replace( - z=surface_first.sag(grid_first) - ) - result.outputs.direction = na.Cartesian3dVectorArray(0, 0, 1) - component_variable = "direction" - - def zfunc(xy: na.AbstractCartesian2dVectorArray): - return np.sqrt(1 - np.square(xy.length)) - - elif na.unit(grid_first).is_equivalent(u.dimensionless_unscaled): - result.outputs.direction = na.Cartesian3dVectorArray( - x=grid_first.x, - y=grid_first.y, - z=np.sqrt(1 - np.square(grid_first.length)), - ) - result.outputs.position = na.Cartesian3dVectorArray() * u.mm - component_variable = "position" - - 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) - - if na.unit(grid_last).is_equivalent(u.mm): - component_target = "position" - elif na.unit(grid_last).is_equivalent(u.dimensionless_unscaled): - component_target = "direction" - else: - raise ValueError(f"unrecognized output grid unit, {na.unit(grid_last)}") - - # 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 - # point for systems of any physical scale. - scale = np.maximum( - grid_last.x.ptp(), - grid_last.y.ptp(), - ) - max_abs_error = 1e-9 * np.maximum( - scale, - 1 * na.unit_normalized(scale), - ) - - variables = getattr(result.outputs, component_variable) - - function = functools.partial( - self._ray_error, - rays=result.outputs, + result.outputs = self._shoot_rays( subsystem=subsystem, + wavelength=wavelength_input, + grid_first=grid_first, grid_last=grid_last, - component_variable=component_variable, - component_target=component_target, - zfunc=zfunc, ) - # The default perturbation used by `na.jacobian` is an absolute - # 1e-10 in the units of the variable, which is below the - # floating-point noise floor of the raytrace for variables - # measured in physical units, yielding a Jacobian made of noise. - # Use a perturbation proportional to the scale of the problem - # instead. - 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) - - try: - root = na.optimize.root_newton( - function=function, - guess=na.Cartesian2dVectorArray( - x=variables.x, - y=variables.y, - ), - jacobian=jacobian, - max_abs_error=max_abs_error, - ) - except ValueError as e: # pragma: nocover - raise ValueError( - f"Could not solve for the rays connecting the stop " - f"surfaces {surface_first.name!r} and " - f"{surface_last.name!r}. " - f"If the field stop is only partially reachable at a " - f"single wavelength (e.g. a spectrograph sensor), " - f"consider marking the object surface, with an angular " - f"(dimensionless, sine of the half-angle) aperture, as " - f"the field stop instead." - ) from e - - variables.x = root.x - variables.y = root.y - variables.z = zfunc(root) - return result def _calc_rayfunction_stops( @@ -627,8 +714,8 @@ def _calc_rayfunction_stops( wavelength_input: na.ScalarLike, axis_pupil_stop: str, axis_field_stop: str, - samples_pupil_stop: int = 101, - samples_field_stop: int = 101, + samples_pupil_stop: int = 21, + samples_field_stop: int = 21, ) -> optika.rays.RayFunctionArray: surfaces = self.surfaces_all @@ -659,7 +746,8 @@ def _calc_rayfunction_stops( rays = obj.transformation.inverse(rays) where = rays.direction @ obj.sag.normal(rays.position) > 0 - result.outputs.direction[where] = -result.outputs.direction[where] + rays.direction[where] = -rays.direction[where] + result.outputs = rays # If the first stop is the object surface, the solved variable is the # position and the direction retains only the field-stop axis, so @@ -672,14 +760,18 @@ def _calc_rayfunction_stops( 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) - return result _axis_pupil_stop: ClassVar[str] = "_stop_pupil" _axis_field_stop: ClassVar[str] = "_stop_field" + @property + def _axis_stops(self) -> tuple[str, str]: + """ + Tuple of :attr:`_axis_pupil_stop` and :attr:`_axis_field_stop`. + """ + return (self._axis_field_stop, self._axis_pupil_stop) + @functools.cached_property def rayfunction_stops(self) -> optika.rays.RayFunctionArray: """ @@ -695,55 +787,331 @@ def rayfunction_stops(self) -> optika.rays.RayFunctionArray: samples_field_stop=21, ) + def _fit_vs_wavelength( + self, + a: na.AbstractCartesian2dVectorArray, + ) -> na.PolynomialFitFunctionArray: + """ + Fit an array of values, as a function of wavelength, + using at most a second-order polynomial. + + Parameters + ---------- + a + The array of values to fit. + """ + rays = self.rayfunction_stops.outputs + + axis_wavelength = self.axis_wavelength_ + if axis_wavelength: + (axis,) = axis_wavelength + degree = int(np.minimum(2, rays.wavelength.shape[axis] - 1)) + else: + axis = None + degree = 0 + + return na.PolynomialFitFunctionArray.from_degree( + inputs=rays.wavelength, + outputs=a, + degree=degree, + axis_polynomial=axis, + ) + + @property + def _field_min(self) -> na.PolynomialFitFunctionArray: + rays = self.rayfunction_stops.outputs + if self.object_is_at_infinity: + a = rays.direction.xy.min(self._axis_stops) + else: + a = rays.position.xy.min(self._axis_stops) + return self._fit_vs_wavelength(a) + + @property + def _field_max(self) -> na.PolynomialFitFunctionArray: + rays = self.rayfunction_stops.outputs + if self.object_is_at_infinity: + a = rays.direction.xy.max(self._axis_stops) + else: + a = rays.position.xy.max(self._axis_stops) + return self._fit_vs_wavelength(a) + @property def field_min(self) -> na.AbstractCartesian2dVectorArray: """ The lower left corner of this optical system's field of view. """ - axis = (self._axis_field_stop, self._axis_pupil_stop) + wavelength = self.grid_input.wavelength if self.object_is_at_infinity: - angles = optika.angles(self.rayfunction_stops.outputs.direction) - return angles.min(axis) - else: - return self.rayfunction_stops.outputs.position.xy.min(axis) + # `optika.angles` reverses the sign of each direction cosine, so the + # minimum-angle corner of the field of view is the maximum-cosine + # corner: evaluate the `_field_max` fit here (and `_field_min` in + # `field_max`). + return optika.angles(self._field_max(wavelength).outputs) + return self._field_min(wavelength).outputs @property def field_max(self) -> na.AbstractCartesian2dVectorArray: """ The upper right corner of this optical system's field of view. """ - axis = (self._axis_field_stop, self._axis_pupil_stop) + wavelength = self.grid_input.wavelength if self.object_is_at_infinity: - angles = optika.angles(self.rayfunction_stops.outputs.direction) - return angles.max(axis) - else: - return self.rayfunction_stops.outputs.position.xy.max(axis) + return optika.angles(self._field_min(wavelength).outputs) + return self._field_max(wavelength).outputs + + def _denormalize_field( + self, + wavelength: u.Quantity | na.AbstractScalar, + field: na.AbstractCartesian2dVectorArray, + ) -> na.AbstractCartesian2dVectorArray: + """ + Convert normalized field coordinates to physical units. + + If :attr:`object_is_at_infinitiy`, the result is in terms of + direction cosines, otherwise the result is in length units. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + field + An array of normalized 2D field coordinates. + """ + field_min = self._field_min(wavelength).outputs + field_max = self._field_max(wavelength).outputs + field_ptp = field_max - field_min + + result = field_ptp * (field + 1) / 2 + field_min + + return result @property - def pupil_min(self) -> na.AbstractCartesian2dVectorArray: + def field_(self) -> na.AbstractCartesian2dVectorArray: + """ + Normalized version of ``grid_input.field``. + """ + wavelength = self.grid_input.wavelength + field = self.grid_input.field + + if na.unit_normalized(field).is_equivalent(u.dimensionless_unscaled): + + field = self._denormalize_field(wavelength, field) + + if self.object_is_at_infinity: + field = optika.angles(field) + + return field + + def _calc_rayfunction_pupil( + self, + wavelength: u.Quantity | na.AbstractScalar, + field: na.AbstractCartesian2dVectorArray, + # samples_field: int = 5, + samples_pupil: int = 21, + ) -> optika.rays.RayFunctionArray: """ - The lower left corner of this optical system's entrance pupil in - physical units. + Trace a grid of rays from the object plane to the wire of the pupil + stop, in order to measure how the entrance pupil depends on field. + + Unlike :attr:`rayfunction_stops`, which samples the field only on the + *border* of the field stop, this samples the field on a coarse grid + across the interior of the object's aperture (the field of view), so + that the field dependence of the entrance pupil (the pupil distortion) + can be captured. Only the *wire* of the pupil stop is sampled, since its + image on the entrance pupil is the boundary of the entrance pupil, whose + minimum and maximum are all that is needed. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + field + The field coordinates in physical units. + samples_pupil + The number of samples along the wire of the pupil stop. """ - axis = (self._axis_field_stop, self._axis_pupil_stop) + subsystem = self.surfaces_all[: self.index_pupil_stop + 1] + + pupil_stop = subsystem[~0] + + aperture_pupil = pupil_stop.aperture + + # the wire (border) of the pupil stop + wire = np.moveaxis( + a=aperture_pupil.wire(num=samples_pupil), + source="wire", + destination=self._axis_pupil_stop, + ) + pupil = wire.xy + + # solve for the launch coordinate on the object plane (the entrance-pupil + # footprint) connecting each field to each point of the pupil-stop wire + rays = self._shoot_rays( + subsystem=subsystem, + wavelength=wavelength, + grid_first=field, + grid_last=pupil, + ) + + obj = subsystem[0] + if obj.transformation is not None: + rays = obj.transformation.inverse(rays) + + return optika.rays.RayFunctionArray( + inputs=optika.vectors.ObjectVectorArray( + wavelength=wavelength, + field=field, + ), + outputs=rays, + ) + + @functools.cached_property + def _rayfunction_pupil(self) -> optika.rays.RayFunctionArray: + """ + Trace a grid of rays from the object plane to the wire of the pupil + stop, using the default wavelength and field grids. + """ + return self._calc_rayfunction_pupil( + wavelength=self.grid_input.wavelength, + field=self.field_, + ) + + def _fit_vs_wavelength_and_field( + self, + a: na.AbstractCartesian2dVectorArray, + ) -> na.PolynomialFitFunctionArray: + """ + Fit an array of values, as a function of wavelength and field coordinate, + using at most a second-order polynomial. + + Parameters + ---------- + a + The array of values to fit. + """ + rays = self._rayfunction_pupil.outputs + + axis_wavelength = self.axis_wavelength_ + axis_field = self.axis_field_ + + if axis_wavelength: + axis = axis_wavelength + components = ("wavelength",) + else: + axis = () + components = () + + axis = axis + axis_field + components = components + ("field.x", "field.y") + if self.object_is_at_infinity: - return self.rayfunction_stops.outputs.position.xy.min(axis) + field = rays.direction.xy else: - angles = optika.angles(self.rayfunction_stops.outputs.direction) - return angles.min(axis) + field = rays.position.xy + + inputs = optika.vectors.SceneVectorArray(rays.wavelength, field) + + shape = inputs.shape + shape = [shape[ax] for ax in axis] + degree = int(np.minimum(2, min(shape) - 1)) + + return na.PolynomialFitFunctionArray.from_degree( + inputs=inputs, + outputs=a, + degree=degree, + components=components, + axis_polynomial=axis, + ) @property - def pupil_max(self): - """ - The upper right corner of this optical system's entrance pupil in - physical units. - """ - axis = (self._axis_field_stop, self._axis_pupil_stop) + def _pupil_min(self) -> na.PolynomialFitFunctionArray: + rays = self._rayfunction_pupil.outputs + if self.object_is_at_infinity: + a = rays.position.xy.min(self._axis_pupil_stop) + else: + a = rays.direction.xy.min(self._axis_pupil_stop) + return self._fit_vs_wavelength_and_field(a) + + @property + def _pupil_max(self) -> na.PolynomialFitFunctionArray: + rays = self._rayfunction_pupil.outputs if self.object_is_at_infinity: - return self.rayfunction_stops.outputs.position.xy.max(axis) + a = rays.position.xy.max(self._axis_pupil_stop) else: - angles = optika.angles(self.rayfunction_stops.outputs.direction) - return angles.max(axis) + a = rays.direction.xy.max(self._axis_pupil_stop) + return self._fit_vs_wavelength_and_field(a) + + # @property + # def pupil_min(self) -> na.AbstractCartesian2dVectorArray: + # """ + # The lower left corner of this optical system's entrance pupil in + # physical units. + # """ + # axis = (self._axis_field_stop, self._axis_pupil_stop) + # if self.object_is_at_infinity: + # return self.rayfunction_stops.outputs.position.xy.min(axis) + # else: + # angles = optika.angles(self.rayfunction_stops.outputs.direction) + # return angles.min(axis) + # + # @property + # def pupil_max(self): + # """ + # The upper right corner of this optical system's entrance pupil in + # physical units. + # """ + # axis = (self._axis_field_stop, self._axis_pupil_stop) + # if self.object_is_at_infinity: + # return self.rayfunction_stops.outputs.position.xy.max(axis) + # else: + # angles = optika.angles(self.rayfunction_stops.outputs.direction) + # return angles.max(axis) + + def _denormalize_pupil( + self, + wavelength: u.Quantity | na.AbstractScalar, + field: na.AbstractCartesian2dVectorArray, + pupil: na.AbstractCartesian2dVectorArray, + ) -> na.AbstractCartesian2dVectorArray: + """ + Convert normalized field coordinates to physical units. + + If :attr:`object_is_at_infinitiy`, the result is in terms of + direction cosines, otherwise the result is in length units. + + Parameters + ---------- + wavelength + The wavelength of the incident light in vacuum. + field + An array of 2D field coordinates in terms of direction cosines. + """ + x = optika.vectors.SceneVectorArray(wavelength, field) + pupil_min = self._pupil_min(x).outputs + pupil_max = self._pupil_max(x).outputs + pupil_ptp = pupil_max - pupil_min + + result = pupil_ptp * (pupil + 1) / 2 + pupil_min + + return result + + @property + def pupil_(self) -> na.AbstractCartesian2dVectorArray: + """ + Normalized version of ``grid_input.pupil``. + """ + wavelength = self.grid_input.wavelength + field = self.field_ + pupil = self.grid_input.pupil + + if na.unit_normalized(field).is_equivalent(u.dimensionless_unscaled): + + pupil = self._denormalize_pupil(wavelength, field, pupil) + + if not self.object_is_at_infinity: + pupil = optika.angles(pupil) + + return pupil def _denormalize_grid( self, @@ -752,41 +1120,63 @@ def _denormalize_grid( normalized_pupil: bool = True, ) -> optika.vectors.ObjectVectorArray: - if (not normalized_field) and (not normalized_pupil): - return grid - - axis_field = self._axis_field_stop - axis_pupil = self._axis_pupil_stop - - rayfunction_stops = self._calc_rayfunction_stops( - wavelength_input=grid.wavelength, - axis_pupil_stop=axis_pupil, - axis_field_stop=axis_field, - samples_pupil_stop=21, - samples_field_stop=21, - ) - - result = grid.copy_shallow() - - object_is_at_infinity = self.object_is_at_infinity - if object_is_at_infinity: - field = optika.angles(rayfunction_stops.outputs.direction) - pupil = rayfunction_stops.outputs.position.xy - else: - field = rayfunction_stops.outputs.position.xy - pupil = optika.angles(rayfunction_stops.outputs.direction) + grid = grid.copy_shallow() if normalized_field: - min_field = field.min(axis=(axis_field, axis_pupil)) - ptp_field = field.ptp(axis=(axis_field, axis_pupil)) - result.field = ptp_field * (result.field + 1) / 2 + min_field + field = self._denormalize_field( + wavelength=grid.wavelength, + field=grid.field, + ) + grid.field = field + if self.object_is_at_infinity: + grid.field = optika.angles(grid.field) if normalized_pupil: - min_pupil = pupil.min(axis=(axis_field, axis_pupil)) - ptp_pupil = pupil.ptp(axis=(axis_field, axis_pupil)) - result.pupil = ptp_pupil * (result.pupil + 1) / 2 + min_pupil - - return result + grid.pupil = self._denormalize_pupil( + wavelength=grid.wavelength, + field=field, + pupil=grid.pupil, + ) + if not self.object_is_at_infinity: + grid.pupil = optika.angles(grid.pupil) + + return grid + + # if (not normalized_field) and (not normalized_pupil): + # return grid + # + # axis_field = self._axis_field_stop + # axis_pupil = self._axis_pupil_stop + # + # rayfunction_stops = self._calc_rayfunction_stops( + # wavelength_input=grid.wavelength, + # axis_pupil_stop=axis_pupil, + # axis_field_stop=axis_field, + # samples_pupil_stop=21, + # samples_field_stop=21, + # ) + # + # result = grid.copy_shallow() + # + # object_is_at_infinity = self.object_is_at_infinity + # if object_is_at_infinity: + # field = optika.angles(rayfunction_stops.outputs.direction) + # pupil = rayfunction_stops.outputs.position.xy + # else: + # field = rayfunction_stops.outputs.position.xy + # pupil = optika.angles(rayfunction_stops.outputs.direction) + # + # if normalized_field: + # min_field = field.min(axis=(axis_field, axis_pupil)) + # ptp_field = field.ptp(axis=(axis_field, axis_pupil)) + # result.field = ptp_field * (result.field + 1) / 2 + min_field + # + # if normalized_pupil: + # min_pupil = pupil.min(axis=(axis_field, axis_pupil)) + # ptp_pupil = pupil.ptp(axis=(axis_field, axis_pupil)) + # result.pupil = ptp_pupil * (result.pupil + 1) / 2 + min_pupil + # + # return result def _calc_rayfunction_input( self, @@ -905,8 +1295,6 @@ def raytrace( rays = result.outputs if intensity is not None: rays.intensity = intensity - if self.transformation is not None: - rays = self.transformation.inverse(rays) surfaces = self.surfaces_all @@ -980,8 +1368,12 @@ def rayfunction( rayfunction = raytrace[{axis: ~0}] rays = rayfunction.outputs - if self.sensor.transformation is not None: - rays = self.sensor.transformation.inverse(rays) + # `raytrace` returns rays in the object frame, and the sensor sits at + # its composed position there, so invert the composed sensor + # transformation to reach the sensor's local frame. + sensor = self.surfaces_all[~0] + if sensor.transformation is not None: + rays = sensor.transformation.inverse(rays) rayfunction.outputs = rays @@ -1469,15 +1861,8 @@ def plot( """ surfaces = self.surfaces_all - transformation_self = self.transformation kwargs_plot = self.kwargs_plot - if transformation is not None: - if transformation_self is not None: - transformation = transformation @ transformation_self - else: - transformation = transformation_self - if kwargs_plot is not None: kwargs = kwargs | kwargs_plot @@ -1654,12 +2039,6 @@ def _write_to_dxf( **kwargs, ) -> None: - if self.transformation is not None: - if transformation is not None: - transformation = transformation @ self.transformation - else: - transformation = self.transformation - surfaces = self.surfaces_all for surface in surfaces: