Exact walk-off collar sum, the modal exponent-budget reason, and an exact per-band loss recipe

changelog.d/per-band-loss-interval.fixed.md

@@ -0,0 +1 @@
+The boundaries guide reads a band's absorbed energy from saved planes with an exact interval formula in place of the trapezoid, which overcounted a strongly absorbed band's loss by up to twice; the medium's own absorption is now its own share (#573).

changelog.d/walkoff-periodic-sum.fixed.md

@@ -0,0 +1 @@
+The temporal collar's walk-off integral, read by `warn_walkoff` and by `absorber(grid, walkoff=...)`, is now the rectangle sum over the periodic time axis, which is exact for the taper, rather than a trapezoid that left out one sample; a collar sized for a drift now removes exactly the stated `exp(-2Γ)` of its power, and the warning fires at the exact threshold (#490).

docs/dev/kernel-contract.md

@@ -447,7 +447,10 @@ bounds read off the operator (`translate.radial.radial_rate_bounds`): a
 constant table's exact extremes of `Re D`, and, for a varying medium, a sufficient bound
 on each interval of the rate's Chebyshev series (`|T_j| ≤ 1`) for a z table or for an
 expansion's scalar weights against its terms. The refusal names the radial operator's
-attenuation and the admissible step. The modal kinds are not checked.
+attenuation and the admissible step. The modal kinds are not checked: the default
+`LossPolicy` clamps a mode's power attenuation on `D` to `2 × 708 / L`, so no step within
+the run forms an exponent past 708; gain, a raised `max_alpha` and the absorber's rates,
+added after the clamp, are outside that bound.
 
 The imaginary part of Φ is also the unitary phase for the noise field. Modal Φ
 includes the real endpoint term `-1/2 [ln w(z) - ln w(0)]`. It vanishes identically

docs/guide/numerics/boundaries.md

@@ -330,7 +330,8 @@ and nothing else. So:
   long drift through a narrow collar is better served by a wider window.
 
 A linear packet walking through this collar at the stated drift leaves `exp(−2Γ)` of its
-power, a millionth, to within the one sample the sizing's integral leaves out of the period.
+power, a millionth: the sizing sums the rate over the periodic axis, which is exact for this
+taper, and the measured exponent agrees with it to a few parts in `10⁹`.
 The sizing assumes a steady drift over the run, as `warn_walkoff` does. Walk-off that starts
 late, after fission or as a soliton decelerates, crosses faster than the mean and is absorbed
 less, so the [monitor](#watching-the-windows-edges) stays the measurement. The default
@@ -582,12 +583,27 @@ the spatial routes' default, and about half off the in-stage term, on the same g
 light can lose most of its band while the run's energy moves by `10⁻⁵`, and the band's loss is
 what says the window is too short. The rate is diagonal in `ω`, so each bin loses
 `α(ω) S(ω, z)`, with `S` the `beam_spectrum`, integrated in `z`. That statistic is taken at saved
-planes and not [between them](../results/statistics.md#between-the-planes), so the integral runs
-over the planes, with an error of about `(α Δz)²/12` of a bin's loss by the trapezoid. Dividing
-by what the band kept and lost together counts light the run generated part-way as well as
-light it launched. `absorbed_power` has the user's share as its `custom` channel, but over the
-whole beam: only this product reads it per band. A lossy medium's own absorption is a further
-channel: fused silica absorbs at these wavelengths, and what it takes appears in neither count.
+planes and not [between them](../results/statistics.md#between-the-planes), so the integral has
+to be read from the planes.
+
+- **Not by the trapezoid.** Its error, about `(a Δz)²/12` of a bin's loss for a bin falling at
+  the rate `a`, is small only while `a Δz` is. The bins this rate exists for have `a Δz` large:
+  light already in one is gone within a fraction of the interval, and the trapezoid charges it
+  at the plane's value over the whole of it. On a 5 µm slab of fused silica with a 1 fs launch
+  at 260 nm it overcounted a band's loss by 1.5 to 2 times, and converged slowly as planes
+  were added.
+- **Between two planes.** A bin obeys `dS/dz = G − a S`, with `a` every rate that drains it,
+  the absorber's and the medium's, and `G` whatever else feeds it. Taking `G` constant over the
+  interval, the two planes fix it:
+  `G Δz = x (S₁ − S₀ e⁻ˣ)/(1 − e⁻ˣ)` with `x = a Δz`, and the drains removed `G Δz − (S₁ − S₀)`.
+  That is exact for a bin that only decays and for one fed at a steady rate, and each rate takes
+  its share `α/a` of it.
+
+Dividing by what the band kept and lost together counts light the run generated part-way as
+well as light it launched. `absorbed_power` has the user's share as its `custom` channel, but
+over the whole beam: only this product reads it per band. A lossy medium's own absorption is a
+further channel: fused silica absorbs at these wavelengths, so its rate is one of the drains,
+and what it takes is its own share, in neither the walk-off's count nor the collar's.
 
 ```python
 planes = np.linspace(length / 32, length, 32)
@@ -611,14 +627,25 @@ output = Simulation.from_free(
 # (n_z, n_omega), J/(rad/s) per metre of y
 spectrum = np.asarray(output["stats/beam_spectrum"])
 z = np.asarray(output["z"])
-lost = np.trapezoid(walkoff_rate * spectrum, z, axis=0)
-collar_lost = np.trapezoid(
-    (boundary.spectral_rate - walkoff_rate) * spectrum, z, axis=0
-)
+glass = np.zeros(grid.n_omega)  # the medium's own power loss, 1/m
+glass[active] = 2 * grid.omega[active] * np.imag(medium.index(grid.omega[active])) / C
+drain = boundary.spectral_rate + glass  # every rate that drains a bin
+x = drain * np.diff(z)[:, np.newaxis]
+# x / (1 - exp(-x)), which is 1 where nothing drains
+weight = np.divide(x, -np.expm1(-x), out=np.ones_like(x), where=x > 0)
+fed = weight * (spectrum[1:] - np.exp(-x) * spectrum[:-1])  # G dz between the planes
+drained = np.sum(fed - np.diff(spectrum, axis=0), axis=0)
+per_rate = np.divide(drained, drain, out=np.zeros_like(drained), where=drain > 0)
+lost = walkoff_rate * per_rate
+collar_lost = (boundary.spectral_rate - walkoff_rate) * per_rate
+glass_lost = glass * per_rate
 wavelength = 2 * np.pi * C / np.where(grid.omega > 0, grid.omega, np.inf)
 band = (wavelength >= 150e-9) & (wavelength <= 190e-9)
 share = lost[band].sum() / (
-    lost[band].sum() + collar_lost[band].sum() + spectrum[-1, band].sum()
+    lost[band].sum()
+    + collar_lost[band].sum()
+    + glass_lost[band].sum()
+    + spectrum[-1, band].sum()
 )
 print(f"the walk-off rate took {share:.2f} of the 150-190 nm band")
 ```

docs/guide/numerics/linear-operator.md

@@ -216,6 +216,13 @@ is never clamped to zero. A non-finite index raises rather than passing through.
 itself, before $\Phi$ is integrated — clamping the accumulated exponent instead would make
 the answer depend on the step.
 
+The same clamp keeps a modal run inside the exponent budget that a separable or radial
+run is checked against ([free-space model](../free-space/model.md#exact-longitudinal-operator)): no step within the
+fibre's length forms a factor past $\mathrm{e}^{708}$, below the overflow at
+$\ln(\mathrm{DBL\_MAX}) \approx 709.78$, so a replayed step cannot overflow and
+translation does not check it. Gain, a `max_alpha` raised above the default and the
+absorber's rates, which are added to `D` after the clamp, are outside that bound.
+
 ## The operator of a dispersion-specified fibre
 
 ```{eval-rst}

src/solprop/boundaries.py

@@ -579,8 +579,10 @@ def _walkoff_collar_length(grid: TimeGrid, walkoff: float) -> float:
         \ell_c = \frac{L}{\Delta\tau\,\Gamma}\,\frac{A}{2}\int (1 - W)\,dt.
 
     The integral is taken as :func:`warn_walkoff` takes it, so that the two agree on any
-    grid. For two edges of width ``w`` the Planck taper's symmetry makes it ``w``,
-    giving ``ell_c = 6.3 w L / (Gamma dtau)``.
+    grid: the rectangle sum over the samples, exact on the periodic axis, where a
+    trapezoid would drop the one sample that closes the period. For two edges of width
+    ``w`` the Planck taper's symmetry makes it ``w``, giving
+    ``ell_c = 6.3 w L / (Gamma dtau)``.
 
     - ``walkoff``: the drift over ``grid.length``, s.
     - Returns: the collar length in metres.
@@ -592,7 +594,7 @@ def _walkoff_collar_length(grid: TimeGrid, walkoff: float) -> float:
         )
         raise GridError(msg)
     unit_rate = temporal_collar_rate(grid.t, grid.t_collar, 1.0)
-    integral = float(np.trapezoid(unit_rate, dx=grid.dt))
+    integral = float(np.sum(unit_rate)) * grid.dt
     length = (
         0.5 * integral * grid.length / (MIN_COLLAR_ATTENUATION * walkoff)
     ) * _WALKOFF_ROUNDING_MARGIN
@@ -689,7 +691,7 @@ def absorber(
     >>> bool(np.array_equal(walking.spectral_rate, boundary.spectral_rate))
     True
     >>> round(walking.in_stage_length / boundary.in_stage_length, 3)
-    1.023
+    1.025
     """
     if walkoff is not None and not temporal:
         msg = (
@@ -1372,9 +1374,9 @@ def warn_walkoff(
     if group_delay_span <= 0.0:
         return
     rate = temporal_collar_rate(grid.t, grid.t_collar, length)
-    crossing = (
-        0.5 * float(np.trapezoid(rate, dx=grid.dt)) * grid.length / group_delay_span
-    )
+    # The rectangle sum, not a trapezoid: the axis is periodic, and the sum is exact
+    # for the taper the packet crosses.
+    crossing = 0.5 * float(np.sum(rate)) * grid.dt * grid.length / group_delay_span
     if crossing >= MIN_COLLAR_ATTENUATION:
         return
     warnings.warn(

src/solprop/translate/problem.py

@@ -1430,8 +1430,16 @@ def _check_exponent_budget(self, stepper: StepperSettings) -> None:
         the budget only at its floor ``dz_min``; above it, its controller rejects a
         non-finite trial and shrinks. The radial kind carries a dense operator and
         forms the same inverse factor from it, so it is bounded by
-        :func:`solprop.translate.radial.radial_rate_bounds` instead (SEP 22, Q16);
-        the modal kinds are not checked.
+        :func:`solprop.translate.radial.radial_rate_bounds` instead (SEP 22, Q16).
+
+        The modal kinds are not checked because they need no check: a mode's loss
+        reaches ``D`` through :func:`solprop.linear.dispersion_operator`, whose default
+        :class:`solprop.linear.LossPolicy` clamps the power attenuation to
+        ``2 * MAX_ATTENUATION_EXPONENT / length``, so no step within the run's length
+        forms an amplitude exponent past 708, below ``log(DBL_MAX)``. Outside that
+        reason, and so unchecked here: gain, which the policy lets through; a
+        ``max_alpha`` raised above the default limit; and the absorber's rates, which
+        are added to ``D`` after the clamp.
         """
         operator = self.linear.separable
         geometry = self.transform.geometry

tests/test_boundaries.py

@@ -454,9 +454,12 @@ def test_a_malformed_absorber_is_rejected(field_grid):
 
 
 def crossing_nepers(grid, group_delay_span, boundary_n=DEFAULT_BOUNDARY_N):
-    """The field attenuation a walking packet accumulates crossing the collar."""
+    """The field attenuation a walking packet accumulates crossing the collar.
+
+    The rate summed over the periodic axis, which the packet crosses whole.
+    """
     rate = absorber(grid, boundary_N=boundary_n).temporal_rate
-    return 0.5 * float(np.trapezoid(rate, dx=grid.dt)) * grid.length / group_delay_span
+    return 0.5 * float(np.sum(rate)) * grid.dt * grid.length / group_delay_span
 
 
 def test_the_walk_off_warning_fires_when_the_absorber_is_too_weak(field_grid):
@@ -468,6 +471,32 @@ def test_the_walk_off_warning_fires_when_the_absorber_is_too_weak(field_grid):
         warn_walkoff(field_grid, span)
 
 
+@pytest.mark.parametrize("grid_name", GRIDS)
+def test_the_walk_off_warning_reads_the_periodic_sum_of_the_collar(grid_name, request):
+    """The threshold is the closed form ``6.3 w L / (Gamma ell)``, not one sample short.
+
+    Over both edges the rate integrates to ``12.6 w / ell`` exactly on the periodic
+    axis; a trapezoid over the samples would read ``w - dt`` and fire at a span
+    ``dt / w`` short of the threshold.
+    """
+    grid = request.getfixturevalue(grid_name)
+    width = edge_width(grid)
+    reference = grid.length / DEFAULT_BOUNDARY_N
+    threshold = (
+        0.5
+        * TEMPORAL_ALPHA_LENGTH
+        * width
+        * grid.length
+        / (reference * MIN_COLLAR_ATTENUATION)
+    )
+    assert grid.dt / width > 1e-4
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        warn_walkoff(grid, threshold * (1.0 - 1e-6))
+    with pytest.warns(SolpropWarning, match="wrapping around"):
+        warn_walkoff(grid, threshold * (1.0 + 1e-6))
+
+
 def test_the_walk_off_warning_is_silent_when_the_absorber_is_strong_enough(field_grid):
     span = 100e-15
     assert crossing_nepers(field_grid, span) > MIN_COLLAR_ATTENUATION
@@ -554,20 +583,15 @@ def edge_width(grid):
 def test_a_walkoff_sets_the_length_the_derivation_gives(field_grid):
     """``ell_c = 6.3 w L / (Gamma dtau)``, on the integral ``warn_walkoff`` reads.
 
-    The integral of ``1 - W`` over both edges is ``w`` by the Planck taper's symmetry;
-    the trapezoid on the axis, which stops one sample short of closing the period,
-    reads ``w - dt``. So the bound is that one sample, and the sampled form is exact.
+    The integral of ``1 - W`` over both edges is ``w`` by the Planck taper's symmetry,
+    and the rectangle sum over the periodic axis reads it exactly, so the bound is the
+    rounding margin of the sizing.
     """
     drift = 3e-12
     width = edge_width(field_grid)
     collar = absorber(field_grid, walkoff=drift)
     closed_form = 6.3 * width * field_grid.length / (MIN_COLLAR_ATTENUATION * drift)
-    assert collar.in_stage_length == pytest.approx(
-        closed_form, rel=2.0 * field_grid.dt / width
-    )
-    assert collar.in_stage_length == pytest.approx(
-        closed_form * (width - field_grid.dt) / width, rel=1e-6
-    )
+    assert collar.in_stage_length == pytest.approx(closed_form, rel=1e-8)
     assert float(collar.temporal_rate.max()) == pytest.approx(
         TEMPORAL_ALPHA_LENGTH / collar.in_stage_length, rel=1e-3
     )
@@ -815,11 +839,12 @@ def test_a_walkoff_collar_removes_light_of_the_stated_walk_off(backend, stretch)
     """A packet of the stated drift leaves a millionth of its power, and no more.
 
     The packet crosses both edges rigidly, so its power falls by exactly the integral
-    of the rate over the ring divided by the speed: ``exp(-2 Gamma)`` of the sizing,
-    corrected by the one sample the sizing's trapezoid leaves out of the periodic sum.
+    of the rate over the ring divided by the speed: ``exp(-2 Gamma)``, which the sizing
+    reads from the same periodic sum.
     A drift of half the collar's or twice its stated value shows the exponent is
     linear in the speed, and that the light has no other way out. The bound on the
-    exponent is the controller's ``rtol``; measured 6e-9 on both backends.
+    exponent is the controller's ``rtol``; measured 2.4e-9 at the stated drift on
+    ``numpy_ref``.
     """
     if backend == "native":
         pytest.importorskip("solprop_native")
@@ -830,7 +855,7 @@ def test_a_walkoff_collar_removes_light_of_the_stated_walk_off(backend, stretch)
     collar = absorber(grid, walkoff=drift / stretch)
     closed_period = np.sum(temporal_collar_rate(grid.t, grid.t_collar, 1.0)) * grid.dt
     expected = closed_period * grid.length / (collar.in_stage_length * drift)
-    assert expected == pytest.approx(2.0 * MIN_COLLAR_ATTENUATION / stretch, rel=0.02)
+    assert expected == pytest.approx(2.0 * MIN_COLLAR_ATTENUATION / stretch, rel=1e-8)
     power = transmitted_power(grid, collar, drift, backend=backend)
     assert -math.log(power) == pytest.approx(expected, rel=1e-7)
 
@@ -846,8 +871,7 @@ def test_without_a_temporal_collar_the_same_packet_loses_nothing():
 def test_the_split_scheme_removes_the_same_light():
     """The split's factors telescope, within its crossing limit, to the same exponent.
 
-    The split's first-order error in the step and the sizing's left-out sample, 0.7 %
-    here, set the bound; measured 0.5 %.
+    The split's first-order error in the step sets the bound; measured 0.23 %.
     """
     grid = EnvelopeGrid.from_limits(1.0, 800e-9, (600e-9, 1200e-9), 2e-12, thg=False)
     drift = 0.8 * grid.n_t * grid.dt

tests/test_spectral_walkoff_recipe.py

@@ -294,6 +294,56 @@ def test_the_per_band_loss_read_from_planes_is_the_loss(backend):
     )
 
 
+def interval_drained(spectrum, z, drain):
+    """What every drain took from each bin, read exactly between the planes.
+
+    The guide's formula: with ``dS/dz = G - a S`` and ``G`` constant over an interval,
+    ``G dz = x (S1 - S0 exp(-x)) / (1 - exp(-x))``, ``x = a dz``, and the drains took
+    ``G dz - (S1 - S0)``. Returns that sum over the intervals, divided by ``a``, so a
+    rate times it is the rate's share.
+    """
+    x = drain * np.diff(z)[:, np.newaxis]
+    weight = np.divide(x, -np.expm1(-x), out=np.ones_like(x), where=x > 0)
+    fed = weight * (spectrum[1:] - np.exp(-x) * spectrum[:-1])
+    drained = np.sum(fed - np.diff(spectrum, axis=0), axis=0)
+    return np.divide(drained, drain, out=np.zeros_like(drained), where=drain > 0)
+
+
+def test_the_guide_s_interval_formula_is_exact_on_few_planes():
+    """The bins the rate is for fall far between planes; the formula still reads them.
+
+    On four planes ``a dz`` reaches 2.7 nepers at the launch's bins, where the
+    trapezoid overcounts. The interval formula is exact for a bin that only decays, so
+    the controller's ``rtol`` and rounding are the bound, and four planes and
+    thirty-two agree to it.
+    """
+    rate = recipe_rate(math.log(1 / SURVIVAL) / LENGTH)
+    glass = np.zeros(GRID.n_omega)
+    active = GRID.active_mask
+    glass[active] = (
+        2 * GRID.omega[active] * np.imag(MEDIUM.index(GRID.omega[active])) / C
+    )
+    whole = rate + glass
+    field = launch([(carrier_at(ONSET + 0.5 * PROTECTED), 1e-9, 0.0)])
+    readings = []
+    for planes in (4, 32):
+        output = run(spectral_only(rate), field, planes=planes)
+        spectrum = np.asarray(output["stats/beam_spectrum"])
+        z = np.asarray(output["z"])
+        readings.append(rate * interval_drained(spectrum, z, whole))
+    few, many = readings
+    weak = (spectrum[0] > 1e-3 * spectrum[0].max()) & (rate > 0.0)
+    assert np.count_nonzero(weak) > 3
+    assert float(np.max(whole[weak])) * LENGTH / 4 > 2.0
+    exact = (
+        spectrum[0, weak]
+        * (rate[weak] / whole[weak])
+        * -np.expm1(-whole[weak] * LENGTH)
+    )
+    np.testing.assert_allclose(few[weak], exact, rtol=1e-6)
+    np.testing.assert_allclose(few[weak], many[weak], rtol=1e-6)
+
+
 def test_a_band_can_lose_most_of_itself_while_the_run_s_energy_barely_moves():
     """Why the check is per band: the total reads ``1e-5`` while the band reads tens of
     per cent."""