nereids_physics/

ikeda_carpenter.rs

//! Ikeda–Carpenter analytical moderator resolution model.
//!
//! A third instrument-resolution model alongside the analytical Gaussian
//! ([`crate::resolution::ResolutionParams`]) and the Monte-Carlo tabulated
//! kernel ([`crate::resolution::TabulatedResolution`]). It exists to settle a
//! methodological dispute about the VENUS instrument resolution: one camp
//! trusts the MC-simulated tabulated kernel (UDR/FTS file); the instrument
//! scientist distrusts the unproven MC and prefers an analytical
//! Ikeda–Carpenter moderator model. NEREIDS implements IC as a first-class
//! model so all three can be cross-validated against synthetic loop-closure
//! and real VENUS data.
//!
//! # Physics — the Ikeda–Carpenter pulse
//!
//! Reference: S. Ikeda & J. M. Carpenter, *Nucl. Instrum. Methods* **A239**
//! (1985) 536–544. The neutron emission-time distribution from a pulsed
//! spallation moderator is
//!
//! ```text
//!   I(τ) = (1−R)·g₃(τ;α)  +  R·[g₃(·;α) ⊛ β e^{−β·}](τ),    τ ≥ 0
//! ```
//!
//! where the prompt (slowing-down) term is a Gamma/Erlang density of shape 3:
//!
//! ```text
//!   g₃(τ;α) = α³ τ² e^{−ατ} / 2          (mode 2/α, mean 3/α, ∫ = 1)
//! ```
//!
//! and the delayed (storage) term is `g₃` convolved with an exponential of
//! rate β. The convolution has the closed form
//!
//! ```text
//!   g₃(·;α) ⊛ β e^{−β·} (τ) = β (α/γ)³ [ e^{−βτ} − e^{−ατ}(1 + γτ + ½γ²τ²) ],  γ = α−β
//! ```
//!
//! SAMMY-RPI χ²+double-exp moderator form). Both terms are individually
//! unit-area, so `I` is unit-area for any α,β>0 and 0≤R≤1, with first moment
//! `⟨τ⟩ = 3/α + R/β`. The asymmetry — sharp rise, long tail toward larger τ
//! (later TOF, lower apparent energy) — is the physical origin of the
//! asymmetric MC kernel.
//!
//! ## Parameters and their energy dependence
//!
//! - `α(E)` [1/µs]: fast moderation/leakage rate; sets the prompt width.
//!   Leading epithermal scaling `α ∝ √E` (Mantid `α = 1/(α₀+α₁λ)`, λ ∝ 1/√E).
//! - `β(E)` [1/µs]: slow storage rate; sets the delayed tail. Constant by
//!   default; optionally energy-dependent through an [`EnergyLaw`].
//! - `R(E)`, 0 ≤ R ≤ 1: storage mixing fraction; `R ≈ exp(−E_meV/κ)` → **R→0 in
//!   the 1–200 eV resonance regime**, so IC there is dominated by the
//!   one-parameter prompt Gamma(3, α(E)) term.
//!
//! Parameters are **fixed** when the instrument scientist provides them, or
//! **fit** from a known calibration foil otherwise (general case).
//!
//! ## Time → energy kernel and centering
//!
//! For a flight path `L`, a resonance at `E_r` has nominal TOF
//! `t_r = TOF_FACTOR·L/√E_r`. A neutron with emission delay τ arrives at TOF
//! `t_r+τ`, apparent energy `E' = (TOF_FACTOR·L/(t_r+τ))²` — that is the kernel's
//! *definition* (its positive-τ tail is delayed emission). Sampling `I(τ)` on
//! a τ-grid and mapping to TOF-offsets yields exactly the `(offset, weight)`
//! kernel representation that [`crate::resolution::TabulatedResolution`]
//! consumes — so IC rides the *same* verified broadening machinery, whose
//! *application* is the convolution gather: the broadened value at measured TOF
//! `t` reads theory at `t−τ` (a neutron measured at `t` with delay τ really flew
//! `t−τ`; see `resolution::broaden` and SAMMY `udr/mudr4.f90` `Ud_Convolute`).
//! At apply time
//! `interpolated_kernel` blends the two bracketing reference kernels as a
//! width-normalized shape blend (offsets scaled to the geometrically
//! interpolated width, shapes merged on the union grid) — unequal point counts
//! from IC's per-kernel tail trimming no longer trigger a nearest-reference
//! fallback, so the between-reference width follows the physical power law
//! smoothly. IC also synthesizes a dense reference
//! grid (default 64 energies), so the between-reference error is negligible. The
//! kernel keeps its offsets on the **emission clock** (0 = pulse start), which
//! is where the Ikeda–Carpenter function itself places them; a loaded UDR file
//! is anchored on its peak instead and `interpolated_kernel` re-centers
//! neither. Because the IC pulse is right-skewed, its *mean* lags its mode by
//! ~1/α(E) in TOF, so the centroid — and even the minimum — of a broadened
//! resonance shifts toward **lower apparent energy** by an α(E)-dependent amount
//! (order 1e-2 eV, ~1e-3 relative, for α≈1.5 in the eV regime; larger toward
//! lower energy). So it *does* move the broadened dip off the nominal energy,
//! by a small amount.
//!
//! For the prompt-only law `α(E) = a0·√E + a1` with `R≈0`, the lag `~1/α(E)` is
//! **exactly** the `1/√E` basis of a flight-path (`L_scale`) error *iff* `a1 = 0`:
//! then the centroid offset scales as `c/√E`, and an `L → L(1+δ)` change shifts
//! every TOF by `δ·K·L/√E`, also `∝ 1/√E`, so the two are degenerate. With `a1 ≠ 0`
//! the lag `1/(a0√E + a1)` only *approximately* follows that basis (and the
//! storage term `R/β`, negligible in the eV regime, adds a further small
//! departure). To leading order, then, the lag is **confounded with the energy
//! scale**, and is handled by the SHARED `(t0, L_scale)` energy scale, not by any
//! per-family knob:
//! - **Run-time fitting** fits the `t0`/`L_scale` energy-scale, which absorbs the
//!   constant-`L` part of the lag exactly (same basis); only the *shape*
//!   (skew/tail) of the asymmetry is not absorbable by position.
//! - **Resolution calibration** (the `nereids-fitting` calibrator) **pins** the
//!   energy scale by default — a pure shape/width fit. It can optionally fit a
//!   SHARED `(t0, L_scale)` under a metrology prior (`with_position_prior`) for
//!   joint energy-scale / identifiability work. A *free, per-family* position knob
//!   is deliberately NOT used: because the lag is the same basis as `L_scale`, a
//!   free position lets a wrong (symmetric) family imitate the asymmetric shift
//!   and erodes the model-selection χ² (the discriminator is then only the
//!   position-independent skew/tail). See `nereids-fitting`'s
//!   `free_l_scale_absorbs_asymmetric_lag_and_erodes_discrimination`.
//!
//! `ic_centering_shifts_broadened_symmetric_dip_with_alpha` quantifies the bare
//! mode→centroid shift; the calibrator's `fit_t0_recovers_injected_energy_scale_shift`
//! checks the energy-scale fit recovers an injected offset. Re-centering the kernel
//! on its centroid (a future convention change spanning the UDR path too) would
//! remove the shift at the source.
//!
//! ## Optional instrument convolutions
//!
//! The full instrument function folds the moderator with a proton-burst
//! (Gaussian σ) and a chopper/channel (triangle, FWHM) term. Both are optional
//! here (`None` ⇒ omitted). Detector-time bin probabilities take the triangle
//! only and refuse a burst. They also take a detector tail, which records a
//! fraction of the counts late by an exponential delay; the kernel table
//! refuses one.
//!
//! **Provenance of the triangle (`channel_fwhm_us`).** SAMMY broadens for the
//! accelerator burst either as a Gaussian of FWHM `DELTAG` (SAMMY Manual R8
//! Sec. III.C.1.a, eq. III C1 a.12) or as a square pulse of width `BURST`
//! (Sec. III.C.2.a). At SNS the proton pulse delivered to the target is shaped
//! by the accumulator ring (Proton Storage Ring, PSR) into an approximately
//! triangular ~700 ns base — FWHM ≈ 350 ns — which is what the VENUS tabulated
//! FTS kernel header records as "folded triang FWHM 350 ns PSR". NEREIDS folds
//! that PSR triangle via `channel_fwhm_us` (symmetric triangle, half-base =
//! FWHM, `triangle_kernel`). Note: a tabulated file whose header says the
//! triangle is already "folded" in must NOT be double-counted against an IC
//! model that also applies it (the `nereids-fitting` calibrator therefore
//! applies its `psr_fwhm_ns` fold to the IC family only, never to
//! tabulated/UDR kernels).

use std::ops::RangeInclusive;
use std::sync::Arc;

use crate::resolution::{ResolutionParseError, TOF_FACTOR, TabulatedResolution};

/// de Broglie wavelength factor: λ (Å) = `LAMBDA_ANGSTROM_FACTOR` / √(E in eV).
///
/// `λ = h/√(2·m_n·E)`; with CODATA 2018 `h`, `m_n` and the 2019-SI eV this is
/// `0.285993 Å·√eV`. Used only by [`EnergyLaw::InverseLambda`]; its precise
/// value folds into the fitted `α₀,α₁`, so high precision is not load-bearing.
const LAMBDA_ANGSTROM_FACTOR: f64 = 0.285_993;

/// Rates below this (1/µs) are clamped to keep the pulse well-defined.
const MIN_RATE: f64 = 1e-9;

/// The detector tail's mean delays `τ_D`, in µs, that the pulse accepts.
pub const DETECTOR_TAIL_DECAY_US: RangeInclusive<f64> = 1e-9..=1e9;

/// Relative-to-peak weight below which kernel tails are trimmed.
const TRIM_REL: f64 = 1e-7;

/// Default number of log-spaced reference energies in a synthesized table.
pub const DEFAULT_N_ENERGIES: usize = 64;

/// Default number of τ-samples spanning the prompt core of each kernel.
pub const DEFAULT_N_TAU: usize = 600;

/// Minimum accepted `n_tau` (τ-samples across the prompt core). Doubles as
/// the module's prompt-core **resolution floor**: when [`MAX_TAU_SAMPLES`]
/// widens the τ-step (see [`tau_geometry`]), the step may never exceed
/// `fast_reach / (MIN_N_TAU − 1)` — the coarsest prompt sampling
/// [`IkedaCarpenter::new`] has ever accepted as valid via its `n_tau ≥ 8`
/// check.
const MIN_N_TAU: usize = 8;

/// Minimum samples per side SPANNING a sampled channel triangle
/// (`dtau ≤ FWHM / 3`, half-base = FWHM). At the exactly-admitted boundary
/// `dtau = FWHM/3` the per-side samples sit at `{FWHM/3, 2FWHM/3, FWHM}` —
/// triangle weights `{2/3, 1/3, 0}` — so each side carries ≥ 2 strictly
/// interior (nonzero) samples while the endpoint lands ON the triangle
/// zero; the sampled fold is distinctly non-delta (discrete variance
/// 4·FWHM²/27 ≈ 89 % of the analytic FWHM²/6). Coarser steps degenerate
/// the discrete triangle toward the exact delta `[0, 1, 0]`, silently
/// erasing a requested fold — [`tau_geometry`] rejects that instead
/// (strictly: `capped_step > FWHM/3`).
const TRI_MIN_SAMPLES_PER_SIDE: f64 = 3.0;

const GAUSS_REACH_SIGMAS: f64 = 8.0;

/// Prompt-tail reach in e-folds: the τ-grid spans `FAST_REACH_E_FOLDS / α`
/// so the prompt Gamma(3) tail `τ²e^{−ατ}` is < ~1e-8 of peak at the edge.
const FAST_REACH_E_FOLDS: f64 = 18.0;

/// Slow/storage-tail reach in e-folds (`SLOW_REACH_E_FOLDS / β` when storage
/// is active). Sits AT the trim horizon: `e^{−16} ≈ 1.1e-7 ≈` [`TRIM_REL`]
/// (`ln(1/TRIM_REL) ≈ 16.1`), so the reach cannot be truncated shorter
/// without discarding tail weight the [`TRIM_REL`] trim would keep — no
/// hidden approximation lives in this constant.
const SLOW_REACH_E_FOLDS: f64 = 16.0;

/// Storage fractions below this are treated as "no storage tail" when sizing
/// the τ-grid. DELIBERATELY two decades below [`TRIM_REL`], not derived from
/// it: whether a slow tail actually survives the [`TRIM_REL`] trim depends on
/// the α/β contrast (the trim threshold is relative to the prompt-dominated
/// peak, and the slow term's peak weight scales with both R and β/α), so no
/// single constant derived from [`TRIM_REL`] is exact for every admitted
/// (α, β) — treating a larger R as absent would be a hidden approximation.
/// The margin's only consequence is conservatism: for R ∈ (1e-9, ~1e-7) the
/// τ-grid is sized — and the [`MAX_TAU_SAMPLES`] cap gate applied — for a
/// slow tail whose weight the trim then discards anyway, spending samples
/// (or rejecting a configuration) for physics that cannot appear in the
/// kernel. It can never drop tail weight the trim would have kept.
const R_NEGLIGIBLE: f64 = 1e-9;

/// Cap on the τ-sample count spanning the PULSE BODY (`[0, τ_max]`) of one
/// synthesized kernel. The τ-step is anchored to the prompt core and refined
/// to resolve active folds (see [`tau_geometry`]), so a long storage tail
/// (β ≪ α with R > 0) grows the sample COUNT rather than the step; this cap
/// bounds that growth (CPU/memory) by widening the step to
/// `tau_max / (MAX_TAU_SAMPLES − 1)` — but never past the resolution floor
/// (prompt core: `fast_reach / (MIN_N_TAU − 1)`; triangle: FWHM /
/// [`TRI_MIN_SAMPLES_PER_SIDE`]; burst: σ). A parameter/grid combination
/// whose floor cannot be met within the cap is REJECTED loudly by
/// [`IkedaCarpenter::new`] instead of silently under-sampled. The actual
/// guarantee is on the pulse body only: the symmetric burst/channel fold
/// margin (± `GAUSS_REACH_SIGMAS·σ + FWHM` at the resolved step) adds its
/// samples ON TOP of the cap, so the final grid can exceed
/// `MAX_TAU_SAMPLES` by the margin sample count.
const MAX_TAU_SAMPLES: usize = 8192;

#[inline]
fn h_over_cube_taylor(u: f64) -> f64 {
    exponential_tail_series(u, 3)
}

#[inline]
fn exponential_tail_series(u: f64, shift: u32) -> f64 {
    let mut sum = 0.0;
    let mut power = 0.5;
    for n in 0..64_u32 {
        let next = sum + power / f64::from(n + shift);
        if next == sum {
            break;
        }
        sum = next;
        power *= -u / f64::from(n + 1);
    }
    sum
}

/// Ikeda–Carpenter moderator emission density `I(τ)`.
///
/// `τ` in µs, rates `α,β` in 1/µs, mixing `r ∈ [0,1]`. Returns 0 for `τ < 0`.
/// Unit-area over `τ ∈ [0,∞)`; first moment `3/α + r/β`. **NaN-free** for all
/// `α,β > 0` (including `α ≈ β` and `β ≫ α`): the slow/storage term is evaluated
/// from the bounded bracket `e^{−βτ} − e^{−ατ}(1+u+½u²)` (both exponentials ≤ 1
/// for `τ,α,β > 0`, so no `e^{|u|}` overflow), falling back to a Taylor limit
/// near `u = α−β·τ → 0` where that bracket cancels.
#[must_use]
pub fn ic_pulse(alpha: f64, beta: f64, r: f64, tau: f64) -> f64 {
    // Same domain contract as ic_cdf: non-finite parameters return NaN so
    // garbage in stays visible instead of flooring into a plausible pulse.
    if !alpha.is_finite() || !beta.is_finite() || !r.is_finite() {
        return f64::NAN;
    }
    if !tau.is_finite() || tau < 0.0 {
        return 0.0;
    }
    let alpha = alpha.max(MIN_RATE);
    let beta = beta.max(MIN_RATE);
    let at = alpha * tau;
    // prompt Gamma(3,α): α³τ²/2 · e^{−ατ} = (α/2)(ατ)² e^{−ατ}
    let fast = 0.5 * alpha * at * at * (-at).exp();
    if r <= 0.0 {
        return fast;
    }
    (1.0 - r) * fast + r * storage_density(alpha, beta, tau)
}

fn storage_density(alpha: f64, beta: f64, tau: f64) -> f64 {
    // slow/storage term = β(α/γ)³[e^{−βτ} − e^{−ατ}(1+u+½u²)], u = γτ = (α−β)τ.
    let u = (alpha - beta) * tau;
    let coeff = beta * alpha.powi(3) * tau.powi(3);
    if u.abs() < 0.05 {
        // α ≈ β: bracket/u³ → e^{−βτ}·h(u)/u³ (Taylor); avoids 0/0 cancellation.
        coeff * (-beta * tau).exp() * h_over_cube_taylor(u)
    } else {
        // Bounded form: both exponentials are ≤ 1, so β ≫ α (u ≪ 0) cannot
        // overflow (the old `e^{−βτ}·h(u)` factored an `e^{|u|}` → 0·∞ = NaN).
        let bracket = (-beta * tau).exp() - (-alpha * tau).exp() * (1.0 + u + 0.5 * u * u);
        coeff * bracket / (u * u * u)
    }
}

fn ic_pulse_slope(alpha: f64, beta: f64, r: f64, tau: f64) -> f64 {
    if !tau.is_finite() || tau < 0.0 {
        return 0.0;
    }
    let alpha = alpha.max(MIN_RATE);
    let beta = beta.max(MIN_RATE);
    let at = alpha * tau;
    let prompt = 0.5 * alpha * at * at * (-at).exp();
    let prompt_slope = 0.5 * alpha * alpha * at * (2.0 - at) * (-at).exp();
    if r <= 0.0 {
        return prompt_slope;
    }
    (1.0 - r) * prompt_slope + r * beta * (prompt - storage_density(alpha, beta, tau))
}

/// Cumulative probability of the Ikeda–Carpenter moderator pulse.
///
/// Returns `P(U <= tau)` for moderator delay `U`. `tau` is in µs and rates
/// are in 1/µs. The prompt term is the Gamma(3, α) cumulative distribution;
/// the storage term is the cumulative distribution of Gamma(3, α) plus an
/// independent exponential delay with rate β.
///
/// The expression is evaluated without subtracting nearly equal exponentials
/// when `α ≈ β`. This is the bin-integral companion to [`ic_pulse`].
///
/// Domain contract: non-finite `alpha`, `beta`, or `r` returns NaN so garbage
/// in stays visible; finite non-positive rates are floored to `MIN_RATE`
/// (1e-9 µs⁻¹) and `r <= 0` disables the storage term, matching
/// [`ic_pulse`]. A NaN `tau`
/// returns 0 — `tau` is the integration coordinate, not a parameter, and a
/// non-arriving coordinate contributes no mass.
#[must_use]
pub fn ic_cdf(alpha: f64, beta: f64, r: f64, tau: f64) -> f64 {
    if !alpha.is_finite() || !beta.is_finite() || !r.is_finite() {
        return f64::NAN;
    }
    if tau.is_nan() || tau <= 0.0 {
        return 0.0;
    }
    if tau == f64::INFINITY {
        return 1.0;
    }

    let alpha = alpha.max(MIN_RATE);
    let beta = beta.max(MIN_RATE);
    let at = alpha * tau;
    let fast = gamma3_cdf(at);
    if r <= 0.0 {
        return fast;
    }

    // Far in the tails the general expressions overflow — `at³` or `u³`
    // reach f64::INFINITY and produce inf·0 / inf/inf = NaN — while the
    // limits are exact to full f64 precision, so return them directly.
    // For at ≥ 1e100 the prompt CDF is 1 and the storage correction is
    // exp(−βτ) with relative error ≤ 3βτ/at; for βτ ≥ 1e100 the storage
    // delay is instantaneous and the correction is 0.
    const CDF_TAIL_LIMIT: f64 = 1.0e100;
    if at >= CDF_TAIL_LIMIT {
        return (1.0 - r * (-beta * tau).exp()).clamp(0.0, 1.0);
    }
    if beta * tau >= CDF_TAIL_LIMIT {
        return fast;
    }

    // The storage CDF is the prompt CDF minus the positive survival
    // correction below. Written this way, the α=β limit is Gamma(4, α).
    let u = (alpha - beta) * tau;
    let correction = if u.abs() < 0.05 {
        at.powi(3) * (-beta * tau).exp() * h_over_cube_taylor(u)
    } else {
        // Bounded form: neither exponential can overflow even when β >> α.
        let bracket = (-beta * tau).exp() - (-alpha * tau).exp() * (1.0 + u + 0.5 * u * u);
        at.powi(3) * bracket / u.powi(3)
    };
    (fast - r * correction).clamp(0.0, 1.0)
}

fn ic_cdf_slopes(alpha: f64, beta: f64, r: f64, tau: f64) -> [f64; 3] {
    if !alpha.is_finite() || !beta.is_finite() || !r.is_finite() {
        return [f64::NAN; 3];
    }
    if tau.is_nan() || tau <= 0.0 || tau == f64::INFINITY {
        return [0.0; 3];
    }
    let at = alpha * tau;
    let fast = 0.5 * tau * at * at * (-at).exp();
    const CDF_TAIL_LIMIT: f64 = 1.0e100;
    if at >= CDF_TAIL_LIMIT {
        let storage = (-beta * tau).exp();
        return [0.0, r * tau * storage, -storage];
    }
    if beta * tau >= CDF_TAIL_LIMIT {
        return [fast, 0.0, 0.0];
    }
    let u = (alpha - beta) * tau;
    let (correction, slope) = if u.abs() < 0.05 {
        let storage = (-beta * tau).exp();
        (
            at.powi(3) * storage * h_over_cube_taylor(u),
            -tau * storage * exponential_tail_series(u, 4),
        )
    } else {
        let prompt = (-alpha * tau).exp();
        let bracket = (-beta * tau).exp() - prompt * (1.0 + u + 0.5 * u * u);
        (
            at.powi(3) * bracket / u.powi(3),
            tau * (prompt / (2.0 * u) - 3.0 * bracket / u.powi(4)),
        )
    };
    let moved = at.powi(3) * slope;
    [
        fast - r * (3.0 * correction / alpha + moved),
        r * (tau * correction + moved),
        -correction,
    ]
}

const GAUSS_LEGENDRE_16: [(f64, f64); 8] = [
    (0.095_012_509_837_637_44, 0.189_450_610_455_068_59),
    (0.281_603_550_779_258_9, 0.182_603_415_044_923_6),
    (0.458_016_777_657_227_37, 0.169_156_519_395_002_62),
    (0.617_876_244_402_643_8, 0.149_595_988_816_576_65),
    (0.755_404_408_355_003_1, 0.124_628_971_255_533_95),
    (0.865_631_202_387_831_6, 0.095_158_511_682_493),
    (0.944_575_023_073_232_6, 0.062_253_523_938_647_61),
    (0.989_400_934_991_649_9, 0.027_152_459_411_754_055),
];

fn fold_nodes(alpha: f64, fwhm_us: f64, delay_us: f64, mut node: impl FnMut(f64, f64)) {
    let end = delay_us.min(fwhm_us);
    let finest = 2.0 / alpha.max(MIN_RATE);
    for (lo, hi) in [(-fwhm_us, end.min(0.0)), (0.0, end)] {
        let (mut right, mut width) = (hi, finest);
        while right > lo {
            let left = (right - width).max(lo);
            let (centre, half) = (0.5 * (left + right), 0.5 * (right - left));
            for &(position, weight) in &GAUSS_LEGENDRE_16 {
                for s in [centre - half * position, centre + half * position] {
                    node(s, half * weight);
                }
            }
            right = left;
            width *= 2.0;
        }
    }
}

fn triangle_fold<const N: usize>(
    alpha: f64,
    fwhm_us: f64,
    delay_us: f64,
    g: impl Fn(f64) -> [f64; N],
) -> [f64; N] {
    let mut total = [0.0; N];
    fold_nodes(alpha, fwhm_us, delay_us, |s, weight| {
        for (sum, value) in total.iter_mut().zip(g(delay_us - s)) {
            *sum += weight * (fwhm_us - s.abs()) / (fwhm_us * fwhm_us) * value;
        }
    });
    total
}

fn folded_slopes(pulse: &Recorded, fwhm_us: f64, delay_us: f64) -> [f64; 4] {
    let mut slopes = [0.0; 4];
    fold_nodes(pulse.alpha, fwhm_us, delay_us, |s, weight| {
        let x = delay_us - s;
        let t = (s / fwhm_us).abs();
        let triangle = weight * (1.0 - t) / fwhm_us;
        let values = pulse.slopes(x);
        for (slope, value) in slopes.iter_mut().zip(&values[..3]) {
            *slope += triangle * value;
        }
        let moment = (1.0 / 6.0 - 0.5 * t * t + t.powi(3) / 3.0) / (2.0 * fwhm_us);
        slopes[3] += weight * moment * values[3];
    });
    slopes
}

/// The pulse at one energy as the detector records it: rates `α` and `β`,
/// storage fraction `R`, and the detector tail's fraction `f` and rate
/// `1/τ_D`.
#[derive(Debug, Clone, Copy)]
struct Recorded {
    alpha: f64,
    beta: f64,
    r: f64,
    tail: Option<(f64, f64)>,
}

impl Recorded {
    fn delayed(&self) -> Option<(f64, f64)> {
        self.tail.filter(|&(fraction, _)| fraction != 0.0)
    }

    fn cdf(&self, tau: f64) -> f64 {
        let plain = ic_cdf(self.alpha, self.beta, self.r, tau);
        match self.delayed() {
            None => plain,
            Some((f, rate)) => ((1.0 - f) * plain + f * self.late(rate, tau)[0]).clamp(0.0, 1.0),
        }
    }

    fn density(&self, tau: f64) -> f64 {
        let plain = ic_pulse(self.alpha, self.beta, self.r, tau);
        match self.delayed() {
            None => plain,
            Some((f, rate)) => (1.0 - f) * plain + f * self.late(rate, tau)[5],
        }
    }

    /// `[∂/∂α, ∂/∂β, ∂/∂R]` of [`Self::cdf`], and the slope of
    /// [`Self::density`].
    fn slopes(&self, tau: f64) -> [f64; 4] {
        let [a, b, c] = ic_cdf_slopes(self.alpha, self.beta, self.r, tau);
        let plain = [a, b, c, ic_pulse_slope(self.alpha, self.beta, self.r, tau)];
        match self.delayed() {
            None => plain,
            Some((f, rate)) => {
                let late = self.late(rate, tau);
                let late = [late[1], late[2], late[3], late[6]];
                std::array::from_fn(|n| (1.0 - f) * plain[n] + f * late[n])
            }
        }
    }

    /// `[∂/∂f, ∂/∂τ_D]` of [`Self::cdf`], 0 without a tail.
    fn tail_slopes(&self, tau: f64) -> [f64; 2] {
        let Some((f, rate)) = self.tail else {
            return [0.0; 2];
        };
        let late = self.late(rate, tau);
        [
            late[0] - ic_cdf(self.alpha, self.beta, self.r, tau),
            -f * rate * rate * late[4],
        ]
    }

    fn late(&self, rate: f64, tau: f64) -> [f64; 7] {
        delayed(self.alpha, self.beta, self.r, rate, tau)
    }
}

/// The relative distance from `β` within which [`delayed`] interpolates in
/// the tail's rate.
const NEAR_RATE: f64 = 1e-3;

/// The pulse delayed by an exponential of rate `rate`, at delay `tau`:
/// `[G, ∂G/∂α, ∂G/∂β, ∂G/∂R, ∂G/∂rate, G′, G″]` with
/// `G = S_λ + Rλ(S_β − S_λ)/(λ − β)`, `λ` the rate and `S_x` the cumulative
/// probability of Gamma(3, α) plus an exponential delay of rate `x`.  Within
/// [`NEAR_RATE`] of `β`, each is the cubic through its values at the rates
/// `β(1 ± NEAR_RATE)` and `β(1 ± 2·NEAR_RATE)`.
fn delayed(alpha: f64, beta: f64, r: f64, rate: f64, tau: f64) -> [f64; 7] {
    let gap = (rate - beta) / (NEAR_RATE * beta);
    if gap.abs() >= 1.0 {
        return delayed_apart(alpha, beta, r, rate, tau);
    }
    let mut total = [0.0; 7];
    for (node, weight) in [
        (-2.0, (gap * gap - 1.0) * (gap - 2.0) / -12.0),
        (-1.0, (gap + 2.0) * (gap - 1.0) * (gap - 2.0) / 6.0),
        (1.0, (gap + 2.0) * (gap + 1.0) * (gap - 2.0) / -6.0),
        (2.0, (gap + 2.0) * (gap * gap - 1.0) / 12.0),
    ] {
        let values = delayed_apart(alpha, beta, r, beta * (1.0 + node * NEAR_RATE), tau);
        for (sum, value) in total.iter_mut().zip(values) {
            *sum += weight * value;
        }
    }
    total
}

fn delayed_apart(alpha: f64, beta: f64, r: f64, rate: f64, tau: f64) -> [f64; 7] {
    if tau.is_nan() || tau <= 0.0 {
        return [0.0; 7];
    }
    if tau == f64::INFINITY {
        return [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
    }
    let alpha = alpha.max(MIN_RATE);
    let at = alpha * tau;
    let (whole, prompt) = (gamma3_cdf(at), 0.5 * alpha * at * at * (-at).exp());
    let storage = |x: f64| {
        let [by_alpha, by_rate, lost] = ic_cdf_slopes(alpha, x, 1.0, tau);
        let density = -x * lost;
        [
            whole + lost,
            by_alpha,
            by_rate,
            density,
            x * (prompt - density),
        ]
    };
    let (at_beta, at_rate) = (storage(beta), storage(rate));
    let gap = rate - beta;
    let apart: [f64; 5] = std::array::from_fn(|n| (at_beta[n] - at_rate[n]) / gap);
    [
        at_rate[0] + r * rate * apart[0],
        at_rate[1] + r * rate * apart[1],
        r * rate * (at_beta[2] + apart[0]) / gap,
        rate * apart[0],
        at_rate[2] + r * apart[0] - r * rate * (at_rate[2] + apart[0]) / gap,
        at_rate[3] + r * rate * apart[3],
        at_rate[4] + r * rate * apart[4],
    ]
}

/// Each bin's difference of `g` between its upper and lower edge, `g` taken
/// at the edge's delay after `arrival_us`.
fn bin_differences<const N: usize>(
    edges_us: &[f64],
    arrival_us: f64,
    g: impl Fn(f64) -> [f64; N],
) -> [Vec<f64>; N] {
    let at_edges: Vec<[f64; N]> = edges_us.iter().map(|edge| g(edge - arrival_us)).collect();
    std::array::from_fn(|n| at_edges.windows(2).map(|w| w[1][n] - w[0][n]).collect())
}

/// Gamma(3, rate=1) cumulative distribution at dimensionless time `x`.
fn gamma3_cdf(x: f64) -> f64 {
    if x <= 0.0 {
        return 0.0;
    }
    // Beyond x ≈ 750 the survival term exp(−x)·(1+x+x²/2) underflows to an
    // exact zero long before the polynomial can overflow (x² reaches
    // f64::INFINITY only at x ~ 1.3e154, where 0·inf would be NaN).
    if x >= 750.0 {
        return 1.0;
    }
    if x < 0.5 {
        // Integral of x² exp(-x) / 2 as an alternating series. The direct
        // `1 - exp(-x)(1+x+x²/2)` loses most digits for small x.
        let mut n = 0_u32;
        let mut term = x.powi(3) / 6.0;
        let mut sum = term;
        loop {
            let n_f = f64::from(n);
            term *= -x * (n_f + 3.0) / ((n_f + 1.0) * (n_f + 4.0));
            let next = sum + term;
            n += 1;
            if next == sum || n >= 128 {
                return next.clamp(0.0, 1.0);
            }
            sum = next;
        }
    }
    (1.0 - (-x).exp() * (1.0 + x + 0.5 * x * x)).clamp(0.0, 1.0)
}

/// Energy-dependence law for an Ikeda–Carpenter parameter.
///
/// A small closed set so the *fixed-or-fit* cases share one representation:
/// fixed ⇒ [`Const`](EnergyLaw::Const); fit ⇒ a parametric law whose
/// coefficients are the fit variables.
#[derive(Debug, Clone, PartialEq)]
pub enum EnergyLaw {
    /// Energy-independent constant.
    Const(f64),
    /// `a0·√(E[eV]) + a1` — leading epithermal scaling of the fast rate `α(E)`.
    SqrtE { a0: f64, a1: f64 },
    /// Mantid IC form `1/(a0 + a1·λ)`, λ (Å) = `LAMBDA_ANGSTROM_FACTOR`/√E.
    /// Behaves as `α ∝ √E` at low E and saturates to `1/a0` at high E.
    InverseLambda { a0: f64, a1: f64 },
    /// `exp(−E[meV]/kappa)` — storage fraction `R(E)`, → 0 in the eV regime.
    ExpMilliEv { kappa: f64 },
    /// Linear in `√E` from `low` at `span_ev.0` to `high` at `span_ev.1` eV,
    /// and between `low` and `high` at every energy.
    SqrtEBetween {
        low: f64,
        high: f64,
        span_ev: (f64, f64),
    },
}

impl EnergyLaw {
    /// Evaluate the law at `energy_ev` (eV). Non-positive energy yields the
    /// `E→0` limit where well-defined (and a clamped value otherwise).
    #[must_use]
    pub fn eval(&self, energy_ev: f64) -> f64 {
        let e = energy_ev.max(0.0);
        match *self {
            EnergyLaw::Const(c) => c,
            EnergyLaw::SqrtE { a0, a1 } => a0 * e.sqrt() + a1,
            EnergyLaw::InverseLambda { a0, a1 } => {
                let denom = inverse_lambda_denom(a0, a1, e);
                if denom.abs() < MIN_RATE {
                    1.0 / MIN_RATE
                } else {
                    1.0 / denom
                }
            }
            EnergyLaw::ExpMilliEv { kappa } => {
                if kappa.abs() < MIN_RATE {
                    0.0
                } else {
                    (-(e * 1000.0) / kappa).exp()
                }
            }
            EnergyLaw::SqrtEBetween {
                low,
                high,
                span_ev: (e_low, e_high),
            } => {
                let share = (e.sqrt() - e_low.sqrt()) / (e_high.sqrt() - e_low.sqrt());
                (low + (high - low) * share).clamp(low.min(high), low.max(high))
            }
        }
    }

    /// True when the law is numerically singular at `energy_ev`: the raw
    /// `InverseLambda` denominator lies inside the ±[`MIN_RATE`] window that
    /// [`EnergyLaw::eval`] floors away. The floor keeps a fit trial step from
    /// dividing by zero mid-optimization, but a *configured* law inside that
    /// window is a mathematically undefined rate, not a large one — without
    /// this check the floor converts an undefined (or tiny-negative)
    /// denominator into a plausible huge positive rate before the range
    /// validation below can see it.
    #[must_use]
    pub(crate) fn is_singular_at(&self, energy_ev: f64) -> bool {
        match *self {
            EnergyLaw::InverseLambda { a0, a1 } => {
                inverse_lambda_denom(a0, a1, energy_ev.max(0.0)).abs() < MIN_RATE
            }
            // κ = 0 is undefined and a tiny NEGATIVE κ is a divergent law
            // (exp(+E/|κ|)); eval's floor maps both to 0.0, which is the
            // legitimate κ → 0⁺ limit only for positive κ.
            EnergyLaw::ExpMilliEv { kappa } => (-MIN_RATE..=0.0).contains(&kappa),
            _ => false,
        }
    }
}

/// Raw `InverseLambda` denominator `a0 + a1·λ(E)` — shared by [`EnergyLaw::eval`]
/// and the singularity check so the two can never disagree on the window.
fn inverse_lambda_denom(a0: f64, a1: f64, e: f64) -> f64 {
    let lambda = if e > 0.0 {
        LAMBDA_ANGSTROM_FACTOR / e.sqrt()
    } else {
        f64::INFINITY
    };
    a0 + a1 * lambda
}

/// Parameters of the Ikeda–Carpenter resolution model.
#[derive(Debug, Clone, PartialEq)]
pub struct IkedaCarpenterParams {
    /// Fast (slowing-down) rate `α(E)`, 1/µs. Must evaluate to at least
    /// 1e-9.
    pub alpha: EnergyLaw,
    /// Slow (storage) rate `β(E)`, 1/µs. Must evaluate to at least 1e-9.
    pub beta: EnergyLaw,
    /// Storage mixing fraction `R(E)`, 0 ≤ R ≤ 1.
    pub r: EnergyLaw,
    /// Optional proton-burst Gaussian standard deviation (µs); refused in
    /// detector time.
    pub burst_sigma_us: Option<f64>,
    /// Optional chopper/channel triangle FWHM (µs).
    pub channel_fwhm_us: Option<f64>,
    /// Optional detector tail; detector time only, refused by the kernel
    /// table.
    pub detector_tail: Option<DetectorTail>,
}

/// A detector that records a fraction `f` of the counts late by an
/// exponential delay of mean `τ_D`, and the rest at their arrival.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DetectorTail {
    /// `f`, within 0–1.
    pub fraction: f64,
    /// `τ_D` in µs, from 1e-9 to 1e9.
    pub decay_us: f64,
}

impl IkedaCarpenterParams {
    /// A pure-moderator parameter set (no burst, no channel) with constant
    /// rates — the simplest case, useful for fixed-parameter or unit-test use.
    #[must_use]
    pub fn constant(alpha: f64, beta: f64, r: f64) -> Self {
        Self {
            alpha: EnergyLaw::Const(alpha),
            beta: EnergyLaw::Const(beta),
            r: EnergyLaw::Const(r),
            burst_sigma_us: None,
            channel_fwhm_us: None,
            detector_tail: None,
        }
    }

    fn tail_rate(&self) -> Result<Option<(f64, f64)>, ResolutionParseError> {
        let Some(tail) = self.detector_tail else {
            return Ok(None);
        };
        if !((0.0..=1.0).contains(&tail.fraction)
            && DETECTOR_TAIL_DECAY_US.contains(&tail.decay_us))
        {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "a detector tail needs a fraction within 0–1 and a decay from {:e} to {:e} µs, \
                 got {tail:?}",
                DETECTOR_TAIL_DECAY_US.start(),
                DETECTOR_TAIL_DECAY_US.end()
            )));
        }
        Ok(Some((tail.fraction, 1.0 / tail.decay_us)))
    }

    fn recorded_at(&self, energy_ev: f64) -> Result<Recorded, ResolutionParseError> {
        self.validate_probe_energy(energy_ev)?;
        Ok(Recorded {
            alpha: self.alpha.eval(energy_ev),
            beta: self.beta.eval(energy_ev),
            r: self.r.eval(energy_ev),
            tail: self.tail_rate()?,
        })
    }

    fn triangle_fwhm_us(&self) -> Result<f64, ResolutionParseError> {
        if self.burst_sigma_us.unwrap_or(0.0) != 0.0 {
            return Err(ResolutionParseError::InvalidFormat(
                "a Gaussian burst is not supported in detector time; model the proton \
                 pulse with channel_fwhm_us"
                    .to_string(),
            ));
        }
        let fwhm_us = self.channel_fwhm_us.unwrap_or(0.0);
        if !(fwhm_us.is_finite() && fwhm_us >= 0.0) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "the triangle's FWHM must be finite and 0 or more, got {fwhm_us} µs"
            )));
        }
        Ok(fwhm_us)
    }

    /// The first and last delay, in µs after the nominal arrival, of a
    /// neutron of `energy_ev`, outside which [`IkedaCarpenter::detector_bin_probabilities`]
    /// gives it less than [`NEGLIGIBLE_ARRIVAL_PROBABILITY`] chance of
    /// arriving: the triangle's half-base before 0, and that half-base after
    /// the delay where the chance of a later record without the triangle
    /// falls to that chance.
    ///
    /// # Errors
    /// [`ResolutionParseError::InvalidFormat`] when `energy_ev` is not
    /// positive and finite, a law is singular or out of range there, the
    /// pulse has a Gaussian burst, or its detector tail is out of range.
    pub fn delays_us(&self, energy_ev: f64) -> Result<(f64, f64), ResolutionParseError> {
        let pulse = self.recorded_at(energy_ev)?;
        let h = self.triangle_fwhm_us()?;
        Ok((-h, tail_delay(&pulse) + h))
    }

    /// Probability that a neutron of `true_energy_ev` whose nominal arrival,
    /// before the pulse's delay, is `arrival_us` is recorded in each bin of
    /// `detector_time_edges_us`; as [`IkedaCarpenter::detector_bin_probabilities`].
    ///
    /// # Errors
    /// As [`IkedaCarpenter::detector_bin_probabilities`], with `arrival_us` finite.
    pub fn bin_probabilities_at(
        &self,
        true_energy_ev: f64,
        arrival_us: f64,
        detector_time_edges_us: &[f64],
    ) -> Result<Vec<f64>, ResolutionParseError> {
        let (pulse, h) = self.at_arrival(true_energy_ev, arrival_us, detector_time_edges_us)?;
        let cdf = |x: f64| {
            if h == 0.0 {
                pulse.cdf(x)
            } else {
                triangle_fold(pulse.alpha, h, x, |tau| [pulse.cdf(tau)])[0].clamp(0.0, 1.0)
            }
        };
        let cdfs: Vec<f64> = detector_time_edges_us
            .iter()
            .map(|edge| cdf(edge - arrival_us))
            .collect();
        Ok(cdfs.windows(2).map(|w| (w[1] - w[0]).max(0.0)).collect())
    }

    /// The derivative of each of [`Self::bin_probabilities_at`] with respect
    /// to `arrival_us`, per µs: the pulse density, folded with the triangle,
    /// at the bin's lower edge minus that at its upper edge.
    ///
    /// # Errors
    /// As [`Self::bin_probabilities_at`].
    pub fn bin_arrival_slopes_at(
        &self,
        true_energy_ev: f64,
        arrival_us: f64,
        detector_time_edges_us: &[f64],
    ) -> Result<Vec<f64>, ResolutionParseError> {
        let (pulse, h) = self.at_arrival(true_energy_ev, arrival_us, detector_time_edges_us)?;
        let density = |x: f64| {
            if h == 0.0 {
                pulse.density(x)
            } else {
                triangle_fold(pulse.alpha, h, x, |tau| [pulse.density(tau)])[0]
            }
        };
        let densities: Vec<f64> = detector_time_edges_us
            .iter()
            .map(|edge| density(edge - arrival_us))
            .collect();
        Ok(densities.windows(2).map(|w| w[0] - w[1]).collect())
    }

    /// The derivatives of each of [`Self::bin_probabilities_at`] with
    /// respect to `α`, `β`, `R` and the square of the triangle's FWHM `h`,
    /// each evaluated at `true_energy_ev`: `[∂P/∂α, ∂P/∂β, ∂P/∂R, ∂P/∂h²]`,
    /// one value per bin.  The probabilities are even in `h` and smooth in
    /// `h²`, so `∂P/∂h²` is finite and generally not 0 at `h = 0`.
    ///
    /// # Errors
    /// As [`Self::bin_probabilities_at`].
    pub fn bin_pulse_slopes_at(
        &self,
        true_energy_ev: f64,
        arrival_us: f64,
        detector_time_edges_us: &[f64],
    ) -> Result<[Vec<f64>; 4], ResolutionParseError> {
        let (pulse, h) = self.at_arrival(true_energy_ev, arrival_us, detector_time_edges_us)?;
        let slopes = |x: f64| {
            if h == 0.0 {
                let [a, b, c, d] = pulse.slopes(x);
                [a, b, c, d / 12.0]
            } else {
                folded_slopes(&pulse, h, x)
            }
        };
        Ok(bin_differences(detector_time_edges_us, arrival_us, slopes))
    }

    /// The derivatives of each of [`Self::bin_probabilities_at`] with
    /// respect to the detector tail's `f` and `τ_D`, each per its unit:
    /// `[∂P/∂f, ∂P/∂τ_D]`, one value per bin.
    ///
    /// # Errors
    /// As [`Self::bin_probabilities_at`], and
    /// [`ResolutionParseError::InvalidFormat`] for a pulse without a
    /// detector tail.
    pub fn bin_tail_slopes_at(
        &self,
        true_energy_ev: f64,
        arrival_us: f64,
        detector_time_edges_us: &[f64],
    ) -> Result<[Vec<f64>; 2], ResolutionParseError> {
        let (pulse, h) = self.at_arrival(true_energy_ev, arrival_us, detector_time_edges_us)?;
        if pulse.tail.is_none() {
            return Err(ResolutionParseError::InvalidFormat(
                "the pulse has no detector tail".to_string(),
            ));
        }
        let slopes = |x: f64| {
            if h == 0.0 {
                pulse.tail_slopes(x)
            } else {
                triangle_fold(pulse.alpha, h, x, |tau| pulse.tail_slopes(tau))
            }
        };
        Ok(bin_differences(detector_time_edges_us, arrival_us, slopes))
    }

    fn at_arrival(
        &self,
        true_energy_ev: f64,
        arrival_us: f64,
        detector_time_edges_us: &[f64],
    ) -> Result<(Recorded, f64), ResolutionParseError> {
        let pulse = self.recorded_at(true_energy_ev)?;
        if !arrival_us.is_finite() {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "the nominal arrival must be finite, got {arrival_us} µs"
            )));
        }
        if detector_time_edges_us.len() < 2
            || detector_time_edges_us.iter().any(|x| !x.is_finite())
            || detector_time_edges_us.windows(2).any(|w| w[0] >= w[1])
        {
            return Err(ResolutionParseError::InvalidFormat(
                "detector time edges must contain at least two finite, strictly increasing values"
                    .to_string(),
            ));
        }
        Ok((pulse, self.triangle_fwhm_us()?))
    }

    fn validate_probe_energy(&self, energy_ev: f64) -> Result<(), ResolutionParseError> {
        if !energy_ev.is_finite() || energy_ev <= 0.0 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "true energy must be positive and finite, got {energy_ev}"
            )));
        }
        for (name, law) in [("alpha", &self.alpha), ("beta", &self.beta), ("R", &self.r)] {
            if law.is_singular_at(energy_ev) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Ikeda–Carpenter {name}({energy_ev}) law is singular (an \
                     InverseLambda denominator within ±{MIN_RATE} of zero, \
                     or an ExpMilliEv κ in [−{MIN_RATE}, 0])"
                )));
            }
        }
        let alpha = self.alpha.eval(energy_ev);
        let beta = self.beta.eval(energy_ev);
        let r = self.r.eval(energy_ev);
        if !alpha.is_finite() || alpha < MIN_RATE {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter alpha({energy_ev}) must be finite and at least {MIN_RATE} \
                 µs⁻¹, got {alpha}"
            )));
        }
        if !beta.is_finite() || beta < MIN_RATE {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter beta({energy_ev}) must be finite and at least {MIN_RATE} \
                 µs⁻¹, got {beta}"
            )));
        }
        if !r.is_finite() || !(0.0..=1.0).contains(&r) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter R({energy_ev}) must be in [0, 1], got {r}"
            )));
        }
        Ok(())
    }
}

/// Configuration of the energy / time grids used to synthesize the IC kernel
/// table. The energy grid spans the data range densely enough that interref
/// interpolation error is negligible.
#[derive(Debug, Clone)]
pub struct SynthesisGrid {
    /// Lowest reference energy (eV), > 0.
    pub e_min_ev: f64,
    /// Highest reference energy (eV), > `e_min_ev`.
    pub e_max_ev: f64,
    /// Number of log-spaced reference energies (≥ 2).
    pub n_energies: usize,
    /// Number of τ-samples spanning the prompt core of each kernel
    /// (≥ [`MIN_N_TAU`](crate::ikeda_carpenter) = 8). A long storage tail
    /// (β ≪ α with R > 0) grows the per-kernel sample count beyond `n_tau`;
    /// that count is capped at 8192 (`MAX_TAU_SAMPLES`), past which the
    /// τ-step widens — never below the resolution floor (prompt core at the
    /// `n_tau = 8` density, folds at ≥ 3 triangle samples per side / ≥ 1
    /// sample per burst σ): a combination that cannot be resolved within the
    /// cap is rejected by [`IkedaCarpenter::new`]. Active burst/channel folds
    /// add their ±(`GAUSS_REACH_SIGMAS`·σ + FWHM) margin samples on top of
    /// the cap.
    pub n_tau: usize,
}

impl SynthesisGrid {
    /// A sensible default grid: `[e_min, e_max]` log-spaced over
    /// [`DEFAULT_N_ENERGIES`] points with [`DEFAULT_N_TAU`] τ-samples.
    #[must_use]
    pub fn new(e_min_ev: f64, e_max_ev: f64) -> Self {
        Self {
            e_min_ev,
            e_max_ev,
            n_energies: DEFAULT_N_ENERGIES,
            n_tau: DEFAULT_N_TAU,
        }
    }
}

/// Analytical Ikeda–Carpenter resolution model.
///
/// Synthesizes a [`TabulatedResolution`] at construction and applies it through
/// the same broadening path as a Monte-Carlo file. Cloning is cheap (the
/// synthesized table is the only large field). Re-synthesize (construct anew)
/// when fitting changes the parameters.
#[derive(Debug, Clone)]
pub struct IkedaCarpenter {
    params: IkedaCarpenterParams,
    flight_path_m: f64,
    /// Shared: rebinding the flight path must not copy it.
    ref_energies: Arc<Vec<f64>>,
    n_tau: usize,
    tabulated: TabulatedResolution,
}

impl IkedaCarpenter {
    /// Build the model, synthesizing its kernel table over `grid`.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] for a non-positive
    /// flight path, a degenerate grid (`n_energies < 2`, `n_tau < 8`,
    /// `e_min ≤ 0`, `e_max ≤ e_min`), an `α(E)` or `β(E)` below 1e-9 µs⁻¹ at a
    /// reference energy, a parameter/grid
    /// combination whose τ-grid cannot resolve the prompt core and requested
    /// folds within the `MAX_TAU_SAMPLES` cap at some reference energy (see
    /// `tau_geometry` — remedy: larger `β`, `R = 0`, or a wider/disabled
    /// fold), a detector tail, or if the synthesized kernels fail
    /// [`TabulatedResolution::from_kernels`] validation.
    pub fn new(
        params: IkedaCarpenterParams,
        flight_path_m: f64,
        grid: &SynthesisGrid,
    ) -> Result<Self, ResolutionParseError> {
        if !flight_path_m.is_finite() || flight_path_m <= 0.0 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Flight path must be a positive finite number, got {flight_path_m}"
            )));
        }
        if grid.n_energies < 2 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "n_energies must be >= 2, got {}",
                grid.n_energies
            )));
        }
        if grid.n_tau < MIN_N_TAU {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "n_tau must be >= {MIN_N_TAU}, got {}",
                grid.n_tau
            )));
        }
        if !grid.e_min_ev.is_finite()
            || !grid.e_max_ev.is_finite()
            || grid.e_min_ev <= 0.0
            || grid.e_max_ev <= grid.e_min_ev
        {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Require finite 0 < e_min < e_max, got [{}, {}]",
                grid.e_min_ev, grid.e_max_ev
            )));
        }
        let ln_lo = grid.e_min_ev.ln();
        let ln_hi = grid.e_max_ev.ln();
        let denom = (grid.n_energies - 1) as f64;
        let ref_energies: Vec<f64> = (0..grid.n_energies)
            .map(|i| (ln_lo + (i as f64 / denom) * (ln_hi - ln_lo)).exp())
            .collect();

        // Reject singular rate laws before the range checks: a near-zero
        // InverseLambda denominator is an undefined configuration, and eval's
        // numerical floor would otherwise convert it into a plausible huge
        // positive rate that the α > 0 / β > 0 checks below cannot distinguish
        // from genuine physics.
        for (name, law) in [("α", &params.alpha), ("β", &params.beta), ("R", &params.r)] {
            if let Some(&bad) = ref_energies.iter().find(|&&e| law.is_singular_at(e)) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Ikeda–Carpenter {name}(E) law is singular at E = {bad} eV \
                     (an InverseLambda denominator within \
                     ±{MIN_RATE} of zero, or an ExpMilliEv κ in [−{MIN_RATE}, 0])"
                )));
            }
        }
        if let Some(&bad) = ref_energies.iter().find(|&&e| {
            let a = params.alpha.eval(e);
            !a.is_finite() || a < MIN_RATE
        }) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter α(E) must be at least {MIN_RATE} µs⁻¹, but α({bad}) = {} is not",
                params.alpha.eval(bad)
            )));
        }
        // β is a rate, so every synthesized reference energy must give a
        // positive finite value. Reject invalid laws instead of clamping them
        // into a different pulse.
        if let Some(&bad) = ref_energies.iter().find(|&&e| {
            let beta = params.beta.eval(e);
            !beta.is_finite() || beta < MIN_RATE
        }) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter β(E) must be at least {MIN_RATE} µs⁻¹, but β({bad}) = {} is not",
                params.beta.eval(bad)
            )));
        }
        // Reject storage-fraction laws that fall outside [0, 1] (synthesis clamps
        // R, masking a mis-specified law); a physical mixing fraction is in [0,1].
        if let Some(&bad) = ref_energies.iter().find(|&&e| {
            let r = params.r.eval(e);
            !r.is_finite() || !(0.0..=1.0).contains(&r)
        }) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Ikeda–Carpenter R(E) must be in [0, 1], but R({bad}) = {} is not",
                params.r.eval(bad)
            )));
        }
        // Reject invalid optional instrument-convolution widths up front; synthesis
        // otherwise masks a negative width via `.abs()` and silently swallows a NaN.
        for (name, v) in [
            ("burst_sigma_us", params.burst_sigma_us),
            ("channel_fwhm_us", params.channel_fwhm_us),
        ] {
            if let Some(x) = v
                && (!x.is_finite() || x < 0.0)
            {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Ikeda–Carpenter {name} must be finite and >= 0, got {x}"
                )));
            }
        }
        if params.detector_tail.is_some() {
            return Err(ResolutionParseError::InvalidFormat(
                "a detector tail is not supported in the kernel table".to_string(),
            ));
        }

        // Synthesis is fallible: a kernel whose τ-grid cannot resolve the
        // requested physics within MAX_TAU_SAMPLES (long slow tail vs a fine
        // fold / fast prompt core) errs loudly here instead of silently
        // degrading — see `tau_geometry`.
        let kernels: Vec<(Vec<f64>, Vec<f64>)> = ref_energies
            .iter()
            .map(|&e| synth_kernel(&params, grid.n_tau, e))
            .collect::<Result<_, _>>()?;

        let tabulated =
            TabulatedResolution::from_kernels(ref_energies.clone(), kernels, flight_path_m)?;

        Ok(Self {
            params,
            flight_path_m,
            ref_energies: Arc::new(ref_energies),
            n_tau: grid.n_tau,
            tabulated,
        })
    }

    /// The synthesized tabulated kernel set (the broadening engine).
    #[must_use]
    pub fn tabulated(&self) -> &TabulatedResolution {
        &self.tabulated
    }

    /// The IC parameters.
    #[must_use]
    pub fn params(&self) -> &IkedaCarpenterParams {
        &self.params
    }

    /// Flight-path length (m).
    #[must_use]
    pub fn flight_path_m(&self) -> f64 {
        self.flight_path_m
    }

    /// The same pulse read against a different flight path.
    ///
    /// The synthesized kernels are emission-time distributions built from
    /// α(E), β(E) and R(E) — no flight path appears in the synthesis. The
    /// flight path enters only the TOF↔energy map and the nominal arrival
    /// time, so rebinding it is exact and does not resynthesize the table.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] if `flight_path_m` is
    /// not positive and finite.
    pub fn with_flight_path(&self, flight_path_m: f64) -> Result<Self, ResolutionParseError> {
        Ok(Self {
            params: self.params.clone(),
            flight_path_m,
            ref_energies: Arc::clone(&self.ref_energies),
            n_tau: self.n_tau,
            tabulated: self.tabulated.with_flight_path(flight_path_m)?,
        })
    }

    /// Reference energies (eV, ascending) the table was synthesized on.
    #[must_use]
    pub fn ref_energies(&self) -> &[f64] {
        &self.ref_energies
    }

    /// The pulse as the counts path reads it, over the energies the table was
    /// synthesized on.
    #[must_use]
    pub fn detector_pulse(&self) -> DetectorPulse {
        DetectorPulse {
            params: self.params.clone(),
            energy_span_ev: (
                self.ref_energies[0],
                self.ref_energies[self.ref_energies.len() - 1],
            ),
            n_tau: self.n_tau,
        }
    }

    /// Evaluate the (burst/channel-folded) IC kernel at one energy.
    ///
    /// Returns ascending TOF-offsets (µs, 0 = pulse start) and peak-normalized
    /// weights (max = 1), matching the [`TabulatedResolution`] storage
    /// convention.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] when the requested
    /// energy or an energy-dependent rate/fraction is non-physical, or when
    /// the τ-grid cannot resolve the prompt core and requested folds within
    /// `MAX_TAU_SAMPLES` at this energy. Construction validates every
    /// *reference* energy, but a probe outside `[e_min, e_max]` can still leave
    /// the physical or resolvable region.
    pub fn kernel_at(&self, energy_ev: f64) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
        self.params.validate_probe_energy(energy_ev)?;
        synth_kernel(&self.params, self.n_tau, energy_ev)
    }

    /// Evaluate the physical source pulse at one true neutron energy.
    ///
    /// Returns sampled moderator-delay coordinates in µs and peak-normalized
    /// densities. Unlike [`Self::kernel_at`], this method does not move the
    /// pulse mode to zero. With no symmetric proton/channel fold, the delay is
    /// causal and starts at zero. A symmetric fold may extend the sampled
    /// support below zero relative to its stated time origin.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] if `energy_ev` is not a
    /// positive finite true energy, if an energy law is unphysical at that
    /// energy, or if the requested pulse cannot be resolved by the configured
    /// sampling limits.
    pub fn source_pulse_at(
        &self,
        energy_ev: f64,
    ) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
        self.params.validate_probe_energy(energy_ev)?;
        synth_source_pulse(&self.params, self.n_tau, energy_ev)
    }

    /// Probability that a neutron of known true energy is recorded in each
    /// supplied detector-time bin.
    ///
    /// `detector_time_edges_us` are the actual measured bin edges. The nominal
    /// arrival time is
    /// `timing_offset_us + TOF_FACTOR * flight_path_m / sqrt(true_energy_ev)`.
    /// `timing_offset_us` represents the shared clock/detector offset; it does
    /// not absorb or remove the moderator pulse's physical mode.
    ///
    /// The returned vector has one entry per adjacent edge pair and is not
    /// renormalized to the supplied window: bins that do not cover the full
    /// pulse correctly sum to less than one.
    ///
    /// Each probability is a difference of the analytical IC CDF, folded with
    /// the proton-pulse triangle when `channel_fwhm_us` is set.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] unless the true energy
    /// is physical, the timing offset is finite, at least two finite bin
    /// edges are supplied in strictly increasing order, and the pulse has no
    /// Gaussian burst.
    pub fn detector_bin_probabilities(
        &self,
        true_energy_ev: f64,
        detector_time_edges_us: &[f64],
        timing_offset_us: f64,
    ) -> Result<Vec<f64>, ResolutionParseError> {
        if !timing_offset_us.is_finite() {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "timing_offset_us must be finite, got {timing_offset_us}"
            )));
        }
        self.params.bin_probabilities_at(
            true_energy_ev,
            timing_offset_us + TOF_FACTOR * self.flight_path_m / true_energy_ev.sqrt(),
            detector_time_edges_us,
        )
    }
}

/// An Ikeda–Carpenter pulse in detector time: its laws, the energies they hold
/// over, and the number of samples across its prompt core used to find its
/// rise.
#[derive(Debug, Clone)]
pub struct DetectorPulse {
    params: IkedaCarpenterParams,
    energy_span_ev: (f64, f64),
    n_tau: usize,
}

impl DetectorPulse {
    /// The pulse with laws `params` over `energy_span_ev`, its rise found on
    /// `n_tau` samples across the prompt core.
    ///
    /// # Errors
    /// [`ResolutionParseError::InvalidFormat`] unless `0 < low < high` are
    /// finite, `n_tau` is at least 8, the triangle's FWHM is finite and 0 or
    /// more, the pulse has no Gaussian burst, a detector tail has its
    /// fraction within 0–1 and its decay from 1e-9 to 1e9 µs, and the
    /// laws are valid at both ends of the span: `α` and `β` finite and at
    /// least 1e-9 µs⁻¹, `R` within 0–1.
    pub fn new(
        params: IkedaCarpenterParams,
        energy_span_ev: (f64, f64),
        n_tau: usize,
    ) -> Result<Self, ResolutionParseError> {
        params.triangle_fwhm_us()?;
        params.tail_rate()?;
        let (low, high) = energy_span_ev;
        if !(low.is_finite() && high.is_finite() && 0.0 < low && low < high) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "the energy span must satisfy 0 < low < high, got {energy_span_ev:?}"
            )));
        }
        params.validate_probe_energy(low)?;
        params.validate_probe_energy(high)?;
        if n_tau < MIN_N_TAU {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "n_tau must be >= {MIN_N_TAU}, got {n_tau}"
            )));
        }
        Ok(Self {
            params,
            energy_span_ev,
            n_tau,
        })
    }

    /// The pulse's laws.
    #[must_use]
    pub fn params(&self) -> &IkedaCarpenterParams {
        &self.params
    }

    /// `(low, high)`, the energies in eV the laws hold over.
    #[must_use]
    pub fn energy_span_ev(&self) -> (f64, f64) {
        self.energy_span_ev
    }

    /// The samples across the prompt core used to find the pulse's rise.
    #[must_use]
    pub fn n_tau(&self) -> usize {
        self.n_tau
    }

    /// The time in µs from the first sample of the pulse at `energy_ev` to
    /// its peak.
    ///
    /// # Errors
    /// [`ResolutionParseError::InvalidFormat`] when `energy_ev` is not
    /// positive and finite, a law is singular or out of range there, or the
    /// pulse has a Gaussian burst.
    pub fn rise_us(&self, energy_ev: f64) -> Result<f64, ResolutionParseError> {
        self.params.validate_probe_energy(energy_ev)?;
        self.params.triangle_fwhm_us()?;
        sampled_rise_us(&self.params, self.n_tau, energy_ev)
    }
}

/// The chance of a later arrival below which an Ikeda–Carpenter pulse's
/// tail is treated as ended.
pub const NEGLIGIBLE_ARRIVAL_PROBABILITY: f64 = 1e-7;

fn tail_delay(pulse: &Recorded) -> f64 {
    let later = |tau: f64| 1.0 - pulse.cdf(tau);
    let late = pulse.delayed().map_or(f64::INFINITY, |(_, rate)| rate);
    let mut high = 1.0 / pulse.alpha.min(pulse.beta).min(late).max(MIN_RATE);
    while later(high) > NEGLIGIBLE_ARRIVAL_PROBABILITY {
        high *= 2.0;
    }
    let mut low = 0.0;
    while high - low > f64::EPSILON * high {
        let middle = 0.5 * (low + high);
        if later(middle) > NEGLIGIBLE_ARRIVAL_PROBABILITY {
            low = middle;
        } else {
            high = middle;
        }
    }
    high
}

fn requested_step(params: &IkedaCarpenterParams, n_tau: usize, alpha: f64) -> (f64, f64, String) {
    let fast_reach = FAST_REACH_E_FOLDS / alpha;
    let mut dtau_req = fast_reach / (n_tau as f64 - 1.0);
    let mut floor = fast_reach / (MIN_N_TAU as f64 - 1.0);
    let mut fold_desc = String::new();
    if let Some(fwhm) = params.channel_fwhm_us
        && fwhm > 0.0
    {
        let tri_floor = fwhm / TRI_MIN_SAMPLES_PER_SIDE;
        dtau_req = dtau_req.min(tri_floor);
        floor = floor.min(tri_floor);
        fold_desc.push_str(&format!(", channel triangle FWHM = {fwhm} µs"));
    }
    if let Some(sigma) = params.burst_sigma_us
        && sigma > 0.0
    {
        dtau_req = dtau_req.min(sigma);
        floor = floor.min(sigma);
        fold_desc.push_str(&format!(", burst σ = {sigma} µs"));
    }

    (dtau_req, floor, fold_desc)
}

/// τ-grid geometry for one kernel: `(dtau, tau_max, margin)`, or a
/// descriptive error when no exact sampled representation fits the cap.
///
/// The step is anchored to the PROMPT core — `n_tau` samples across the fast
/// Gamma(3) pulse (`fast_reach / (n_tau − 1)`) — and REFINED to resolve any
/// requested instrument fold (triangle: ≥ [`TRI_MIN_SAMPLES_PER_SIDE`]
/// samples per side, i.e. `dtau ≤ FWHM/TRI_MIN_SAMPLES_PER_SIDE`; Gaussian burst: `dtau ≤ σ`, i.e.
/// ≥ [`GAUSS_REACH_SIGMAS`] samples per ±`GAUSS_REACH_SIGMAS`·σ side). A longer storage tail
/// (β ≪ α, R > 0) extends the SAMPLE COUNT (`j_hi ∝ tau_max/dtau`) instead
/// of the step; [`MAX_TAU_SAMPLES`] bounds that count by widening the step —
/// but never past the resolution FLOOR (`fast_reach / (MIN_N_TAU − 1)` for
/// the prompt core, the fold minima above for folds). A combination whose
/// floor cannot be met within the cap has no faithful sampled representation
/// here, so it is rejected loudly: a capped step above the fold width would
/// degenerate the sampled triangle to an exact delta `[0,1,0]` (the fold
/// silently vanishes), and a capped step above the prompt scale steps OVER
/// the prompt pulse entirely (probe: α = 250, β = 0.02, R = 0.1 loses the
/// prompt's 0.9 weight share).
///
/// Equals `max(fast_reach/(n_tau−1), tau_max/(MAX_TAU_SAMPLES−1))` whenever
/// no fold is finer than the prompt design step and the capped step stays at
/// or below the floor.
fn tau_geometry(
    params: &IkedaCarpenterParams,
    n_tau: usize,
    alpha: f64,
    beta: f64,
    r: f64,
) -> Result<(f64, f64, f64), String> {
    let tau_max = tau_reach(alpha, beta, r);
    let (dtau_req, floor, fold_desc) = requested_step(params, n_tau, alpha);

    let capped_step = tau_max / (MAX_TAU_SAMPLES as f64 - 1.0);
    if capped_step > floor {
        return Err(format!(
            "the {MAX_TAU_SAMPLES}-sample τ-grid cap cannot resolve the requested physics: \
             the pulse spans τ_max = {tau_max:.3} µs (α = {alpha:.4} µs⁻¹, β = {beta:.4} µs⁻¹, \
             R = {r:.3}{fold_desc}), forcing a τ-step of {capped_step:.4} µs — above the \
             finest-feature resolution floor of {floor:.4} µs. Increase β (shorter storage \
             tail), set R = 0 (drop the storage term), or widen/disable the burst/channel fold"
        ));
    }
    Ok((dtau_req.max(capped_step), tau_max, margin_of(params)))
}

fn tau_reach(alpha: f64, beta: f64, r: f64) -> f64 {
    let fast_reach = FAST_REACH_E_FOLDS / alpha;
    if r > R_NEGLIGIBLE {
        fast_reach.max(SLOW_REACH_E_FOLDS / beta)
    } else {
        fast_reach
    }
}

/// Symmetric τ-grid margin for the burst/channel folds:
/// ±([`GAUSS_REACH_SIGMAS`]·σ + FWHM), the folds' full reach. These samples
/// come ON TOP of [`MAX_TAU_SAMPLES`] (the cap governs the pulse body only).
fn margin_of(params: &IkedaCarpenterParams) -> f64 {
    params
        .burst_sigma_us
        .map_or(0.0, |s| GAUSS_REACH_SIGMAS * s)
        + params.channel_fwhm_us.unwrap_or(0.0)
}

/// Synthesize one `(offsets, weights)` kernel for `energy_ev` from the IC
/// parameters: sample `I(τ)`, fold in burst + channel, trim negligible tails,
/// peak-normalize.
///
/// Offsets are on the emission clock: `0` is when the pulse begins, which is
/// where the Ikeda–Carpenter function puts it. A loaded SAMMY UDR file is
/// anchored differently and never passes through here.
///
/// # Errors
/// [`ResolutionParseError::InvalidFormat`] when [`tau_geometry`] cannot
/// resolve the prompt core and requested folds within [`MAX_TAU_SAMPLES`].
fn synth_kernel(
    params: &IkedaCarpenterParams,
    n_tau: usize,
    energy_ev: f64,
) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
    synth_source_pulse(params, n_tau, energy_ev)
}

/// Synthesize one physical-time source pulse without moving its mode.
fn synth_source_pulse_density(
    params: &IkedaCarpenterParams,
    n_tau: usize,
    energy_ev: f64,
) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
    let alpha = params.alpha.eval(energy_ev).max(MIN_RATE);
    let beta = params.beta.eval(energy_ev).max(MIN_RATE);
    let r = params.r.eval(energy_ev).clamp(0.0, 1.0);

    let geometry = tau_geometry(params, n_tau, alpha, beta, r).map_err(|msg| {
        ResolutionParseError::InvalidFormat(format!(
            "Ikeda–Carpenter kernel at E = {energy_ev} eV: {msg}"
        ))
    })?;
    sampled_density(params, (alpha, beta, r), geometry, energy_ev)
}

fn sampled_rise_us(
    params: &IkedaCarpenterParams,
    n_tau: usize,
    energy_ev: f64,
) -> Result<f64, ResolutionParseError> {
    let alpha = params.alpha.eval(energy_ev).max(MIN_RATE);
    let beta = params.beta.eval(energy_ev).max(MIN_RATE);
    let r = params.r.eval(energy_ev).clamp(0.0, 1.0);
    let fast_reach = FAST_REACH_E_FOLDS / alpha;
    let margin = margin_of(params);
    let (dtau_req, _, _) = requested_step(params, n_tau, alpha);
    let dtau = dtau_req.max((fast_reach + 2.0 * margin) / (MAX_TAU_SAMPLES as f64 - 1.0));
    let (times, densities) = sampled_density(
        params,
        (alpha, beta, r),
        (dtau, fast_reach, margin),
        energy_ev,
    )?;
    Ok(times[argmax(&densities)] - times[0])
}

fn sampled_density(
    params: &IkedaCarpenterParams,
    (alpha, beta, r): (f64, f64, f64),
    (dtau, tau_max, margin): (f64, f64, f64),
    energy_ev: f64,
) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
    // Extend the grid to slightly negative τ so a symmetric burst/channel can
    // spread the leading edge correctly (the moderator pulse itself is 0 there).
    let j_lo: isize = -((margin / dtau).ceil() as isize);
    let j_hi: isize = ((tau_max + margin) / dtau).ceil() as isize;

    let taus: Vec<f64> = (j_lo..=j_hi).map(|j| j as f64 * dtau).collect();
    let mut weights: Vec<f64> = taus.iter().map(|&t| ic_pulse(alpha, beta, r, t)).collect();

    let sampled_area = dtau
        * (0.5 * weights[0]
            + weights[1..weights.len() - 1].iter().sum::<f64>()
            + 0.5 * weights[weights.len() - 1]);
    if !sampled_area.is_finite() || sampled_area <= 0.0 {
        return Err(ResolutionParseError::InvalidFormat(format!(
            "Ikeda–Carpenter pulse at E = {energy_ev} eV has zero sampled area"
        )));
    }

    if let Some(sigma) = params.burst_sigma_us
        && sigma > 0.0
    {
        weights = convolve_same(&weights, &gaussian_kernel(dtau, sigma));
    }
    if let Some(fwhm) = params.channel_fwhm_us
        && fwhm > 0.0
    {
        weights = convolve_same(&weights, &triangle_kernel(dtau, fwhm));
    }

    let peak_idx = argmax(&weights);
    let peak_val = weights[peak_idx].max(f64::MIN_POSITIVE);

    // Trim tails below TRIM_REL of peak, keeping one guard sample each side so
    // the convolution's neighbor-difference trapezoid widths stay defined.
    let thresh = TRIM_REL * peak_val;
    let lo = weights
        .iter()
        .position(|&w| w > thresh)
        .map_or(0, |i| i.saturating_sub(1));
    let hi = weights
        .iter()
        .rposition(|&w| w > thresh)
        .map_or(weights.len() - 1, |i| (i + 1).min(weights.len() - 1));

    let offsets: Vec<f64> = (lo..=hi).map(|j| taus[j]).collect();
    let densities: Vec<f64> = (lo..=hi).map(|j| weights[j]).collect();
    Ok((offsets, densities))
}

/// Synthesize the public peak-normalized source-pulse representation.
fn synth_source_pulse(
    params: &IkedaCarpenterParams,
    n_tau: usize,
    energy_ev: f64,
) -> Result<(Vec<f64>, Vec<f64>), ResolutionParseError> {
    let (offsets, densities) = synth_source_pulse_density(params, n_tau, energy_ev)?;
    let peak = densities
        .iter()
        .copied()
        .fold(0.0_f64, f64::max)
        .max(f64::MIN_POSITIVE);
    let weights = densities.into_iter().map(|value| value / peak).collect();
    Ok((offsets, weights))
}

/// Index of the maximum element (first on ties). Slice is non-empty by
/// construction in [`synth_kernel`].
fn argmax(xs: &[f64]) -> usize {
    let mut best = 0;
    let mut best_v = xs[0];
    for (i, &x) in xs.iter().enumerate().skip(1) {
        if x > best_v {
            best_v = x;
            best = i;
        }
    }
    best
}

/// Symmetric, unit-sum Gaussian kernel sampled on a `dtau`-spaced grid out to
/// ±[`GAUSS_REACH_SIGMAS`]·σ.
fn gaussian_kernel(dtau: f64, sigma: f64) -> Vec<f64> {
    let half = ((GAUSS_REACH_SIGMAS * sigma / dtau).ceil() as isize).max(1);
    let mut k: Vec<f64> = (-half..=half)
        .map(|j| {
            let t = j as f64 * dtau / sigma;
            (-0.5 * t * t).exp()
        })
        .collect();
    normalize_sum(&mut k);
    k
}

/// Symmetric, unit-sum triangle kernel of FWHM `fwhm` (half-base = FWHM;
/// variance FWHM²/6) sampled on a `dtau`-spaced grid.
fn triangle_kernel(dtau: f64, fwhm: f64) -> Vec<f64> {
    let a = fwhm; // half-base equals FWHM for a symmetric triangle
    let half = ((a / dtau).ceil() as isize).max(1);
    let mut k: Vec<f64> = (-half..=half)
        .map(|j| (1.0 - (j as f64 * dtau).abs() / a).max(0.0))
        .collect();
    normalize_sum(&mut k);
    k
}

fn normalize_sum(k: &mut [f64]) {
    let s: f64 = k.iter().sum();
    if s > 0.0 {
        for v in k.iter_mut() {
            *v /= s;
        }
    }
}

/// Discrete convolution with a centered symmetric `kernel` (odd length),
/// returning an output the same length as `input` (zero-padded edges).
fn convolve_same(input: &[f64], kernel: &[f64]) -> Vec<f64> {
    let n = input.len();
    let kh = (kernel.len() / 2) as isize;
    let mut out = vec![0.0f64; n];
    for (i, o) in out.iter_mut().enumerate() {
        let mut acc = 0.0;
        for (kk, &kv) in kernel.iter().enumerate() {
            let src = i as isize + (kk as isize - kh);
            if src >= 0 && (src as usize) < n {
                acc += input[src as usize] * kv;
            }
        }
        *o = acc;
    }
    out
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::resolution::{
        ResolutionFunction, apply_resolution, apply_resolution_with_plan, build_resolution_plan,
    };
    use std::sync::Arc;

    /// Trapezoidal integral of `I(τ)` over a fine grid out to many decay times.
    fn pulse_area(alpha: f64, beta: f64, r: f64) -> f64 {
        let tau_max = (18.0 / alpha).max(if r > 0.0 { 18.0 / beta } else { 0.0 });
        let n = 200_000;
        let dt = tau_max / n as f64;
        let mut area = 0.0;
        for i in 0..n {
            let t0 = i as f64 * dt;
            let t1 = (i + 1) as f64 * dt;
            area += 0.5 * (ic_pulse(alpha, beta, r, t0) + ic_pulse(alpha, beta, r, t1)) * dt;
        }
        area
    }

    fn pulse_mean(alpha: f64, beta: f64, r: f64) -> f64 {
        let tau_max = (24.0 / alpha).max(if r > 0.0 { 24.0 / beta } else { 0.0 });
        let n = 400_000;
        let dt = tau_max / n as f64;
        let mut m = 0.0;
        for i in 0..n {
            let t = (i as f64 + 0.5) * dt;
            m += t * ic_pulse(alpha, beta, r, t) * dt;
        }
        m
    }

    #[test]
    fn pulse_is_unit_area() {
        for &(a, b, r) in &[
            (0.5, 0.05, 0.0),
            (1.0, 0.1, 0.3),
            (2.0, 0.2, 0.6),
            (0.8, 0.5, 0.9),
        ] {
            let area = pulse_area(a, b, r);
            assert!(
                (area - 1.0).abs() < 1e-3,
                "area for (α={a},β={b},R={r}) = {area}, expected 1"
            );
        }
    }

    #[test]
    fn pulse_mean_matches_formula() {
        // ⟨τ⟩ = 3/α + R/β
        for &(a, b, r) in &[(1.0, 0.1, 0.0), (1.0, 0.1, 0.4), (2.0, 0.25, 0.7)] {
            let want = 3.0 / a + r / b;
            let got = pulse_mean(a, b, r);
            assert!(
                (got - want).abs() / want < 2e-3,
                "mean (α={a},β={b},R={r}) = {got}, expected {want}"
            );
        }
    }

    #[test]
    fn pulse_mode_is_two_over_alpha_for_pure_fast() {
        // r=0 ⇒ Gamma(3,α): mode at τ=2/α.
        let alpha = 1.3;
        let mode = 2.0 / alpha;
        let here = ic_pulse(alpha, 0.1, 0.0, mode);
        for d in [-0.3, -0.1, 0.1, 0.3] {
            assert!(ic_pulse(alpha, 0.1, 0.0, mode + d) <= here + 1e-12);
        }
    }

    #[test]
    fn pulse_alpha_equals_beta_is_finite_gamma4() {
        // α=β limit: slow term → α⁴τ³/6·e^{−ατ} (Gamma(4)). r=1 isolates it.
        let a = 1.0;
        let tau = 2.5;
        let got = ic_pulse(a, a, 1.0, tau);
        let want = a.powi(4) * tau.powi(3) / 6.0 * (-a * tau).exp();
        assert!(got.is_finite());
        assert!(
            (got - want).abs() < 1e-9,
            "α=β pulse {got} != Gamma(4) {want}"
        );
        // Near-degenerate (β just below α) must also be finite and close.
        let near = ic_pulse(a, a - 1e-9, 1.0, tau);
        assert!(near.is_finite() && (near - want).abs() < 1e-6);
    }

    #[test]
    fn pulse_is_finite_for_beta_much_greater_than_alpha() {
        // old e^{−βτ}·h(u) form factored an e^{|u|} that overflowed → 0·∞ = NaN).
        for &tau in &[0.0, 1.0, 50.0, 400.0, 1000.0] {
            let v = ic_pulse(0.05, 4.0, 0.5, tau);
            assert!(
                v.is_finite() && v >= 0.0,
                "ic_pulse(0.05,4,0.5,{tau}) = {v}"
            );
        }
        // The synthesized kernel must contain no non-finite entries.
        let p = IkedaCarpenterParams {
            alpha: EnergyLaw::Const(0.05),
            beta: EnergyLaw::Const(4.0),
            r: EnergyLaw::Const(0.5),
            burst_sigma_us: None,
            channel_fwhm_us: None,
            detector_tail: None,
        };
        let (offs, wts) = synth_kernel(&p, 600, 1.0).unwrap();
        assert!(offs.iter().chain(wts.iter()).all(|v| v.is_finite()));
    }

    #[test]
    fn pulse_is_nonnegative_and_zero_before_t0() {
        assert_eq!(ic_pulse(1.0, 0.1, 0.5, -0.5), 0.0);
        for i in 0..200 {
            let t = i as f64 * 0.1;
            assert!(ic_pulse(1.0, 0.1, 0.5, t) >= 0.0);
        }
    }

    #[test]
    fn energy_law_eval() {
        assert_eq!(EnergyLaw::Const(3.0).eval(50.0), 3.0);
        let s = EnergyLaw::SqrtE { a0: 0.2, a1: 0.1 };
        assert!((s.eval(100.0) - (0.2 * 10.0 + 0.1)).abs() < 1e-12);
        // InverseLambda: α grows with E (λ shrinks).
        let il = EnergyLaw::InverseLambda { a0: 0.1, a1: 0.5 };
        assert!(il.eval(200.0) > il.eval(5.0));
        // R = exp(−E_meV/κ) → ~0 at eV-scale energies, ~1 at sub-meV.
        let rr = EnergyLaw::ExpMilliEv { kappa: 25.0 };
        assert!(rr.eval(10.0) < 1e-6); // 10 eV
        assert!(rr.eval(0.001) > 0.9); // 1 meV
    }

    /// The synthesized table and the model that produced it place the pulse at
    /// the same instant.
    ///
    /// Both are evaluated on the detector clock for the same `timing_offset_us`,
    /// so their mean arrival and their shape must agree.
    #[test]
    fn synthesized_table_agrees_with_the_model_it_came_from() {
        let model = IkedaCarpenter::new(
            IkedaCarpenterParams {
                alpha: EnergyLaw::SqrtE { a0: 0.35, a1: 0.05 },
                beta: EnergyLaw::Const(0.25),
                r: EnergyLaw::Const(0.15),
                burst_sigma_us: None,
                channel_fwhm_us: None,
                detector_tail: None,
            },
            25.0,
            &SynthesisGrid {
                e_min_ev: 1.0,
                e_max_ev: 100.0,
                n_energies: 16,
                n_tau: 512,
            },
        )
        .expect("valid IC");
        let table = model.tabulated();

        // Between reference energies, so the width blend is exercised: it
        // scales offsets about 0, which is the emission instant.
        for energy in [3.7_f64, 17.3, 55.0] {
            let tof = TOF_FACTOR * 25.0 / energy.sqrt();
            let edges: Vec<f64> = (0..=3000).map(|i| tof - 5.0 + i as f64 * 0.01).collect();
            let mean = |p: &[f64]| -> f64 {
                let mass: f64 = p.iter().sum();
                p.iter()
                    .enumerate()
                    .map(|(i, w)| w * 0.5 * (edges[i] + edges[i + 1]))
                    .sum::<f64>()
                    / mass
            };
            let from_model = model
                .detector_bin_probabilities(energy, &edges, 0.0)
                .expect("model evaluates");
            let from_table = table
                .detector_bin_probabilities(energy, &edges, 0.0)
                .expect("table evaluates");
            let peak = from_model.iter().copied().fold(0.0_f64, f64::max);
            let shape = from_model
                .iter()
                .zip(&from_table)
                .map(|(a, b)| (a - b).abs())
                .fold(0.0_f64, f64::max);
            assert!(
                shape < 0.01 * peak,
                "E={energy}: table differs from the analytic pulse by {:.2}% of peak",
                100.0 * shape / peak
            );
            let gap = mean(&from_model) - mean(&from_table);
            assert!(
                gap.abs() < 0.01,
                "E={energy}: model puts the pulse at {:.4} µs, its own table at \
                 {:.4} µs ({gap:+.4} µs apart)",
                mean(&from_model),
                mean(&from_table)
            );
        }
    }

    #[test]
    fn kernel_tail_points_to_positive_offset() {
        // Asymmetry: longer/heavier tail toward +offset (later TOF, lower E).
        let p = IkedaCarpenterParams::constant(1.0, 0.1, 0.2);
        let (offsets, weights) = synth_kernel(&p, 600, 10.0).unwrap();
        let max_pos = offsets.iter().cloned().fold(f64::MIN, f64::max);
        let min_neg = offsets.iter().cloned().fold(f64::MAX, f64::min);
        assert!(
            max_pos > min_neg.abs(),
            "expected longer +offset tail: +{max_pos} vs −{}",
            min_neg.abs()
        );
        let pos_w: f64 = offsets
            .iter()
            .zip(&weights)
            .filter(|(o, _)| **o > 0.0)
            .map(|(_, w)| *w)
            .sum();
        let neg_w: f64 = offsets
            .iter()
            .zip(&weights)
            .filter(|(o, _)| **o < 0.0)
            .map(|(_, w)| *w)
            .sum();
        assert!(pos_w > neg_w, "expected more weight at +offset");
    }

    #[test]
    fn higher_energy_gives_narrower_kernel() {
        // α(E) ∝ √E ⇒ prompt width 1/α shrinks with E ⇒ smaller TOF support.
        let p = IkedaCarpenterParams {
            alpha: EnergyLaw::SqrtE { a0: 0.3, a1: 0.0 },
            beta: EnergyLaw::Const(0.1),
            r: EnergyLaw::Const(0.0),
            burst_sigma_us: None,
            channel_fwhm_us: None,
            detector_tail: None,
        };
        let support = |e: f64| {
            let (o, _) = synth_kernel(&p, 600, e).unwrap();
            o.iter().cloned().fold(f64::MIN, f64::max) - o.iter().cloned().fold(f64::MAX, f64::min)
        };
        assert!(support(100.0) < support(5.0));
    }

    #[test]
    fn synthesize_builds_valid_ascending_table() {
        let p = IkedaCarpenterParams::constant(1.0, 0.1, 0.2);
        let grid = SynthesisGrid {
            e_min_ev: 0.5e-3,
            e_max_ev: 1000.0,
            n_energies: 32,
            n_tau: 400,
        };
        let ic = IkedaCarpenter::new(p, 25.0, &grid).expect("synthesis");
        assert_eq!(ic.ref_energies().len(), 32);
        assert_eq!(ic.tabulated().ref_energies().len(), 32);
        for w in ic.ref_energies().windows(2) {
            assert!(w[1] > w[0]);
        }
    }

    #[test]
    fn rejects_bad_config() {
        let p = IkedaCarpenterParams::constant(1.0, 0.1, 0.2);
        let bad = SynthesisGrid {
            e_min_ev: 1.0,
            e_max_ev: 0.5,
            n_energies: 16,
            n_tau: 100,
        };
        assert!(IkedaCarpenter::new(p.clone(), 25.0, &bad).is_err());
        assert!(IkedaCarpenter::new(p, -1.0, &SynthesisGrid::new(1.0, 10.0)).is_err());
        // A parameter law that yields α(E) ≤ 0 is rejected, not silently clamped.
        let neg_alpha = IkedaCarpenterParams {
            alpha: EnergyLaw::Const(-1.0),
            ..IkedaCarpenterParams::constant(1.0, 0.1, 0.0)
        };
        assert!(IkedaCarpenter::new(neg_alpha, 25.0, &SynthesisGrid::new(1.0, 100.0)).is_err());
        // Negative / non-finite burst or channel widths are rejected up front
        // (not `.abs()`-masked or NaN-swallowed during synthesis).
        for bad_width in [-1.0, f64::NAN, f64::INFINITY] {
            let neg_burst = IkedaCarpenterParams {
                burst_sigma_us: Some(bad_width),
                ..IkedaCarpenterParams::constant(1.0, 0.1, 0.0)
            };
            assert!(
                IkedaCarpenter::new(neg_burst, 25.0, &SynthesisGrid::new(1.0, 100.0)).is_err(),
                "burst_sigma_us={bad_width} should be rejected"
            );
            let neg_chan = IkedaCarpenterParams {
                channel_fwhm_us: Some(bad_width),
                ..IkedaCarpenterParams::constant(1.0, 0.1, 0.0)
            };
            assert!(
                IkedaCarpenter::new(neg_chan, 25.0, &SynthesisGrid::new(1.0, 100.0)).is_err(),
                "channel_fwhm_us={bad_width} should be rejected"
            );
        }
    }

    #[test]
    fn broadens_constant_to_constant() {
        // An area-normalized kernel preserves a flat spectrum.
        let p = IkedaCarpenterParams::constant(1.0, 0.1, 0.3);
        let ic = IkedaCarpenter::new(p, 25.0, &SynthesisGrid::new(1.0, 200.0)).unwrap();
        let energies: Vec<f64> = (0..400).map(|i| 1.0 + i as f64 * 0.5).collect();
        let spectrum = vec![0.7f64; energies.len()];
        let res = ResolutionFunction::IkedaCarpenter(Arc::new(ic));
        let out = apply_resolution(&energies, &spectrum, &res).unwrap();
        // Interior points (away from grid edges where the kernel is clipped).
        for v in &out[40..energies.len() - 40] {
            assert!((v - 0.7).abs() < 1e-3, "flat broadening drifted: {v}");
        }
    }

    #[test]
    fn ic_centering_shifts_broadened_symmetric_dip_with_alpha() {
        // The IC kernel anchors its MODE at offset 0, but the right-skewed pulse's
        // intensity centroid lags the mode by ~1/α in TOF (module docstring).
        // Broadening a symmetric-in-energy resonance therefore moves the dip
        // minimum off the nominal energy by an α-dependent amount. This guards the
        // documented bias (the loop-closure calibration tests cannot see it) and
        // pins it against a silent "fix" (e.g. re-centering on the centroid).
        const E0: f64 = 20.0;
        // Absorption-weighted centroid of the broadened dip — a robust position
        // estimator (the bare minimum is fragile under the wide low-α kernel).
        fn dip_centroid_energy(alpha: f64) -> f64 {
            // Constant-α kernel so the shape (hence the mode→mean lag) is uniform.
            let p = IkedaCarpenterParams::constant(alpha, 0.1, 0.0);
            let ic = IkedaCarpenter::new(p, 25.0, &SynthesisGrid::new(1.0, 200.0)).unwrap();
            let res = ResolutionFunction::IkedaCarpenter(Arc::new(ic));
            // Fine uniform grid; symmetric (in energy) Gaussian dip centered at E0.
            let energies: Vec<f64> = (0..4000).map(|i| 10.0 + i as f64 * 0.005).collect();
            let spectrum: Vec<f64> = energies
                .iter()
                .map(|&e| 1.0 - 0.8 * (-((e - E0) / 0.1).powi(2)).exp())
                .collect();
            let out = apply_resolution(&energies, &spectrum, &res).unwrap();
            let (mut num, mut den) = (0.0, 0.0);
            for (&e, &t) in energies.iter().zip(&out) {
                if (e - E0).abs() <= 2.0 {
                    let a = (1.0 - t).max(0.0); // absorption weight
                    num += e * a;
                    den += a;
                }
            }
            num / den
        }
        let signed_small_alpha = dip_centroid_energy(0.8) - E0;
        let signed_large_alpha = dip_centroid_energy(2.0) - E0;
        // The shift is toward LOWER apparent energy: the delayed-emission
        // tail gathers theory from earlier TOF (higher E), so the theory
        // dip at E0 surfaces at measured energies below E0. A positive
        // shift here means the kernel was applied time-mirrored.
        assert!(
            signed_small_alpha < 0.0 && signed_large_alpha < 0.0,
            "centering shift must be toward lower energy: \
             α=0.8 {signed_small_alpha:+e}, α=2.0 {signed_large_alpha:+e}"
        );
        let shift_small_alpha = signed_small_alpha.abs();
        let shift_large_alpha = signed_large_alpha.abs();
        // The bias is real (resolvable at the 1e-3 eV level)…
        assert!(
            shift_small_alpha > 1e-3,
            "centering shift vanished: {shift_small_alpha}"
        );
        // …and shrinks with increasing α (the ~1/α scaling of the mode→mean lag).
        assert!(
            shift_small_alpha > 1.3 * shift_large_alpha,
            "shift should scale ~1/α: α=0.8 {shift_small_alpha} vs α=2.0 {shift_large_alpha}"
        );
    }

    #[test]
    fn plan_path_matches_direct_path() {
        let p = IkedaCarpenterParams::constant(1.2, 0.15, 0.25);
        let ic = IkedaCarpenter::new(p, 25.0, &SynthesisGrid::new(1.0, 200.0)).unwrap();
        let res = ResolutionFunction::IkedaCarpenter(Arc::new(ic));
        let energies: Vec<f64> = (0..300).map(|i| 1.0 + i as f64 * 0.6).collect();
        // A localized dip to exercise the convolution non-trivially.
        let spectrum: Vec<f64> = energies
            .iter()
            .map(|&e| if (e - 90.0).abs() < 5.0 { 0.2 } else { 1.0 })
            .collect();
        let direct = apply_resolution(&energies, &spectrum, &res).unwrap();
        let plan = build_resolution_plan(&energies, &res).unwrap();
        let planned =
            apply_resolution_with_plan(plan.as_ref(), &energies, &spectrum, &res).unwrap();
        assert_eq!(direct.len(), planned.len());
        for (a, b) in direct.iter().zip(&planned) {
            assert!((a - b).abs() < 1e-12, "plan vs direct mismatch: {a} vs {b}");
        }
    }

    #[test]
    fn tau_step_anchors_to_prompt_core_not_storage_tail() {
        // Regression for the τ-grid fix: with a slow storage tail active
        // (β ≪ α, R > 0) the τ-step must stay anchored to the prompt core —
        // the old `tau_max/(n_tau−1)` step let the tail dilate dtau and
        // degenerate the 0.35 µs channel triangle toward a delta.
        let n_tau = 600;
        let alpha = 2.0; // fast_reach = 18/α = 9 µs
        let prompt_step = (18.0 / alpha) / (n_tau as f64 - 1.0);

        let p = IkedaCarpenterParams {
            channel_fwhm_us: Some(0.35),
            ..IkedaCarpenterParams::constant(alpha, 0.25, 0.5)
        };
        let (offs, wts) = synth_kernel(&p, n_tau, 10.0).unwrap();
        let dtau = offs[1] - offs[0];
        assert!(
            (dtau - prompt_step).abs() < 1e-12,
            "uncapped step {dtau} != prompt-core step {prompt_step}"
        );
        assert!(
            dtau <= prompt_step + 1e-12,
            "prompt-core spacing {dtau} > fast_reach/(n_tau−1) = {prompt_step}"
        );
        // ≥ 3 nonzero triangle samples per side at this step.
        let tri = triangle_kernel(dtau, 0.35);
        let nonzero_per_side = tri.iter().take(tri.len() / 2).filter(|&&v| v > 0.0).count();
        assert!(
            nonzero_per_side >= 3,
            "triangle degenerated: {nonzero_per_side} nonzero samples per side"
        );
        // The tail is still reached: for β = 0.25 the slow tail crosses the
        // TRIM_REL = 1e-7 trim level near τ ≈ 63 µs (e^{−βτ}·slow-amplitude ÷
        // peak = 1e-7 at τ ≈ 62.6 µs), so the trimmed kernel must still span
        // well past 40 µs — a step anchored to τ_max instead of the prompt
        // core would pass a narrow-span check, but a grid that stops short of
        // the storage tail would not survive this one.
        let span = offs.last().unwrap() - offs.first().unwrap();
        let peak = wts.iter().cloned().fold(f64::MIN, f64::max);
        assert!((peak - 1.0).abs() < 1e-12, "weights not peak-normalized");
        assert!(span > 3.0 / alpha, "kernel span {span} lost the pulse body");
        assert!(
            span > 40.0,
            "kernel span {span} µs stops short of the β = 0.25 storage tail \
             (trim horizon ≈ 63 µs)"
        );

        // (b) Extreme admitted tail (β = 0.02 ⇒ slow_reach = 800 µs): the
        // MAX_TAU_SAMPLES cap bites; the step widens to tau_max/(cap−1) but
        // must still resolve the 0.35 µs triangle with ≥ 3 samples per side
        // (the old tail-anchored step was 800/599 ≈ 1.34 µs — a delta).
        let p_ext = IkedaCarpenterParams {
            channel_fwhm_us: Some(0.35),
            ..IkedaCarpenterParams::constant(alpha, 0.02, 0.5)
        };
        let (offs_ext, _) = synth_kernel(&p_ext, n_tau, 10.0).unwrap();
        let dtau_ext = offs_ext[1] - offs_ext[0];
        let capped_step = 800.0 / (MAX_TAU_SAMPLES as f64 - 1.0);
        assert!(
            (dtau_ext - capped_step).abs() < 1e-9,
            "capped step {dtau_ext} != tau_max/(MAX_TAU_SAMPLES−1) = {capped_step}"
        );
        assert!(
            dtau_ext <= 0.35 / 3.0,
            "capped step {dtau_ext} cannot resolve the 0.35 µs triangle"
        );
        let tri_ext = triangle_kernel(dtau_ext, 0.35);
        let nonzero_ext = tri_ext
            .iter()
            .take(tri_ext.len() / 2)
            .filter(|&&v| v > 0.0)
            .count();
        assert!(
            nonzero_ext >= 3,
            "triangle degenerated under cap: {nonzero_ext} nonzero samples per side"
        );
    }

    #[test]
    fn burst_and_channel_broaden_further() {
        // Folding in burst + channel widens the kernel TOF support.
        let base = IkedaCarpenterParams::constant(1.0, 0.1, 0.0);
        let (o0, _) = synth_kernel(&base, 600, 10.0).unwrap();
        let support0 = o0.iter().cloned().fold(f64::MIN, f64::max)
            - o0.iter().cloned().fold(f64::MAX, f64::min);
        let folded = IkedaCarpenterParams {
            burst_sigma_us: Some(0.3),
            channel_fwhm_us: Some(0.35),
            ..IkedaCarpenterParams::constant(1.0, 0.1, 0.0)
        };
        let (o1, _) = synth_kernel(&folded, 600, 10.0).unwrap();
        let support1 = o1.iter().cloned().fold(f64::MIN, f64::max)
            - o1.iter().cloned().fold(f64::MAX, f64::min);
        assert!(support1 > support0, "folded {support1} !> bare {support0}");
    }

    /// Sum-weighted variance of a sampled `(offsets, weights)` kernel on a
    /// uniform τ-grid. Normalization-independent (divides by Σw).
    fn kernel_variance(offsets: &[f64], weights: &[f64]) -> f64 {
        let w_sum: f64 = weights.iter().sum();
        let mean: f64 = offsets.iter().zip(weights).map(|(o, w)| o * w).sum::<f64>() / w_sum;
        offsets
            .iter()
            .zip(weights)
            .map(|(o, w)| (o - mean).powi(2) * w)
            .sum::<f64>()
            / w_sum
    }

    #[test]
    fn unresolvable_tau_grid_is_rejected_loudly() {
        // Review F1, probe 1: β = 0.005 ⇒ slow reach 16/β = 3200 µs ⇒
        // capped step 3200/8191 ≈ 0.39 µs > the 0.35 µs triangle FWHM — the
        // sampled triangle would be the exact delta [0, 1, 0] and the
        // requested fold would silently vanish from the kernel. Must be a
        // loud construction error, not silent physics degradation.
        let p1 = IkedaCarpenterParams {
            channel_fwhm_us: Some(0.35),
            ..IkedaCarpenterParams::constant(2.0, 0.005, 0.5)
        };
        let err = IkedaCarpenter::new(p1.clone(), 25.0, &SynthesisGrid::new(1.0, 100.0))
            .expect_err("a delta-degenerate fold must be rejected");
        let msg = format!("{err:?}");
        assert!(
            msg.contains("cannot resolve") && msg.contains("Increase"),
            "error must diagnose the cap and name remedies: {msg}"
        );
        // The per-energy synthesis (kernel_at path) refuses identically.
        assert!(synth_kernel(&p1, 600, 10.0).is_err());

        // Probe 2: α = 250 ⇒ the whole prompt pulse spans 18/α ≈ 0.07 µs,
        // less than ONE capped step (800/8191 ≈ 0.098 µs): sampling would
        // step over the prompt term entirely (0.9 of the pulse weight at
        // R = 0.1). Must also be rejected.
        let p2 = IkedaCarpenterParams::constant(250.0, 0.02, 0.1);
        assert!(
            IkedaCarpenter::new(p2, 25.0, &SynthesisGrid::new(1.0, 100.0)).is_err(),
            "a capped step wider than the prompt core must be rejected"
        );
    }

    #[test]
    fn requested_fold_is_never_a_silent_no_op() {
        // Adjacent to probe 1 but resolvable (β = 0.05 ⇒ capped step
        // 320/8191 ≈ 0.039 µs ≤ FWHM/3): when synthesis is Ok and a fold is
        // requested, the folded kernel must genuinely differ from the
        // unfolded one — by approximately the analytic fold variance
        // FWHM²/6 — never by ~0 (the silent delta no-op this guards against).
        let base = IkedaCarpenterParams::constant(2.0, 0.05, 0.5);
        let folded = IkedaCarpenterParams {
            channel_fwhm_us: Some(0.35),
            ..base.clone()
        };
        let (o0, w0) = synth_kernel(&base, 600, 10.0).unwrap();
        let (o1, w1) = synth_kernel(&folded, 600, 10.0).unwrap();
        let dv = kernel_variance(&o1, &w1) - kernel_variance(&o0, &w0);
        let expected = 0.35f64.powi(2) / 6.0;
        assert!(
            dv > 0.5 * expected,
            "fold variance increment {dv} µs² vs analytic {expected} µs²: \
             the requested fold (nearly) vanished"
        );
    }

    #[test]
    fn boundary_step_fwhm_over_3_is_admitted_and_fold_survives() {
        // Pin the exactly-admitted boundary `dtau = FWHM/3` (rejection is
        // STRICT: `capped_step > floor`). Per the TRI_MIN_SAMPLES_PER_SIDE
        // doc, each triangle side then samples
        // at {FWHM/3, 2FWHM/3, FWHM} — weights {2/3, 1/3, 0}: exactly 2
        // strictly interior nonzero samples, endpoint ON the triangle zero —
        // and the fold stays effective (non-delta, second moment ≈ 4FWHM²/27).
        //
        // Route to the boundary: n_tau = MIN_N_TAU = 8 with α = 1 gives a
        // prompt design step 18/7 ≈ 2.57 µs, far coarser than FWHM/3 for a
        // 1 µs triangle, so the fold refinement pins dtau to exactly
        // `fwhm / TRI_MIN_SAMPLES_PER_SIDE`. R = 0 keeps the cap inert
        // (capped step 18/8191 ≪ floor).
        let fwhm = 1.0;
        let p = IkedaCarpenterParams {
            channel_fwhm_us: Some(fwhm),
            ..IkedaCarpenterParams::constant(1.0, 0.1, 0.0)
        };
        let (offs, _) = synth_kernel(&p, MIN_N_TAU, 10.0)
            .expect("the dtau = FWHM/3 boundary must be admitted, not rejected");
        let dtau = offs[1] - offs[0];
        let boundary = fwhm / TRI_MIN_SAMPLES_PER_SIDE;
        assert!(
            (dtau - boundary).abs() < 1e-12,
            "step {dtau} µs is not the FWHM/3 boundary {boundary} µs"
        );

        // The sampled triangle at the boundary: ≥ 7 samples, exactly 2
        // strictly interior nonzero per side, center far below the delta's 1.
        let tri = triangle_kernel(dtau, fwhm);
        assert!(tri.len() >= 7, "boundary triangle too short: {}", tri.len());
        let per_side = tri
            .iter()
            .take(tri.len() / 2)
            .filter(|&&v| v > 1e-9)
            .count();
        assert_eq!(
            per_side, 2,
            "boundary triangle must keep 2 strictly interior nonzero samples \
             per side: {tri:?}"
        );
        let center = tri[tri.len() / 2];
        assert!(
            center < 0.5,
            "center weight {center} — boundary triangle degenerated toward a delta"
        );

        // Fold effectiveness: the discrete triangle's variance is 4·FWHM²/27
        // (samples {0, ±FWHM/3, ±2FWHM/3} with weights ∝ {1, 2/3, 1/3}),
        // ≈ 89 % of the analytic FWHM²/6 — the fold physics survives the
        // 2-interior-sample discretization.
        let tri_offs: Vec<f64> = (0..tri.len())
            .map(|i| (i as f64 - (tri.len() / 2) as f64) * dtau)
            .collect();
        let v = kernel_variance(&tri_offs, &tri);
        let want = 4.0 * fwhm * fwhm / 27.0;
        assert!(
            ((v - want) / want).abs() < 0.02,
            "boundary triangle variance {v} µs² vs discrete analytic {want} µs²"
        );
        assert!(
            v > 0.5 * (fwhm * fwhm / 6.0),
            "boundary fold variance {v} µs² collapsed below half the \
             continuous FWHM²/6 — fold (nearly) vanished"
        );
    }

    #[test]
    fn fold_variance_matches_analytic_oracle() {
        // Independent analytic oracle for the fold convention: a
        // convolution adds second central moments, so the sampled kernel's
        // variance must grow by EXACTLY the fold kernel's analytic variance —
        // FWHM²/6 for the symmetric channel triangle (half-base = FWHM),
        // σ² for the Gaussian burst. Unlike the closed-loop calibration
        // tests, the expectation here comes from analysis, not from the
        // shared synthesis code. Pure-prompt pulse (R = 0) keeps trim /
        // edge-truncation error far below the increments.
        let base = IkedaCarpenterParams::constant(2.0, 0.1, 0.0);
        let (o0, w0) = synth_kernel(&base, 600, 10.0).unwrap();
        let v0 = kernel_variance(&o0, &w0);
        // Sanity: the bare Gamma(3, α) variance is 3/α².
        let v_gamma = 3.0 / (2.0f64 * 2.0);
        assert!(
            (v0 - v_gamma).abs() < 0.01 * v_gamma,
            "bare pulse variance {v0}, Gamma(3) analytic {v_gamma}"
        );

        let fwhm = 0.35;
        let tri = IkedaCarpenterParams {
            channel_fwhm_us: Some(fwhm),
            ..base.clone()
        };
        let (o1, w1) = synth_kernel(&tri, 600, 10.0).unwrap();
        let dv_tri = kernel_variance(&o1, &w1) - v0;
        let want_tri = fwhm * fwhm / 6.0;
        assert!(
            ((dv_tri - want_tri) / want_tri).abs() < 0.02,
            "triangle fold added {dv_tri} µs², analytic FWHM²/6 = {want_tri} µs²"
        );

        let sigma = 0.3;
        let gau = IkedaCarpenterParams {
            burst_sigma_us: Some(sigma),
            ..base.clone()
        };
        let (o2, w2) = synth_kernel(&gau, 600, 10.0).unwrap();
        let dv_gau = kernel_variance(&o2, &w2) - v0;
        let want_gau = sigma * sigma;
        assert!(
            ((dv_gau - want_gau) / want_gau).abs() < 0.02,
            "Gaussian burst added {dv_gau} µs², analytic σ² = {want_gau} µs²"
        );
    }
}