SequentialSystem backward raytrace fails for transmissive pupil stops (FZP, gratings, lenses)

[jacobdparker]

Bug description

SequentialSystem.rayfunction() and any method that uses normalized field/pupil
coordinates (field=1, pupil=1, etc.) produces wrong results when the pupil stop
is a transmissive optic — specifically a Fresnel Zone Plate (FZP) or any surface
whose material is not a mirror.

The internal _calc_rayfunction_stops method performs a backward raytrace to
determine the physical extents of the field and pupil stops, which are then used
to map normalized [-1, 1] coordinates to physical units. For transmissive pupil
stops this backward raytrace gives extents that are off by ~20×, causing all
normalized-coordinate calls (raytrace(), rayfunction(), spot_diagram(),
plot(plot_rays=True)) to trace rays far outside the valid aperture.

Purely reflective systems (e.g. a prime-focus mirror telescope) are not affected.

Root cause

There are two coupled errors in _calc_rayfunction_stops_only and
_calc_rayfunction_stops in optika/systems/_sequential.py.

1. Root-finding skips the grating (_ray_error, line ~198)

_calc_rayfunction_stops_only places trial rays at the pupil stop edge and
root-finds the direction that makes them reach the field stop boundary. The inner
propagation call uses subsystem[1:], which omits subsystem[0] — the FZP itself.
The root-finder therefore finds straight-line geometric angles that reach the sensor
without applying diffraction.

2. Backward trace re-applies the grating (_calc_rayfunction_stops, line ~373)

_calc_rayfunction_stops takes the (already-wrong) directions from step 1 and
propagates them backward through the full subsystem = [FZP, source]. Because the
FZP is transmissive, propagate_rays applies the grating equation again in the
forward direction, compounding the error.

For a mirror this step works correctly: Snell's law with is_mirror=True reflects
the ray direction, which is the natural time-reversal of reflection. Transmissive
surfaces have no equivalent automatic reversal.

Workaround

Build rayfunction_default manually with physical (non-normalized) coordinates and
cache it on the system before calling any normalized-coordinate method:

rays_default = system.rayfunction(
    field=field,          # physical units (deg or mm)
    pupil=pupil,          # physical units (mm or deg)
    normalized_field=False,
    normalized_pupil=False,
)
system.__dict__["rayfunction_default"] = rays_default

Minimal reproduction

import numpy as np
import astropy.units as u
import named_arrays as na
import optika

focal_length = 1000 * u.mm
radius = 70 * u.mm
wavelength = 171 * u.AA

spacing = optika.rulings.HolographicRulingSpacing(
    x1=na.Cartesian3dVectorArray(0, 0, -1) * (1 * u.au).to(u.mm),
    x2=na.Cartesian3dVectorArray(0 * u.mm, 0 * u.mm, focal_length),
    wavelength=wavelength,
    is_diverging_1=True,
    is_diverging_2=False,
)
fzp = optika.surfaces.Surface(
    name="FZP",
    aperture=optika.apertures.CircularAperture(radius),
    rulings=optika.rulings.Rulings(spacing=spacing, diffraction_order=1),
    is_pupil_stop=True,
)
sensor = optika.sensors.ImagingSensor(
    name="sensor",
    width_pixel=13 * u.um,
    axis_pixel=na.Cartesian2dVectorArray("x", "y"),
    num_pixel=na.Cartesian2dVectorArray(512, 512),
    transformation=na.transformations.Cartesian3dTranslation(z=focal_length),
    is_field_stop=True,
)
grid = optika.vectors.ObjectVectorArray(
    wavelength=wavelength,
    field=na.Cartesian2dVectorLinearSpace(start=-1, stop=1,
        axis=na.Cartesian2dVectorArray("fx","fy"), num=5),
    pupil=na.Cartesian2dVectorLinearSpace(start=-1, stop=1,
        axis=na.Cartesian2dVectorArray("px","py"), num=5),
)
system = optika.systems.SequentialSystem(
    surfaces=[fzp], sensor=sensor, grid_input=grid,
)

# Bug: field_max should be ~0.191 deg, pupil_max should be ~70 mm
print(system.field_max)   # gives ~3-4 deg (wrong)
print(system.pupil_max)   # gives ~1400 mm (wrong)

Notes on the fix

The fix requires detecting whether the first stop surface is transmissive
(not surface.material.is_mirror) and:

1. Including the surface in the root-finding propagation so diffraction
   is part of the residual function.
2. Excluding the surface from the backward propagation (the pre-diffraction
   incident directions found by root-finding are already the correct object-space
   directions; re-applying the grating would corrupt them).

An alternative, more architecturally principled approach would be to introduce a
reverse=True mode in propagate_rays that negates the diffraction order (m → −m)
for time-reversal. This requires restructuring the backward-trace algorithm and
propagating a reverse flag through propagators.py, surfaces.py, and
rulings/_rulings.py, and is recommended as a future refactor.