nereids_physics/

resolution.rs

//! Resolution broadening via convolution with instrument resolution function.
//!
//! Convolves theoretical cross-sections (or transmission) with the instrument
//! resolution function to account for finite energy resolution. The resolution
//! function is modeled as a Gaussian with energy-dependent width, optionally
//! combined with an exponential tail, derived from time-of-flight instrument
//! parameters.
//!
//! ## SAMMY Reference
//! - `rsl/mrsl1.f90` — Main RSL resolution broadening routines (Resbrd)
//! - `rsl/mrsl4.f90` — Resolution width calculation (Wdsint, Rolowg)
//! - `rsl/mrsl5.f90` — Exponential tail peak shift (Shftge)
//! - `fnc/exerfc.f90` — Scaled complementary error function
//! - `convolution/DopplerAndResolutionBroadener.cpp` — Xcoef quadrature weights
//! - Manual Section III.C (Resolution Broadening); quadrature Eq. IV B 3.8
//!   (R3-revision numbering — see `compute_xcoef_weights` and the
//!   Gaussian+exponential path in `resolution_broaden_presorted`)
//!
//! ## Width convention — read before supplying a number
//!
//! Every width here is a **W-parameter**, the width appearing in
//! `exp(-x²/W²)`, not a standard deviation. The two differ by √2:
//!
//!   σ = W/√2        FWHM = 2·√(ln 2)·W = 1.6651·W
//!
//! This is SAMMY's convention, and it is the same one the Doppler width Δ_D
//! uses, so the two broadening kernels compose without a conversion. Supplying
//! a 1σ value where a W is expected yields a kernel √2 too narrow — 29 % — and
//! a resolution that is too narrow is absorbed into a fitted temperature that
//! is too high.
//!
//! [`ResolutionParams::from_sigma`] and [`ResolutionParams::from_fwhm`] convert
//! for you. Use them rather than scaling by hand.
//!
//! SAMMY's own inputs are in yet other measures — `Deltag` is a FWHM and
//! `Deltal` is the full width of a rectangular path spread — so reading a
//! SAMMY `.inp` goes through `nereids_endf::sammy::sammy_to_nereids_resolution`
//! (`Deltag/(2√ln2)`, `Deltal/√6`), never straight into these fields.
//!
//! ## Physics
//!
//! For a time-of-flight instrument, the energy resolution is:
//!
//!   (ΔE/E)² = (2·Δt/t)² + (2·ΔL/L)²
//!
//! where t = L/v is the neutron time-of-flight, Δt is the total timing width,
//! and ΔL is the flight path width — both W-parameters, as above, so ΔE is one
//! too. The factor 2 is the kinematic derivative dE/E = 2·dt/t; it is not a
//! width-measure conversion. Since t ∝ 1/√E, the timing contribution gives
//! ΔE ∝ E^(3/2) while the path contribution gives ΔE ∝ E.
//!
//! The broadened cross-section is:
//!
//!   σ_res(E) = ∫ R(E, E') · σ(E') dE'
//!
//! When Deltae = 0, R is a pure Gaussian (Iesopr=1):
//!   R(E, E') = exp(-(E-E')²/Wg²) / (Wg·√π)
//!
//! When Deltae > 0, R is the convolution of a Gaussian with an exponential
//! tail (Iesopr=3):
//!   R(E, E') ∝ exp(2·C·A + C²) · erfc(C + A)
//!
//! where C = Wg/(2·We), A = (E - E')/Wg, Wg = Gaussian width, We = exponential
//! width. This is the analytical result for convolving exp(-x²/Wg²) with
//! exp(-x/We)·H(x).

use nereids_core::constants::{DIVISION_FLOOR, NEAR_ZERO_FLOOR};
use std::f64::consts::SQRT_2;
use std::fmt;
use std::sync::Arc;

/// FWHM of `exp(-x²/W²)` divided by W: `2·√(ln 2)` = 1.6651.
///
/// The conversion between this module's W-parameters and a full width at half
/// maximum. Paired with σ = W/√2 it fixes all three width measures against each
/// other; see the module's width-convention section.
const FWHM_PER_W: f64 = 1.665_109_222_315_395_4;

/// TOF conversion factor: `t (μs) = TOF_FACTOR × L (m) / √(E in eV)`.
///
/// Derived from t = L / √(2E/m_n), converting to microseconds:
///   TOF_FACTOR = 1e6 / √(2 × EV_TO_JOULES / NEUTRON_MASS_KG)
///
/// Uses CODATA 2018 values (both exact in the 2019 SI).
///
/// `pub` so the analytical [`crate::ikeda_carpenter`] model and the
/// `nereids-fitting` resolution calibrator both use the *identical* TOF↔energy
/// constant (the calibrator's position nuisance shifts the grid in TOF) — any
/// drift here would make the IC-vs-tabulated cross-validation unfair.
pub const TOF_FACTOR: f64 = 72.298_254_398_292_8;

/// Errors from resolution broadening operations.
#[derive(Debug, PartialEq)]
pub enum ResolutionError {
    /// The energy grid is not sorted in ascending order.
    UnsortedEnergies,
    /// The energy grid and data arrays have mismatched lengths.
    LengthMismatch { energies: usize, data: usize },
    /// A [`ResolutionPlan`] was passed together with an `energies`
    /// slice that does not match the grid the plan was built for.
    ///
    /// Cheapest-available check hierarchy: length mismatch is caught
    /// first via [`Self::LengthMismatch`] (`plan.len() ==
    /// energies.len()` is necessary but not sufficient); a content
    /// mismatch fires `PlanGridMismatch` with the index of the first
    /// differing element so callers can diagnose silent-staleness
    /// bugs at the cache layer.
    PlanGridMismatch { first_diff_index: usize },
    /// A [`ResolutionMatrix`] was passed together with an `energies`
    /// slice that does not match the grid the matrix was compiled for.
    /// Same semantics as [`Self::PlanGridMismatch`] but for the CSR
    /// path (see [`apply_r`]).
    MatrixGridMismatch { first_diff_index: usize },
    /// [`TabulatedResolution::width_corrected`] was called with invalid
    /// parameters: `s0` must be finite and `> 0`, `e_ref` finite and `> 0`,
    /// and `p` finite. A non-positive `s0` would reverse/collapse the
    /// (ascending) offset ordering the broadening loop assumes.
    InvalidWidthCorrection { s0: f64, p: f64, e_ref: f64 },
}

impl fmt::Display for ResolutionError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::UnsortedEnergies => write!(
                f,
                "energy grid must be sorted in non-descending order for binary search"
            ),
            Self::LengthMismatch { energies, data } => write!(
                f,
                "energy grid length ({}) must match data length ({})",
                energies, data
            ),
            Self::PlanGridMismatch { first_diff_index } => write!(
                f,
                "resolution plan was built for a different energy grid than was \
                 passed to apply_resolution_with_plan (first differing index: {})",
                first_diff_index,
            ),
            Self::MatrixGridMismatch { first_diff_index } => write!(
                f,
                "resolution matrix was compiled for a different energy grid than was \
                 passed to apply_resolution_with_matrix (first differing index: {})",
                first_diff_index,
            ),
            Self::InvalidWidthCorrection { s0, p, e_ref } => write!(
                f,
                "width_corrected requires finite s0 > 0, finite e_ref > 0, and finite p; \
                 got s0={s0}, p={p}, e_ref={e_ref}"
            ),
        }
    }
}

impl std::error::Error for ResolutionError {}

/// Errors from `ResolutionParams` construction.
#[derive(Debug, PartialEq)]
pub enum ResolutionParamsError {
    /// Flight path must be positive and finite.
    InvalidFlightPath(f64),
    /// Timing uncertainty must be non-negative and finite.
    InvalidDeltaT(f64),
    /// Path length uncertainty must be non-negative and finite.
    InvalidDeltaL(f64),
    /// Exponential tail parameter must be non-negative and finite.
    InvalidDeltaE(f64),
}

impl fmt::Display for ResolutionParamsError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidFlightPath(v) => {
                write!(f, "flight_path_m must be positive and finite, got {v}")
            }
            Self::InvalidDeltaT(v) => {
                write!(f, "delta_t_us must be non-negative and finite, got {v}")
            }
            Self::InvalidDeltaL(v) => {
                write!(f, "delta_l_m must be non-negative and finite, got {v}")
            }
            Self::InvalidDeltaE(v) => {
                write!(f, "delta_e_us must be non-negative and finite, got {v}")
            }
        }
    }
}

impl std::error::Error for ResolutionParamsError {}

/// Resolution function parameters for time-of-flight instruments.
#[derive(Debug, Clone, Copy)]
pub struct ResolutionParams {
    /// Flight path length in meters (source to detector).
    flight_path_m: f64,
    /// Total timing width in microseconds, as a W-parameter (σ = W/√2).
    /// Combines moderator pulse width, detector timing, and electronics.
    ///
    /// Not a standard deviation — see the module's width-convention section.
    delta_t_us: f64,
    /// Flight path width in meters, as a W-parameter (σ = W/√2).
    delta_l_m: f64,
    /// Exponential tail parameter (SAMMY Deltae, raw SAMMY units).
    ///
    /// When zero, pure Gaussian broadening is used (SAMMY Iesopr=1).
    /// When positive, the kernel is the convolution of a Gaussian with an
    /// exponential tail (SAMMY Iesopr=3).
    ///
    /// SAMMY Ref: `RslResolutionFunction_M.f90` getCo2, `rsl/mrsl4.f90` Wdsint.
    delta_e_us: f64,
}

impl ResolutionParams {
    /// Create validated resolution parameters.
    ///
    /// `delta_t_us` and `delta_l_m` are W-parameters (σ = W/√2), not standard
    /// deviations. If your instrument numbers are 1σ or FWHM, use
    /// [`from_sigma`](Self::from_sigma) or [`from_fwhm`](Self::from_fwhm)
    /// instead of converting by hand.
    ///
    /// # Arguments
    /// * `flight_path_m` — Flight path length in meters (must be > 0).
    /// * `delta_t_us` — Timing width in microseconds, W-parameter (must be >= 0).
    /// * `delta_l_m` — Flight path width in meters, W-parameter (must be >= 0).
    /// * `delta_e_us` — Exponential tail parameter in SAMMY Deltae units
    ///   (must be >= 0). When 0, pure Gaussian broadening is used.
    ///
    /// # Errors
    /// Returns `ResolutionParamsError::InvalidFlightPath` if `flight_path_m <= 0.0`
    /// or is not finite.
    /// Returns `ResolutionParamsError::InvalidDeltaT` if `delta_t_us < 0.0` or is
    /// not finite.
    /// Returns `ResolutionParamsError::InvalidDeltaL` if `delta_l_m < 0.0` or is
    /// not finite.
    /// Returns `ResolutionParamsError::InvalidDeltaE` if `delta_e_us < 0.0` or is
    /// not finite.
    pub fn new(
        flight_path_m: f64,
        delta_t_us: f64,
        delta_l_m: f64,
        delta_e_us: f64,
    ) -> Result<Self, ResolutionParamsError> {
        if !flight_path_m.is_finite() || flight_path_m <= 0.0 {
            return Err(ResolutionParamsError::InvalidFlightPath(flight_path_m));
        }
        if !delta_t_us.is_finite() || delta_t_us < 0.0 {
            return Err(ResolutionParamsError::InvalidDeltaT(delta_t_us));
        }
        if !delta_l_m.is_finite() || delta_l_m < 0.0 {
            return Err(ResolutionParamsError::InvalidDeltaL(delta_l_m));
        }
        if !delta_e_us.is_finite() || delta_e_us < 0.0 {
            return Err(ResolutionParamsError::InvalidDeltaE(delta_e_us));
        }
        Ok(Self {
            flight_path_m,
            delta_t_us,
            delta_l_m,
            delta_e_us,
        })
    }

    /// Create resolution parameters from standard deviations.
    ///
    /// Instrument metrology is usually quoted as a 1σ jitter, while this type
    /// stores W-parameters. This constructor applies W = σ·√2 to the timing and
    /// flight-path terms, so a measured 1σ can be passed in directly.
    ///
    /// `delta_e_us` is not a Gaussian width — it is SAMMY's exponential tail
    /// parameter Deltae — so it is passed through unchanged.
    ///
    /// # Errors
    /// Same as [`new`](Self::new), on the converted values.
    pub fn from_sigma(
        flight_path_m: f64,
        sigma_t_us: f64,
        sigma_l_m: f64,
        delta_e_us: f64,
    ) -> Result<Self, ResolutionParamsError> {
        Self::new(
            flight_path_m,
            sigma_t_us * SQRT_2,
            sigma_l_m * SQRT_2,
            delta_e_us,
        )
    }

    /// Create resolution parameters from full widths at half maximum.
    ///
    /// Applies W = FWHM/(2·√(ln 2)) to the timing and flight-path terms. This
    /// is the measure SAMMY's `Deltag` uses, so a value taken from a SAMMY
    /// `.inp` timing field can be passed in directly.
    ///
    /// `delta_e_us` is SAMMY's exponential tail parameter, not a Gaussian
    /// width, so it is passed through unchanged.
    ///
    /// # Errors
    /// Same as [`new`](Self::new), on the converted values.
    pub fn from_fwhm(
        flight_path_m: f64,
        fwhm_t_us: f64,
        fwhm_l_m: f64,
        delta_e_us: f64,
    ) -> Result<Self, ResolutionParamsError> {
        Self::new(
            flight_path_m,
            fwhm_t_us / FWHM_PER_W,
            fwhm_l_m / FWHM_PER_W,
            delta_e_us,
        )
    }

    /// Returns the flight path length in meters.
    #[must_use]
    pub fn flight_path_m(&self) -> f64 {
        self.flight_path_m
    }

    /// Total timing width in microseconds, as a W-parameter (σ = W/√2).
    ///
    /// The factor of 2 in [`gaussian_width()`](Self::gaussian_width) comes from
    /// the energy-TOF derivative dE/E = 2·dt/t, not from a width-measure
    /// conversion — the convention is unchanged between time and energy.
    #[must_use]
    pub fn delta_t_us(&self) -> f64 {
        self.delta_t_us
    }

    /// Returns the flight path width in meters, as a W-parameter (σ = W/√2).
    #[must_use]
    pub fn delta_l_m(&self) -> f64 {
        self.delta_l_m
    }

    /// Returns the exponential tail parameter (SAMMY Deltae units).
    #[must_use]
    pub fn delta_e_us(&self) -> f64 {
        self.delta_e_us
    }

    /// Whether the exponential tail is active (Deltae > 0, SAMMY Iesopr=3).
    #[must_use]
    pub fn has_exponential_tail(&self) -> bool {
        self.delta_e_us > NEAR_ZERO_FLOOR
    }

    /// Exponential tail width Widexp(E) in eV.
    ///
    /// SAMMY Ref: `rsl/mrsl4.f90` Wdsint lines 55-56 (Kedxfw=false path):
    ///   `Widexp = E * Co2 * sqrt(E)` where `Co2 = 2·Deltae / (Sm2·Dist)`.
    ///
    /// Combined: `Widexp = 2·Deltae·E^(3/2) / (TOF_FACTOR·L)`.
    #[must_use]
    pub fn exp_width(&self, energy_ev: f64) -> f64 {
        if energy_ev <= 0.0 || self.delta_e_us <= 0.0 {
            return 0.0;
        }
        2.0 * self.delta_e_us * energy_ev.powf(1.5) / (TOF_FACTOR * self.flight_path_m)
    }

    /// Gaussian resolution width W_g(E) in eV — the W of `exp(-x²/W_g²)`.
    ///
    /// Combines timing and flight-path contributions in quadrature:
    ///   W_g² = (2·Δt/t × E)² + (2·ΔL/L × E)²
    ///
    /// where t = TOF_FACTOR × L / √E is the time-of-flight in μs. Δt and ΔL are
    /// W-parameters too, so the convention carries through unchanged; the
    /// standard deviation is `W_g/√2` and the FWHM is [`fwhm`](Self::fwhm).
    #[must_use]
    pub fn gaussian_width(&self, energy_ev: f64) -> f64 {
        if energy_ev <= 0.0 || self.flight_path_m <= 0.0 {
            return 0.0;
        }

        // Timing contribution: W_t = 2 × Δt × E^(3/2) / (TOF_FACTOR × L)
        let timing =
            2.0 * self.delta_t_us * energy_ev.powf(1.5) / (TOF_FACTOR * self.flight_path_m);

        // Path length contribution: W_L = 2 × ΔL × E / L
        let path = 2.0 * self.delta_l_m * energy_ev / self.flight_path_m;

        (timing * timing + path * path).sqrt()
    }

    /// FWHM of the resolution function at energy E, in eV.
    ///
    /// `FWHM = 2·√(ln 2)·W`, the conversion for `exp(-x²/W²)`. For the standard
    /// deviation instead, use `gaussian_width(E)/√2`.
    #[must_use]
    pub fn fwhm(&self, energy_ev: f64) -> f64 {
        FWHM_PER_W * self.gaussian_width(energy_ev)
    }
}

/// Apply Gaussian resolution broadening to cross-section data.
///
/// Convolves the input cross-sections with a Gaussian kernel whose width
/// varies with energy according to the instrument resolution function.
///
/// # Arguments
/// * `energies` — Energy grid in eV (must be sorted ascending).
/// * `cross_sections` — Cross-sections in barns at each energy point.
/// * `params` — Resolution function parameters.
///
/// # Returns
/// Resolution-broadened cross-sections on the same energy grid.
///
/// # Errors
/// Returns [`ResolutionError::LengthMismatch`] if the arrays differ in length,
/// or [`ResolutionError::UnsortedEnergies`] if the energy grid is not sorted
/// in non-descending order.
pub fn resolution_broaden(
    energies: &[f64],
    cross_sections: &[f64],
    params: &ResolutionParams,
) -> Result<Vec<f64>, ResolutionError> {
    validate_inputs(energies, cross_sections)?;
    Ok(resolution_broaden_presorted(
        energies,
        cross_sections,
        params,
    ))
}

/// Check that the energy grid is sorted and that its length matches the data.
fn validate_inputs(energies: &[f64], data: &[f64]) -> Result<(), ResolutionError> {
    if energies.len() != data.len() {
        return Err(ResolutionError::LengthMismatch {
            energies: energies.len(),
            data: data.len(),
        });
    }
    if !energies.windows(2).all(|w| w[0] <= w[1]) {
        return Err(ResolutionError::UnsortedEnergies);
    }
    Ok(())
}

// ─── Xcoef quadrature weights ──────────────────────────────────────────────────

/// Compute SAMMY's 4-point quadrature weights for a non-uniform energy grid.
///
/// Replaces the simple trapezoidal rule `de = (E[j+1] - E[j-1]) / 2` with
/// SAMMY's higher-order scheme from Eq. IV B 3.8 (page 80 of SAMMY manual R3).
///
/// SAMMY Ref: `convolution/DopplerAndResolutionBroadener.cpp`, `setXcoefWeights()`.
///
/// The weights include a correction term x2(k) that accounts for non-uniform
/// grid spacing, providing 4th-order accuracy on smooth grids.
///
/// Note: the returned weights are 12x the quantity in Eq. IV B 3.8. This
/// constant factor cancels during normalization (sum/norm), so the broadened
/// result is independent of the scaling.
fn compute_xcoef_weights(energies: &[f64]) -> Vec<f64> {
    let n = energies.len();
    if n == 0 {
        return vec![];
    }
    if n == 1 {
        return vec![1.0];
    }

    // SAMMY's 4-point quadrature weights (Eq. IV B 3.8, SAMMY Manual R3 p80).
    //
    // Uses a sliding window of 5 consecutive energies E[0..4] to compute
    // coefficients A[0..5] at each grid point k:
    //
    //   A[0] = v1  (k >= 2)
    //   A[1] = 5·v2  (k >= 1)
    //   A[2] = 5·v3  (k < n-1)
    //   A[3] = v4  (k < n-2)
    //   A[4] = (v3² - v1²)/v2   curvature correction  (k >= 2)
    //   A[5] = -(v4² - v2²)/v3  curvature correction  (k >= 1)
    //
    // where v1..v4 are consecutive grid spacings around point k.
    //
    // The result is 12× Eq. IV B 3.8; this constant factor cancels during
    // normalization (sum/norm) in the broadening loop.
    //
    // SAMMY Ref: `convolution/DopplerAndResolutionBroadener.cpp` lines 365-457
    let mut weights = vec![0.0f64; n];

    // Sliding window: e[j] holds energies relative to current k.
    // At loop start for k: e[0]=E[k-2], e[1]=E[k-1], e[2]=E[k],
    //                       e[3]=E[k+1], e[4]=E[k+2]
    // Out-of-bounds positions are 0.0 (matching SAMMY's convention).
    let mut e = [0.0f64; 5];
    e[3] = energies[0];
    if n > 1 {
        e[4] = energies[1];
    }

    for k in 0..n {
        // Shift window left.
        e[0] = e[1];
        e[1] = e[2];
        e[2] = e[3];
        e[3] = e[4];
        e[4] = if k + 2 < n { energies[k + 2] } else { 0.0 };

        let v1 = e[1] - e[0];
        let v2 = e[2] - e[1];
        let v3 = e[3] - e[2];
        let v4 = e[4] - e[3];

        let mut a = [0.0f64; 6];

        if k >= 2 {
            a[0] = v1;
            // Curvature correction: x2(k-2) = (v3² - v1²) / v2
            if v2.abs() > NEAR_ZERO_FLOOR {
                a[4] = (v3 * v3 - v1 * v1) / v2;
            }
        }
        if k >= 1 {
            a[1] = 5.0 * v2;
            // Curvature correction: -x2(k-1) = -(v4² - v2²) / v3
            if v3.abs() > NEAR_ZERO_FLOOR {
                a[5] = -(v4 * v4 - v2 * v2) / v3;
            }
        }
        if k != n - 1 {
            a[2] = 5.0 * v3;
        }
        if k < n.saturating_sub(2) {
            a[3] = v4;
        }

        // Boundary overrides (SAMMY source lines 446-450).
        if k == n.saturating_sub(2) {
            a[5] = 0.0;
        }
        if k == n - 1 {
            a[4] = 0.0;
            a[5] = 0.0;
        }

        weights[k] = a.iter().sum::<f64>();
    }

    weights
}

/// Complementary error function erfc(x) = 1 − erf(x).
///
/// erfc(x) = exp(-x²) · exerfc(x) / √π for x ≥ 0, and
/// erfc(x) = 2 − erfc(−x) for x < 0.
pub(crate) fn erfc(x: f64) -> f64 {
    const SQRT_PI: f64 = 1.772_453_850_905_516;
    if x >= 0.0 {
        (-x * x).exp() * exerfc(x) / SQRT_PI
    } else {
        let xp = -x;
        2.0 - (-xp * xp).exp() * exerfc(xp) / SQRT_PI
    }
}

// ─── Scaled complementary error function ───────────────────────────────────────

/// Compute exp(x²)·erfc(x)·√π, numerically stable for all x.
///
/// SAMMY Ref: `fnc/exerfc.f90`.
///
/// Uses rational approximation for 0 < x ≤ 5.01 and asymptotic expansion
/// (Abramowitz & Stegun 7.1.23) for x > 5.01. For x < 0 it uses
/// erfc(x) = 2 − erfc(−x): exerfc(x) = 2√π·exp(x²) − exerfc(−x).
pub(crate) fn exerfc(x: f64) -> f64 {
    const SQRT_PI: f64 = 1.772_453_850_905_516;
    const TWO_SQRT_PI: f64 = 3.544_907_701_811_032;
    const XMAX: f64 = 5.01;
    // Rational approximation coefficients (from SAMMY's exerfc.f90)
    const A1: f64 = 8.584_076_57e-1;
    const A2: f64 = 3.078_181_93e-1;
    const A3: f64 = 6.383_238_91e-2;
    const A4: f64 = 1.824_050_75e-4;
    const A5: f64 = 6.509_742_65e-1;
    const A6: f64 = 2.294_848_19e-1;
    const A7: f64 = 3.403_018_23e-2;

    if x < 0.0 {
        TWO_SQRT_PI * (x * x).exp() - exerfc(-x)
    } else if x > XMAX {
        asympt(x)
    } else if x > 0.0 {
        let a = (A1 + x * (A2 + x * (A3 - x * A4))) / (1.0 + x * (A5 + x * (A6 + x * A7)));
        let b = SQRT_PI + x * (2.0 - a);
        b / (x * b + 1.0)
    } else {
        SQRT_PI
    }
}

/// Asymptotic expansion of exp(x²)·erfc(x)·√π for large positive x.
///
/// SAMMY Ref: `fnc/exerfc.f90`, Asympt function.
/// Uses Abramowitz & Stegun 7.1.23.
fn asympt(x: f64) -> f64 {
    if x == 0.0 {
        return 0.0;
    }
    let e = 1.0 / x;
    if e == 0.0 {
        return 0.0;
    }
    let b = 1.0 / (x * x);
    let mut a = 1.0;
    let mut c = b * 0.5;
    for n in 1..=40 {
        a -= c;
        c *= -(n as f64 + 0.5) * b;
        if (a - c) == a || (c / a).abs() < 1e-8 {
            break;
        }
    }
    a * e
}

/// Compute the Gaussian+exponential combined kernel weight Z(A, B).
///
/// Returns √π · exp(-A² + B²) · erfc(B), computed via exerfc for stability.
///
/// SAMMY Ref: `rsl/mrsl1.f90` lines 467-484 (Resbrd, Iesopr=3 path).
///
/// When B >= 0: `Z = exp(-A²) · Exerfc(B)`
/// When B < 0:  `Z = 2√π · exp(B² − A²) − exp(−A²) · Exerfc(−B)`, from
///   erfc(B) = 2 − erfc(−B), with exp(B² − A²) taken as one exponential so
///   that it does not underflow where exp(−A²) does.
///
/// SAMMY Ref: `fnc/xxerfc.f90`.
fn gauss_exp_kernel(a: f64, b: f64) -> f64 {
    const TWO_SQRT_PI: f64 = 3.544_907_701_811_032;
    let exp_neg_a2 = (-a * a).exp();
    if b >= 0.0 {
        if exp_neg_a2 == 0.0 {
            return 0.0;
        }
        exp_neg_a2 * exerfc(b)
    } else {
        TWO_SQRT_PI * (b * b - a * a).exp() - exp_neg_a2 * exerfc(-b)
    }
}

/// Compute the energy shift for the Gaussian+exponential kernel peak.
///
/// Finds the peak of the combined kernel relative to E=0 via Newton-Raphson
/// iteration. This centers the convolution window on the kernel maximum.
///
/// SAMMY Ref: `rsl/mrsl5.f90`, Shftge function.
///
/// # Arguments
/// * `c` — Mixing parameter: Widgau / (2·Widexp)
/// * `widgau` — Gaussian resolution width (eV)
///
/// # Returns
/// The energy shift Est (eV) to apply to the measurement energy.
fn shftge(c: f64, widgau: f64) -> f64 {
    const ONE_OVER_SQRT_PI: f64 = 0.564_189_583_547_756_3;
    const SMALL: f64 = 0.01;

    let ax = c;
    let bx = widgau;

    // Initial guess
    let mut x0 = if ax > ONE_OVER_SQRT_PI { ax } else { 0.0 };

    let f0_initial = ax * exerfc(x0) - 1.0;
    let mut f0 = f0_initial;
    let fff = f0;

    for _iter in 0..100 {
        let f = ax * exerfc(x0) - 1.0;
        let xma = x0 - ax;
        let q = 1.0 - 2.0 * x0 * xma;
        let delx = if q.abs() < NEAR_ZERO_FLOOR {
            // q ≈ 0: division would overflow; accept current estimate.
            break;
        } else if xma * xma - q * f > 0.0 {
            let disc = (xma * xma - q * f).sqrt();
            if xma > 0.0 {
                (-xma + disc) / q
            } else {
                (-xma - disc) / q
            }
        } else {
            if xma.abs() < NEAR_ZERO_FLOOR {
                break;
            }
            -f * 0.5 / xma
        };
        let x1 = x0 + delx;
        let shftg = (ax - x1) * bx;
        if (x1 - x0).abs() / x1.abs().max(1.0) < SMALL
            && fff.abs() > NEAR_ZERO_FLOOR
            && (f - f0).abs() / fff.abs() < SMALL
        {
            return shftg;
        }
        f0 = f;
        x0 = x1;
    }

    (ax - x0) * bx
}

/// Threshold for the ratio C = W_g / (2·W_e) above which the exponential
/// tail is negligible and the pure Gaussian PW-linear path is used instead.
///
/// At C = 2.5, erfc(2.5) ≈ 0.0005, so the exp tail contributes <0.05% of the
/// kernel integral.  Using the pure Gaussian path at this threshold introduces
/// negligible systematic error while enabling the more accurate PW-linear
/// integration and adaptive intermediate point insertion.
const EXP_TAIL_NEGLIGIBLE_C: f64 = 2.5;

/// Resolution broadening assuming the energy grid is already validated
/// (sorted ascending, same length as cross_sections).
///
/// For each broadening energy, selects the optimal integration method:
/// - **PW-linear Gaussian** (exact, second-order): when `delta_e == 0` or
///   the ratio C = W_g/(2·W_e) > [`EXP_TAIL_NEGLIGIBLE_C`] (exp tail negligible).
/// - **Combined Gaussian+exp kernel** with SAMMY Xcoef quadrature: when the
///   exponential tail is significant (C ≤ threshold).
///
/// SAMMY Ref: `rsl/mrsl1.f90` Resbrd, `convolution/DopplerAndResolutionBroadener.cpp`
pub(crate) fn resolution_broaden_presorted(
    energies: &[f64],
    cross_sections: &[f64],
    params: &ResolutionParams,
) -> Vec<f64> {
    let n = energies.len();
    if n == 0 {
        return vec![];
    }

    // Precompute Xcoef weights (used only by the combined kernel path).
    // Even if some energies take the PW-linear path, we compute weights for
    // the full grid — cheaper than branching per-energy.
    let xcoef = if params.has_exponential_tail() {
        compute_xcoef_weights(energies)
    } else {
        vec![]
    };
    let n_sigma = 5.0; // Integrate out to 5σ for Gaussian
    let mut broadened = vec![0.0f64; n];

    for i in 0..n {
        let e = energies[i];
        let widgau = params.gaussian_width(e);

        if widgau < NEAR_ZERO_FLOOR {
            broadened[i] = cross_sections[i];
            continue;
        }

        // Per-energy decision: use combined kernel only when the exp tail
        // is significant at THIS energy.
        let widexp = params.exp_width(e);
        let use_combined =
            widexp > NEAR_ZERO_FLOOR && widgau / (2.0 * widexp) <= EXP_TAIL_NEGLIGIBLE_C;

        // Compute integration limits.
        let (e_low, e_high) = if use_combined {
            // SAMMY Ref: mrsl4.f90 lines 57-65
            let wlow = n_sigma * widgau;
            let rwid = widgau / widexp;
            let wup = if rwid <= 1.0 {
                6.25 * widexp
            } else if rwid <= 2.0 {
                n_sigma * (3.0 - rwid) * widgau
            } else {
                n_sigma * widgau
            };
            (e - wlow, e + wup)
        } else {
            (e - n_sigma * widgau, e + n_sigma * widgau)
        };

        let j_lo = energies.partition_point(|&ej| ej < e_low);
        let j_hi = energies.partition_point(|&ej| ej <= e_high);

        if j_hi.saturating_sub(j_lo) <= 1 {
            broadened[i] = cross_sections[i];
            continue;
        }

        let mut sum = 0.0;
        let mut norm = 0.0;

        if use_combined {
            // Combined Gaussian + exponential kernel (SAMMY Iesopr=3)
            // with 4-point Xcoef quadrature weights.
            // SAMMY Ref: mrsl1.f90 lines 455-484
            let c = widgau * 0.5 / widexp;
            let est = shftge(c, widgau);
            let y = c * widgau + e - est;

            for j in j_lo..j_hi {
                let ee = energies[j];
                let a = (e - est - ee) / widgau;
                let b = (y - ee) / widgau;
                let z = gauss_exp_kernel(a, b);
                let wt = xcoef[j] * z;
                sum += wt * cross_sections[j];
                norm += wt;
            }
        } else {
            // Pure Gaussian kernel with piecewise-linear exact integration.
            //
            // For each interval [E_j, E_{j+1}], integrate G(E_i - E') × σ_linear(E')
            // exactly, where G(x) = exp(-x²/W²) / (W√π).
            //
            // Substituting u = (E' - E_i)/W, dE' = W du:
            //   ∫ G × [σ_j + slope×(E'-E_j)] dE'
            //   = (1/√π) ∫ exp(-u²) [σ_j + slope×W×(u - a_j)] du
            //
            // With I₀ = erf(a_{j+1}) - erf(a_j) and
            //      I₁ = (exp(-a_j²) - exp(-a_{j+1}²)) / 2:
            //
            // The normalization integral is I₀/2, so after sum/norm (2 cancels):
            //   sum += σ_j × I₀ + slope × W × (2/√π × I₁ - a_j × I₀)
            //   norm += I₀
            //
            // The factor 2/√π on I₁ comes from the u·exp(-u²) integral
            // needing to match the normalization convention erf(x) = 2/√π ∫ exp(-t²) dt.
            const TWO_OVER_SQRT_PI: f64 = std::f64::consts::FRAC_2_SQRT_PI;
            let inv_w = 1.0 / widgau;
            for j in j_lo..j_hi.saturating_sub(1) {
                let e_j = energies[j];
                let e_j1 = energies[j + 1];
                let h = e_j1 - e_j;
                if h < NEAR_ZERO_FLOOR {
                    continue;
                }

                let a_j = (e_j - e) * inv_w;
                let a_j1 = (e_j1 - e) * inv_w;

                // I₀ = erf(a_{j+1}) - erf(a_j) = erfc(a_j) - erfc(a_{j+1})
                let erfc_aj = erfc(a_j);
                let erfc_aj1 = erfc(a_j1);
                let i0 = erfc_aj - erfc_aj1;

                if i0 < NEAR_ZERO_FLOOR {
                    continue;
                }

                // I₁ = (exp(-a_j²) - exp(-a_{j+1}²)) / 2
                let i1 = ((-a_j * a_j).exp() - (-a_j1 * a_j1).exp()) * 0.5;

                let slope = (cross_sections[j + 1] - cross_sections[j]) / h;

                // σ_j × I₀ + slope × W × (2/√π × I₁ - a_j × I₀)
                sum += cross_sections[j] * i0 + slope * widgau * (TWO_OVER_SQRT_PI * i1 - a_j * i0);
                norm += i0;
            }
        }

        if norm > DIVISION_FLOOR {
            broadened[i] = sum / norm;
        } else {
            broadened[i] = cross_sections[i];
        }
    }

    broadened
}

/// A tabulated resolution function from Monte Carlo instrument simulation.
///
/// Contains reference kernels R(Δt; E_ref) at discrete energies, stored in
/// TOF-offset space (μs). Kernels are interpolated between reference energies
/// and converted from TOF to energy space when applied.
///
/// ## Offset orientation
///
/// Positive `Δt` = delayed emission (the moderator storage tail); the
/// kernel mode sits at `Δt = 0`. At apply time the broadener gathers
/// theory at `t − Δt` (convolution — see [`Self::broaden`]), so the
/// positive-`Δt` tail reads theory from earlier TOF = higher energy and
/// broadened dips acquire their tail toward lower apparent energy.
///
/// ## File Format (VENUS/FTS)
///
/// ```text
/// FTS BL10 case i00dd folded triang FWHM 350 ns PSR   ← header
/// -----                                                 ← separator
///    5.00000e-004   0.00000e+000                        ← energy block start
/// -53.458917835671329 2.051764258257523e-04             ← (tof_offset_μs, weight)
/// ...
///                                                       ← blank line separates blocks
///    1.00000e-003   0.00000e+000                        ← next energy block
/// ...
/// ```
#[derive(Debug, Clone)]
pub struct TabulatedResolution {
    /// Reference energies (eV), sorted ascending.
    ///
    /// Shared: a fit with a free `L_scale` rebinds the flight path on every
    /// forward evaluation, which must cost a reference count, not a copy.
    ref_energies: Arc<Vec<f64>>,
    /// For each reference energy: (tof_offsets_μs, weights) pairs.
    /// Weights are peak-normalized (max=1.0).
    ///
    /// Shared for the same reason as `ref_energies`.
    kernels: Arc<Vec<(Vec<f64>, Vec<f64>)>>,
    /// Flight path length in meters (needed for TOF↔energy conversion).
    flight_path_m: f64,
}

/// Trapezoidal-weighted centroid and RMS width of one kernel block.
///
/// The `dt` weights match the quadrature `broaden_presorted` integrates
/// with (single point → 1.0; edges → one-sided span; interior → half
/// the neighbour span), so these are the moments the broadener
/// effectively applies: zero-weight entries contribute nothing
/// (`tw = 0`, matching the broadener's `w <= 0` skip) and negative
/// weights are rejected at construction, so the integration domains
/// coincide exactly. The centroid pass is byte-identical to the
/// was factored out. Returns `(centroid, sigma)`; `sigma` is `0.0` for
/// a single-point or zero-mass block — callers treat a non-positive or
/// non-finite `sigma` as degenerate.
fn trapezoidal_moments(offsets: &[f64], weights: &[f64]) -> (f64, f64) {
    let n_k = offsets.len();
    let dt_width = |k: usize| -> f64 {
        if n_k <= 1 {
            1.0
        } else if k == 0 {
            offsets[1] - offsets[0]
        } else if k == n_k - 1 {
            offsets[k] - offsets[k - 1]
        } else {
            (offsets[k + 1] - offsets[k - 1]) * 0.5
        }
    };
    let (mut cnum, mut cden) = (0.0, 0.0);
    for (k, (&o, &w)) in offsets.iter().zip(weights).enumerate() {
        let tw = w * dt_width(k).abs();
        cnum += o * tw;
        cden += tw;
    }
    let centroid = if cden > 0.0 { cnum / cden } else { 0.0 };
    let mut m2 = 0.0;
    for (k, (&o, &w)) in offsets.iter().zip(weights).enumerate() {
        let tw = w * dt_width(k).abs();
        m2 += (o - centroid).powi(2) * tw;
    }
    let sigma = if cden > 0.0 { (m2 / cden).sqrt() } else { 0.0 };
    (centroid, sigma)
}

/// Integrate a sampled probability distribution into adjacent requested
/// edges. A one-point kernel is treated as a delta mass at that point; longer
/// kernels are piecewise-linear densities.
///
/// The density is normalized over its complete supplied support. The result
/// is deliberately not renormalized to `edges`: a requested detector window
/// that covers only part of the supplied pulse therefore sums to less than
/// one.
pub(crate) fn piecewise_linear_bin_probabilities(
    times: &[f64],
    weights: &[f64],
    edges: &[f64],
) -> Option<Vec<f64>> {
    if times.is_empty() || times.len() != weights.len() {
        return None;
    }
    // A duplicated (or decreasing) time gives a zero-width segment whose
    // slope is ±inf/NaN inside the CDF interpolation when an edge lands
    // there — and a NaN time passes a pure monotonicity check (every NaN
    // comparison is false) only to underflow `hi - 1` in the CDF closure
    // (`partition_point` returns 0 when the first predicate is false).
    // Validate finiteness and strict ordering of BOTH coordinate arrays up
    // front; NaN weights are caught downstream by the total-mass check.
    if times.iter().any(|t| !t.is_finite()) || times.windows(2).any(|w| w[1] <= w[0]) {
        return None;
    }
    if edges.iter().any(|e| !e.is_finite()) || edges.windows(2).any(|w| w[1] <= w[0]) {
        return None;
    }

    if times.len() == 1 {
        if !times[0].is_finite() || !weights[0].is_finite() || weights[0] <= 0.0 {
            return None;
        }
        let mut probabilities = vec![0.0; edges.len().saturating_sub(1)];
        for (index, edge) in edges.windows(2).enumerate() {
            let in_bin = edge[0] <= times[0]
                && (times[0] < edge[1]
                    || (index + 1 == probabilities.len() && times[0] == edge[1]));
            if in_bin {
                probabilities[index] = 1.0;
                break;
            }
        }
        return Some(probabilities);
    }

    let mut cumulative = Vec::with_capacity(times.len());
    cumulative.push(0.0);
    for i in 0..times.len() - 1 {
        let width = times[i + 1] - times[i];
        let area = 0.5 * (weights[i] + weights[i + 1]) * width;
        cumulative.push(cumulative[i] + area.max(0.0));
    }
    let total = *cumulative.last()?;
    if !total.is_finite() || total <= 0.0 {
        return None;
    }
    let cdf = |x: f64| -> f64 {
        if x <= times[0] {
            return 0.0;
        }
        if x >= times[times.len() - 1] {
            return 1.0;
        }
        let hi = times.partition_point(|&time| time <= x);
        let lo = hi - 1;
        let width = times[hi] - times[lo];
        let dx = x - times[lo];
        let slope = (weights[hi] - weights[lo]) / width;
        let partial = weights[lo] * dx + 0.5 * slope * dx * dx;
        ((cumulative[lo] + partial) / total).clamp(0.0, 1.0)
    };

    Some(
        edges
            .windows(2)
            .map(|edge| (cdf(edge[1]) - cdf(edge[0])).max(0.0))
            .collect(),
    )
}

impl TabulatedResolution {
    /// Reference energies (eV), sorted ascending.
    pub fn ref_energies(&self) -> &[f64] {
        &self.ref_energies
    }

    /// For each reference energy: (tof_offsets_μs, weights) pairs.
    /// Weights are peak-normalized (max=1.0).
    pub fn kernels(&self) -> &[(Vec<f64>, Vec<f64>)] {
        &self.kernels
    }

    /// Flight path length in meters (needed for TOF↔energy conversion).
    pub fn flight_path_m(&self) -> f64 {
        self.flight_path_m
    }

    /// Probability that a neutron of known true energy is recorded in each
    /// supplied detector-time bin.
    ///
    /// The tabulated pulse is selected and interpolated at `true_energy_ev`;
    /// an energy outside the tabulated reference range uses the nearest
    /// reference kernel unchanged (the same clamping the broadening path
    /// applies).  Its offsets are relative to the reference pulse mode, so
    /// the nominal arrival is
    /// `timing_offset_us + TOF_FACTOR * flight_path_m / sqrt(E)`.
    /// `timing_offset_us` is the effective clock/energy-axis offset calibrated
    /// for the measurement; this method does not invent an absolute moderator
    /// emission time that is absent from a mode-centred UDR file.
    ///
    /// The returned vector has one entry per adjacent edge pair and is not
    /// renormalized to the supplied window.  Probability outside the measured
    /// window remains outside it; the quantified acquisition-window loss at
    /// this energy is one minus the sum of the returned vector.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] unless the true energy
    /// is positive and finite, the timing offset is finite, and at least two
    /// finite bin edges are supplied in strictly increasing order.  It also
    /// fails if the interpolated tabulated pulse has zero area.
    pub fn detector_bin_probabilities(
        &self,
        true_energy_ev: f64,
        detector_time_edges_us: &[f64],
        timing_offset_us: f64,
    ) -> Result<Vec<f64>, ResolutionParseError> {
        if !true_energy_ev.is_finite() || true_energy_ev <= 0.0 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "true energy must be positive and finite, got {true_energy_ev}"
            )));
        }
        if !timing_offset_us.is_finite() {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "timing_offset_us must be finite, got {timing_offset_us}"
            )));
        }
        if detector_time_edges_us.len() < 2
            || detector_time_edges_us.iter().any(|edge| !edge.is_finite())
            || detector_time_edges_us
                .windows(2)
                .any(|edge| edge[0] >= edge[1])
        {
            return Err(ResolutionParseError::InvalidFormat(
                "detector time edges must contain at least two finite, strictly increasing values"
                    .to_string(),
            ));
        }

        let nominal_arrival =
            timing_offset_us + TOF_FACTOR * self.flight_path_m / true_energy_ev.sqrt();
        let relative_edges: Vec<f64> = detector_time_edges_us
            .iter()
            .map(|edge| edge - nominal_arrival)
            .collect();
        let (times, weights) = self.interpolated_kernel(true_energy_ev);
        piecewise_linear_bin_probabilities(&times, &weights, &relative_edges).ok_or_else(|| {
            ResolutionParseError::InvalidFormat(format!(
                "tabulated resolution at E = {true_energy_ev} eV has zero sampled \
                 area or a degenerate interpolated kernel"
            ))
        })
    }

    /// Width-corrected copy of this tabulated kernel.
    ///
    /// Shape-preserving instrument-resolution calibration knob: each
    /// reference-energy block's TOF offsets are scaled by
    /// `s(E) = s0 · (E / e_ref)^p` **about the block's intensity centroid**, so
    /// the kernel widens/narrows without moving its centroid — width and position
    /// stay orthogonal (`t0`/`L` handle absolute position). Weights are unchanged;
    /// the apply-time trapezoidal renormalization preserves unit area.
    ///
    /// Exactness note: the orthogonality is exact **at reference
    /// energies**. Between references, `interpolated_kernel`'s
    /// width-normalized blend re-scales each block about the mode
    /// (offset 0), so the applied centroid picks up a second-order
    /// dependence on the width exponent `p` (measured ~1 % of σ for
    /// |p| ≤ 0.1 on widely spaced references) — absorbed by the
    /// jointly fitted `t0` in calibration.
    ///
    /// The pivot is the **trapezoidal-weighted** centroid `Σ o·w·dt / Σ w·dt`,
    /// using the *same* `dt` quadrature weights as the broadening integral (see
    /// [`Self::broaden`]). Because `dt` is itself affine in the offsets, the width
    /// scale multiplies every `dt` by `s`, so the integrated centroid is preserved
    /// exactly on **any** offset grid (uniform or not) — not just on uniform grids
    /// where the trapezoidal and plain centroids happen to coincide.
    ///
    /// `s0 = 1, p = 0` returns a width-identical copy. This is the fittable model
    /// behind the `udr_corr` resolution-calibration family: it trusts the
    /// Monte-Carlo *shape* and calibrates only its width / energy-dependence.
    ///
    /// # Errors
    /// Returns [`ResolutionError::InvalidWidthCorrection`] unless `s0` is finite
    /// and `> 0`, `e_ref` is finite and `> 0`, and `p` is finite. A non-positive
    /// `s0` would reverse/collapse the (ascending) offset ordering the broadening
    /// loop assumes, so it is rejected up front rather than silently clamped.
    pub fn width_corrected(
        &self,
        s0: f64,
        p: f64,
        e_ref: f64,
    ) -> Result<TabulatedResolution, ResolutionError> {
        if !(s0.is_finite() && s0 > 0.0 && e_ref.is_finite() && e_ref > 0.0 && p.is_finite()) {
            return Err(ResolutionError::InvalidWidthCorrection { s0, p, e_ref });
        }
        let kernels = self
            .ref_energies
            .iter()
            .zip(self.kernels.iter())
            .map(|(&e, (offsets, weights))| {
                // The power law can overflow `s` to ±∞ for finite-but-extreme `p`;
                // reject up front rather than build a non-finite kernel directly
                // (which would bypass `from_kernels`' finiteness check).
                let s = s0 * (e / e_ref).powf(p);
                if !(s.is_finite() && s > 0.0) {
                    return Err(ResolutionError::InvalidWidthCorrection { s0, p, e_ref });
                }
                // Pivot about the trapezoidal-weighted centroid (matching the
                // `dt`-weighting in `broaden_presorted`), so the *integrated*
                // centroid is preserved on non-uniform offset grids — not only on
                // uniform grids where this reduces to the plain centroid.
                let (centroid, _) = trapezoidal_moments(offsets, weights);
                let scaled = offsets
                    .iter()
                    .map(|&o| centroid + s * (o - centroid))
                    .collect();
                Ok((scaled, weights.clone()))
            })
            .collect::<Result<Vec<_>, ResolutionError>>()?;
        Ok(TabulatedResolution {
            ref_energies: self.ref_energies.clone(),
            kernels: Arc::new(kernels),
            flight_path_m: self.flight_path_m,
        })
    }

    /// The same kernel table read against a different flight path.
    ///
    /// The stored offsets are times relative to this table's own anchor, and
    /// the flight path enters only the TOF↔energy map they are applied
    /// through. So rebinding it is exact and needs no resynthesis — the kernel
    /// itself is a property of the moderator, not of how far the neutron then
    /// flew.
    ///
    /// This is what a fitted `L_scale` requires: the data's energy grid is
    /// built with `L·L_scale`, and a kernel still reading `L` applies a width
    /// wrong by that same factor.
    ///
    /// # 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> {
        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}"
            )));
        }
        Ok(Self {
            ref_energies: Arc::clone(&self.ref_energies),
            kernels: Arc::clone(&self.kernels),
            flight_path_m,
        })
    }

    /// Kernel support at energy `e_ev`, in eV: the larger of the two distances
    /// in [`Self::gather_bounds_ev`].  Returns `0.0` for non-positive or
    /// non-finite `e_ev`, an empty kernel set, or a non-positive flight path.
    #[must_use]
    pub fn kernel_support_ev(&self, e_ev: f64) -> f64 {
        let (lo, hi) = self.gather_bounds_ev(e_ev);
        (e_ev - lo).max(hi - e_ev).max(0.0)
    }

    /// The lowest and highest energies the kernel at `e_ev` gathers theory
    /// from, as `(low, high)` in eV with `low ≤ e_ev ≤ high`, over the points
    /// [`Self::broaden`] keeps and in its arithmetic.  Returns `(e_ev, e_ev)`
    /// for non-positive or non-finite `e_ev`, an empty kernel set, or a
    /// non-positive flight path.
    #[must_use]
    pub fn gather_bounds_ev(&self, e_ev: f64) -> (f64, f64) {
        if e_ev <= 0.0 || !e_ev.is_finite() || self.kernels.is_empty() || self.flight_path_m <= 0.0
        {
            return (e_ev, e_ev);
        }
        let tof_center = TOF_FACTOR * self.flight_path_m / e_ev.sqrt();
        let (offsets, weights) = self.interpolated_kernel(e_ev);
        let (mut low, mut high) = (e_ev, e_ev);
        for (&dt, &w) in offsets.iter().zip(weights.iter()) {
            if w <= 0.0 {
                continue;
            }
            let tof_prime = tof_center - dt;
            if tof_prime <= 0.0 {
                continue;
            }
            let e_prime = (TOF_FACTOR * self.flight_path_m / tof_prime).powi(2);
            low = low.min(e_prime);
            high = high.max(e_prime);
        }
        (low, high)
    }
}

/// Resolution function: analytical Gaussian, tabulated from Monte Carlo, or
/// analytical Ikeda–Carpenter moderator model.
///
/// The `Tabulated` and `IkedaCarpenter` variants wrap an `Arc` so that cloning
/// (e.g., per-pixel in spatial mapping) is a cheap reference-count bump rather
/// than a deep copy.
///
/// `IkedaCarpenter` synthesizes a [`TabulatedResolution`] at construction and
/// is applied through the *same* per-call convolution path as `Tabulated`
/// (`broaden` / `broaden_presorted` / `plan`) — only the kernel *source* differs
/// (analytic IC pulse vs Monte-Carlo file). This keeps the three-way resolution
/// cross-validation (Gaussian | tabulated-UDR | Ikeda–Carpenter) fair on the
/// reference broadening path. Note: `IkedaCarpenter` does **not** opt into the
/// spatial-map surrogate fast-paths (the scalar/cubature plans gate on
/// `Tabulated`); it falls back to the general path, which is correct but
/// unoptimized — see the resolution-calibration notes for the W6 follow-up.
#[derive(Debug, Clone)]
pub enum ResolutionFunction {
    /// Analytical Gaussian resolution from instrument parameters.
    Gaussian(ResolutionParams),
    /// Tabulated resolution from Monte Carlo instrument simulation.
    Tabulated(Arc<TabulatedResolution>),
    /// Analytical Ikeda–Carpenter moderator resolution model.
    IkedaCarpenter(Arc<crate::ikeda_carpenter::IkedaCarpenter>),
}

/// Widths of Gaussian resolution the broadening limits reach on each side.
///
/// SAMMY Ref: `rsl/mrsl4.f90` `Wdsint`, `Wlow = Wup = Brdlim*Widgau`;
/// `inp/minp06.f` line 212, `Brdlim = 5`.
const BRDLIM: f64 = 5.0;

/// Lowest energy the Gaussian working grid extends to, in eV: the PW-linear
/// quadrature maps points through `1/√E`, which has no value at zero.
const GAUSSIAN_LOW_ENERGY_FLOOR_EV: f64 = 0.001;

impl ResolutionFunction {
    /// The energies a working grid for the data window `energies` has to span,
    /// as `(low, high)` in eV with `low ≤ e_min` and `high ≥ e_max`: SAMMY's
    /// Wdsint limits at the two ends for a Gaussian, the extremes of
    /// [`TabulatedResolution::gather_bounds_ev`] over the grid for a sampled
    /// kernel.  Returns `(0.0, 0.0)` for an empty grid.
    #[must_use]
    pub fn grid_bounds_ev(&self, energies: &[f64]) -> (f64, f64) {
        let (Some(&e_min), Some(&e_max)) = (energies.first(), energies.last()) else {
            return (0.0, 0.0);
        };
        let sampled = |table: &TabulatedResolution| {
            energies.iter().fold((e_min, e_max), |(low, high), &e| {
                let (l, h) = table.gather_bounds_ev(e);
                (low.min(l), high.max(h))
            })
        };
        match self {
            Self::Gaussian(params) => {
                let wg_high = params.gaussian_width(e_max);
                let we_high = params.exp_width(e_max);
                // SAMMY grades the high side by the ratio of the Gaussian core to
                // the one-sided exponential tail (Wdsint, `Rwid`).
                let above = if we_high > 1e-30 {
                    let rwid = wg_high / we_high;
                    if rwid <= 1.0 {
                        6.25 * we_high
                    } else if rwid <= 2.0 {
                        BRDLIM * (3.0 - rwid) * wg_high
                    } else {
                        BRDLIM * wg_high
                    }
                } else {
                    BRDLIM * wg_high
                };
                let below = (BRDLIM * params.gaussian_width(e_min))
                    .min(e_min - GAUSSIAN_LOW_ENERGY_FLOOR_EV)
                    .max(0.0);
                (e_min - below, e_max + above)
            }
            Self::Tabulated(tabulated) => sampled(tabulated),
            Self::IkedaCarpenter(ic) => sampled(ic.tabulated()),
        }
    }

    /// Flight path used to map true neutron energy to detector time.
    pub fn flight_path_m(&self) -> f64 {
        match self {
            Self::Gaussian(params) => params.flight_path_m(),
            Self::Tabulated(tabulated) => tabulated.flight_path_m(),
            Self::IkedaCarpenter(ic) => ic.flight_path_m(),
        }
    }

    /// The same resolution function read against a different flight path.
    ///
    /// A fit that frees `L_scale` evaluates the theory on an energy grid built
    /// with `L·L_scale`. The kernel has to be read against the same flight
    /// path or its width is wrong by that factor — on every family, since all
    /// three convert between energy and detector time through `L`. No variant
    /// resynthesizes: the flight path is not part of what the kernel IS, only
    /// of the map it is applied through.
    ///
    /// # 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> {
        match self {
            Self::Gaussian(params) => Ok(Self::Gaussian(
                ResolutionParams::new(
                    flight_path_m,
                    params.delta_t_us(),
                    params.delta_l_m(),
                    params.delta_e_us(),
                )
                .map_err(|e| ResolutionParseError::InvalidFormat(e.to_string()))?,
            )),
            Self::Tabulated(tabulated) => Ok(Self::Tabulated(Arc::new(
                tabulated.with_flight_path(flight_path_m)?,
            ))),
            Self::IkedaCarpenter(ic) => Ok(Self::IkedaCarpenter(Arc::new(
                ic.with_flight_path(flight_path_m)?,
            ))),
        }
    }

    /// Probability that one neutron of known true energy is recorded in each
    /// supplied detector-time bin.
    ///
    /// The tabulated and Ikeda–Carpenter variants are evaluated directly in
    /// detector time. In particular, the analytical IC variant does not pass
    /// through its legacy synthesized [`TabulatedResolution`] broadening
    /// table, and it refuses a Gaussian burst. The older Gaussian
    /// energy-broadening model has no physical detector-time probability law
    /// and is therefore rejected rather than silently treated as one.
    ///
    /// `timing_offset_us` is convention-dependent and NOT transferable
    /// between variants: a mode-centred tabulated (UDR) kernel places its
    /// pulse mode at the nominal arrival, so the offset must absorb the
    /// calibrated moderator mean delay, while the causal Ikeda–Carpenter
    /// pulse rises from the nominal arrival onward and its offset is a pure
    /// clock/detector shift.  Swapping response models under one calibrated
    /// offset shifts every bin systematically.
    pub fn detector_bin_probabilities(
        &self,
        true_energy_ev: f64,
        detector_time_edges_us: &[f64],
        timing_offset_us: f64,
    ) -> Result<Vec<f64>, ResolutionParseError> {
        match self {
            Self::Tabulated(tabulated) => tabulated.detector_bin_probabilities(
                true_energy_ev,
                detector_time_edges_us,
                timing_offset_us,
            ),
            Self::IkedaCarpenter(ic) => ic.detector_bin_probabilities(
                true_energy_ev,
                detector_time_edges_us,
                timing_offset_us,
            ),
            Self::Gaussian(_) => Err(ResolutionParseError::InvalidFormat(
                "Gaussian energy broadening cannot produce detector-time bin probabilities; use a validated tabulated or Ikeda–Carpenter time response"
                    .to_string(),
            )),
        }
    }
}

/// Pre-built resolution-broadening plan for a specific target energy grid.
///
/// Encodes every quantity that depends only on the target grid, the
/// reference kernel, and the flight path — so applying the plan to a
/// spectrum reduces to a gather + multiply-add loop with no
/// transcendentals, no allocations, and no binary / pointer search.
///
/// Build via [`TabulatedResolution::plan`] — returns a `Result` and
/// validates the sorted-grid precondition that `broaden` enforces.
/// Apply via [`ResolutionPlan::apply`].  One plan is tied to one
/// `(target_energies, ref_energies, flight_path_m)` triple; the plan
/// owns a copy of the target-energy grid so callers cannot apply it to
/// a spectrum that was measured on a *different* grid even when the
/// grid length matches — use [`Self::target_energies`] to verify the
/// grid identity before applying.
///
/// The layout is a flat Struct-of-Arrays (SoA): per-target `(lo_idx,
/// frac, weight)` tuples packed into three parallel `Vec`s, with
/// `starts[i]..starts[i+1]` naming the range for target `i`.  SoA keeps
/// the inner loop memory-access pattern sequential and cache-friendly.
#[derive(Debug, Clone)]
pub struct ResolutionPlan {
    /// Target energy grid the plan was built for (owned copy).
    ///
    /// Stored so `apply()` can verify `spectrum.len() == self.len()`
    /// and expose a cheap grid identity for caller-side caching.
    /// ~28 KB for the VENUS 3471-point grid — negligible compared to
    /// the ~8 MB `lo_idx`/`frac`/`weight` footprint of a full plan.
    target_energies: Vec<f64>,
    /// `starts[i]..starts[i+1]` indexes into `lo_idx`/`frac`/`weight`
    /// for target `i`.  `starts` has length `target_energies.len() + 1`.
    starts: Vec<u32>,
    /// For each valid (target, kernel-point) entry: the lower bracket
    /// index into the target grid (spectrum[lo] + frac * (spectrum[lo+1]
    /// - spectrum[lo])).
    lo_idx: Vec<u32>,
    /// Spectrum-interp fraction in [0, 1].  Set to 0 for degenerate
    /// brackets; the apply-time loop short-circuits `frac == 0.0` so
    /// degenerate entries never touch `spectrum[lo+1]`.  This matches
    /// `broaden_presorted` even when `spectrum[lo+1]` is NaN/±∞.
    frac: Vec<f64>,
    /// Pre-computed per-entry weight (`w * dt_width.abs()`).  Summing
    /// these yields the per-target normalisation.
    weight: Vec<f64>,
    /// Pre-summed `Σ weight` per target (in the same accumulation order
    /// as `broaden_presorted` visits the valid entries).  When `norm <=
    /// DIVISION_FLOOR` the apply path returns `spectrum[i]` directly
    /// — the exact `broaden_presorted` passthrough behaviour.
    norm: Vec<f64>,
}

impl ResolutionPlan {
    /// Number of target energies this plan covers.
    pub fn len(&self) -> usize {
        self.target_energies.len()
    }

    /// True when the plan covers no target energies.
    pub fn is_empty(&self) -> bool {
        self.target_energies.is_empty()
    }

    /// Target energy grid the plan was built for.
    ///
    /// Callers implementing plan caches can compare this against their
    /// current grid to decide whether the plan is still valid.  Using
    /// pointer identity of the returned slice gives an O(1) check when
    /// the grid hasn't moved; slice equality is `O(n)` but catches
    /// cases where the underlying buffer was reallocated.
    pub fn target_energies(&self) -> &[f64] {
        &self.target_energies
    }

    /// Apply the plan to a spectrum on the same target grid the plan
    /// was built for.
    ///
    /// The spectrum length must equal [`Self::len`].  Passing a
    /// spectrum on a different grid that happens to have the same
    /// length is caller error — verify via [`Self::target_energies`]
    /// when in doubt.
    ///
    /// Bit-exact with `broaden_presorted(target_energies, spectrum)`
    /// for finite spectrum values; degenerate-bracket entries
    /// short-circuit the interpolation so the equivalence also holds
    /// when `spectrum[lo+1]` is NaN or ±∞ (the reference path returns
    /// `spectrum[lo]` directly in that case without touching the upper
    /// bracket).
    pub fn apply(&self, spectrum: &[f64]) -> Vec<f64> {
        let n = self.target_energies.len();
        assert_eq!(
            spectrum.len(),
            n,
            "spectrum length ({}) must match plan target-grid length ({})",
            spectrum.len(),
            n,
        );
        if n == 0 {
            return Vec::new();
        }

        let mut result = vec![0.0f64; n];

        // Pre-bind plan slices once per call and pre-slice each
        // target's entry range before the hot loop.  This is a
        // bounds-check-elimination (BCE) refactor — every per-entry
        // index is proven in-bounds by the invariants established in
        // `plan_presorted`, so the inner loop uses `get_unchecked`
        // with SAFETY comments citing those invariants.  The compiler
        // then auto-vectorizes the inner compute where profitable.
        //
        // We deliberately do NOT use explicit 2-wide SIMD here — an
        // experiment via the `wide` crate (commit abandoned;
        // `perf-lessons.md`) showed that 2-wide f64x2 with gather
        // emulation is net-negative on AArch64 Neon vs the compiler's
        // scalar auto-vectorization of the BCE'd inner loop.  On
        // wider targets (x86 AVX2 / AVX-512) a SIMD rewrite could
        // still pay off but is out of scope here.
        //
        // Control flow, accumulation order, and the `frac == 0.0`
        // NaN-safety short-circuit are all preserved exactly so the
        // bit-exact contract with `broaden_presorted` holds for
        // finite AND pathological (NaN, ±∞) spectra.
        let lo_idx = self.lo_idx.as_slice();
        let frac_all = self.frac.as_slice();
        let weight_all = self.weight.as_slice();
        let starts = self.starts.as_slice();
        let norm = self.norm.as_slice();
        let spec = spectrum;

        // Defence-in-depth: debug-only invariant checks right after
        // slice binding, so a future change to `plan_presorted` that
        // silently violates the `unsafe { get_unchecked }` SAFETY
        // claims below fails loudly in debug builds.  Zero release-
        // build cost.
        debug_assert_eq!(starts.len(), n + 1);
        debug_assert_eq!(
            starts.last().copied(),
            Some(lo_idx.len() as u32),
            "plan_presorted invariant: starts.last() must equal lo_idx.len()",
        );
        debug_assert_eq!(lo_idx.len(), frac_all.len());
        debug_assert_eq!(lo_idx.len(), weight_all.len());
        debug_assert_eq!(norm.len(), n);
        debug_assert_eq!(spec.len(), n);

        for i in 0..n {
            let norm_i = norm[i];
            if norm_i <= DIVISION_FLOOR {
                // Passthrough — matches `broaden_presorted`'s
                // `spectrum[i]` fallback for e ≤ 0, empty kernel, or
                // degenerate norm accumulation.
                result[i] = spec[i];
                continue;
            }
            let start = starts[i] as usize;
            let end = starts[i + 1] as usize;
            // Zip-compatible pre-bound slices of exactly `end - start`
            // elements each — the per-j bounds check is elided by the
            // compiler because the slice length bounds the loop.
            let los = &lo_idx[start..end];
            let fracs = &frac_all[start..end];
            let ws = &weight_all[start..end];

            let mut sum = 0.0f64;
            for k in 0..los.len() {
                // SAFETY: `k < los.len()` is guaranteed by the range;
                // `los`, `fracs`, and `ws` all have length `end - start`
                // (same subslice bounds), so each `get_unchecked(k)`
                // read is in-bounds.
                let lo = unsafe { *los.get_unchecked(k) } as usize;
                let frac = unsafe { *fracs.get_unchecked(k) };
                let w = unsafe { *ws.get_unchecked(k) };

                // Degenerate-bracket short-circuit: when the plan
                // built `frac = -0.0` (span < NEAR_ZERO_FLOOR) we skip
                // `spectrum[lo+1]` entirely.  Without this branch,
                // `0.0 * NaN = NaN` would propagate and diverge from
                // the reference `broaden_presorted`, which returns
                // `spectrum[lo]` directly for that case.  Branch is
                // well-predicted (degenerate brackets are rare on
                // real grids) and preserves bit-exactness under
                // pathological spectra.
                //
                // The check MUST use `to_bits()` because the non-
                // degenerate path can legitimately produce
                // `frac == +0.0` when `e_prime == energies[lo]`
                // exactly.  In that case `broaden_presorted` still
                // reads `spectrum[lo+1]` (and propagates NaN if
                // present there), so the short-circuit MUST NOT
                // trigger.  `+0.0 == -0.0` returns `true` but
                // `(+0.0).to_bits() != (-0.0).to_bits()`, so the
                // bit-pattern check disambiguates exactly which
                // semantic `plan_presorted` meant.
                let s = if frac.to_bits() == (-0.0_f64).to_bits() {
                    // SAFETY: `lo < n` by plan invariant.
                    // `plan_presorted` only pushes `lo = bracket_hi - 1`
                    // with `bracket_hi ∈ [1, n - 1]`, so `lo ∈
                    // [0, n - 2]`.  `spec.len() == n` by the
                    // precondition assert at the top of `apply`.
                    unsafe { *spec.get_unchecked(lo) }
                } else {
                    // SAFETY: same `lo ∈ [0, n - 2]` invariant, so
                    // `lo + 1 ∈ [1, n - 1]` is also in-bounds.
                    let s_lo = unsafe { *spec.get_unchecked(lo) };
                    let s_hi = unsafe { *spec.get_unchecked(lo + 1) };
                    s_lo + frac * (s_hi - s_lo)
                };
                // Serial accumulation preserved — no multi-accumulator
                // reassociation, no SIMD lane-wise tree reduce.
                // IEEE-754 addition is not associative; changing the
                // order would break bit-exactness with
                // `broaden_presorted_reference` (and all
                // `*_bit_exact_*` unit tests + the maintainers'
                // real-VENUS bit-exact baseline harness).
                sum += w * s;
            }
            result[i] = sum / norm_i;
        }

        result
    }

    /// Compile this plan into a row-stochastic CSR
    /// [`ResolutionMatrix`].
    ///
    /// The compiled matrix is an explicit sparse representation of
    /// the resolution operator `R` on the plan's target grid.  Each
    /// row sums to 1.0 to machine precision (passthrough rows store
    /// a single `(i, i, 1.0)` entry to match [`ResolutionPlan::apply`]
    /// 's `norm ≤ DIVISION_FLOOR` fallback).
    ///
    /// Degenerate-bracket handling uses the `-0.0` sentinel
    /// convention from `plan_presorted`: if `plan.frac[e]` has the
    /// bit pattern of `-0.0`, the entry contributes `weight / norm`
    /// at column `lo` only (no `lo+1` bracket).  A regular `+0.0`
    /// frac contributes `weight * 1.0 / norm` at `lo` and
    /// `weight * 0.0 / norm = 0.0` at `lo+1` — those zero columns
    /// are retained in CSR with `value = 0.0` to preserve
    /// downstream NaN-safety if the consumer re-multiplies by a
    /// spectrum containing NaN at `lo+1`.
    ///
    /// # Equivalence contract (finite spectra only)
    ///
    /// For a spectrum with **all finite values**, [`apply_r`] on the
    /// compiled matrix produces per-element output within `1e-12`
    /// relative tolerance of [`Self::apply`] on the same spectrum —
    /// not bit-exact, because the CSR matvec sums contributions in
    /// column order while `apply` sums in entry order and IEEE-754
    /// addition is non-associative.  The `1e-12` bound accounts for
    /// accumulation error across the ~82 entries per row on the
    /// 3471-bin VENUS production grid (500 × 2.22e-16 ≈ 1.1e-13 per
    /// row; `1e-12` leaves comfortable headroom).
    ///
    /// # Non-finite and near-overflow spectra
    ///
    /// The equivalence bound does **NOT** extend to spectra with
    /// `NaN` / `±∞` values, **nor to near-f64::MAX overflow
    /// inputs**.  Both divergences trace back to the same
    /// algebraic rewrite:
    ///
    /// * [`Self::apply`] computes each entry as `spec[lo] + frac *
    ///   (spec[lo+1] - spec[lo])`, which can overflow the
    ///   subtraction even for finite inputs (opposite-sign
    ///   f64::MAX → `-∞`).
    /// * The compiled CSR form splits the interp into `(1 - frac) *
    ///   spec[lo] + frac * spec[lo + 1]`, which scales before
    ///   summing and stays finite in the same case.
    ///
    /// For bounded finite Beer-Lambert transmissions (`T ∈ [0, 1]`)
    /// neither divergence can arise; callers who deliberately pass
    /// non-finite or near-overflow spectra (e.g., as debug sentinels
    /// or out-of-range diagnostics) must not rely on cross-API
    /// equivalence.  See `resolution_matrix_nonfinite_contract` and
    /// `resolution_matrix_large_finite_contract` for executable
    /// demonstrations.
    pub fn compile_to_matrix(&self) -> ResolutionMatrix {
        let n = self.target_energies.len();
        let mut row_starts: Vec<u32> = Vec::with_capacity(n + 1);
        row_starts.push(0);
        let mut col_indices: Vec<u32> = Vec::new();
        let mut values: Vec<f64> = Vec::new();

        // Reusable per-row accumulator.  Columns accumulate into a
        // BTreeMap keyed by spectrum index so the final CSR row is
        // emitted in ascending column order — the required CSR
        // invariant and the condition the `apply_r` equivalence
        // bound depends on.
        let mut acc: std::collections::BTreeMap<u32, f64> = std::collections::BTreeMap::new();

        for i in 0..n {
            acc.clear();
            let norm_i = self.norm[i];
            if norm_i <= DIVISION_FLOOR {
                // Passthrough row — matches `apply`'s early return.
                col_indices.push(i as u32);
                values.push(1.0);
                // See u32-overflow `debug_assert!` below — the same
                // bound applies after every `push`.
                debug_assert!(
                    col_indices.len() <= u32::MAX as usize,
                    "CSR row_starts/col_indices u32 overflow: nnz = {}",
                    col_indices.len(),
                );
                row_starts.push(col_indices.len() as u32);
                continue;
            }
            let start = self.starts[i] as usize;
            let end = self.starts[i + 1] as usize;
            for e in start..end {
                let lo = self.lo_idx[e];
                let frac = self.frac[e];
                let w = self.weight[e];
                if frac.to_bits() == (-0.0_f64).to_bits() {
                    // Degenerate bracket — `apply` reads `spec[lo]`
                    // only, so the CSR row contributes only at `lo`.
                    *acc.entry(lo).or_insert(0.0) += w / norm_i;
                } else {
                    // Regular linear-interp entry: `w * ((1 - frac)
                    // * spec[lo] + frac * spec[lo + 1]) / norm_i`.
                    *acc.entry(lo).or_insert(0.0) += w * (1.0 - frac) / norm_i;
                    *acc.entry(lo + 1).or_insert(0.0) += w * frac / norm_i;
                }
            }
            for (&col, &val) in acc.iter() {
                col_indices.push(col);
                values.push(val);
            }
            // Defence-in-depth: a future large-grid caller that
            // accumulates more than u32::MAX entries would silently
            // truncate the `as u32` cast below.  The `plan_presorted`
            // helper already has matching `debug_assert!` guards on
            // its u32 offsets (resolution.rs, `plan_presorted`).
            debug_assert!(
                col_indices.len() <= u32::MAX as usize,
                "CSR row_starts/col_indices u32 overflow: nnz = {}",
                col_indices.len(),
            );
            row_starts.push(col_indices.len() as u32);
        }

        ResolutionMatrix {
            target_energies: self.target_energies.clone(),
            row_starts,
            col_indices,
            values,
        }
    }
}

/// Row-stochastic CSR representation of the resolution operator `R`
/// on a fixed target energy grid.
///
/// Built from a [`ResolutionPlan`] via
/// [`ResolutionPlan::compile_to_matrix`].  Exposed so downstream
/// surrogates (see ) can access the row-local entries
/// `R_{i, j}` directly for LP / quadrature construction.
///
/// Owns a copy of the target energy grid for the same reason
/// [`ResolutionPlan`] does: caller-side grid-identity checks and
/// explicit grid-mismatch errors via
/// [`ResolutionError::MatrixGridMismatch`].
#[derive(Debug, Clone)]
pub struct ResolutionMatrix {
    /// Target energy grid the matrix was compiled for (owned copy).
    target_energies: Vec<f64>,
    /// `row_starts[i]..row_starts[i+1]` indexes into
    /// `col_indices`/`values` for row `i`.  Length `n + 1`.
    row_starts: Vec<u32>,
    /// Column indices in ascending order within each row.
    col_indices: Vec<u32>,
    /// CSR values.  Row `i` sums to 1.0 within machine precision
    /// (passthrough rows store exactly `1.0` at column `i`).
    values: Vec<f64>,
}

impl ResolutionMatrix {
    /// Number of rows (target-grid size) covered by this matrix.
    pub fn len(&self) -> usize {
        self.target_energies.len()
    }

    /// True when the matrix covers no target energies.
    pub fn is_empty(&self) -> bool {
        self.target_energies.is_empty()
    }

    /// Total number of stored entries (structural nnz).
    ///
    /// Regular-bracket entries with `frac == +0.0` retain a
    /// zero-valued contribution at the `lo + 1` column to preserve
    /// NaN-safety under re-application to spectra with NaN at that
    /// column; those stored zeros are counted in this total.
    pub fn nnz(&self) -> usize {
        self.values.len()
    }

    /// Target energy grid the matrix was compiled for.
    pub fn target_energies(&self) -> &[f64] {
        &self.target_energies
    }

    /// CSR row-start offsets.  `row_starts()[i]..row_starts()[i+1]`
    /// names the entry range for row `i`.  Length `len() + 1`.
    pub fn row_starts(&self) -> &[u32] {
        &self.row_starts
    }

    /// CSR column indices.  Sorted ascending within each row.
    pub fn col_indices(&self) -> &[u32] {
        &self.col_indices
    }

    /// CSR values.  Each row sums to 1.0 to machine precision.
    pub fn values(&self) -> &[f64] {
        &self.values
    }
}

/// Apply a compiled [`ResolutionMatrix`] to a spectrum on the same
/// target grid the matrix was compiled for.
///
/// For finite spectra, the output is numerically equivalent to
/// [`ResolutionPlan::apply`] on the same spectrum within `1e-12`
/// relative tolerance per element; not bit-exact, because CSR matvec
/// sums in column order while `ResolutionPlan::apply` sums in entry
/// order.
///
/// # Non-finite and near-overflow inputs
///
/// See [`ResolutionPlan::compile_to_matrix`] for the full contract
/// on `NaN` / `±∞` spectra **and on near-f64::MAX finite spectra** —
/// the equivalence bound does not extend to either.  Production
/// forward models feed Beer-Lambert transmissions (`T ∈ [0, 1]`) so
/// the distinction never arises in practice.
///
/// # Panics
///
/// Panics if `spectrum.len() != matrix.len()`.  Use
/// [`apply_resolution_with_matrix`] for a checked entrypoint that
/// returns [`ResolutionError::LengthMismatch`] instead.
pub fn apply_r(matrix: &ResolutionMatrix, spectrum: &[f64]) -> Vec<f64> {
    let n = matrix.len();
    assert_eq!(
        spectrum.len(),
        n,
        "spectrum length ({}) must match matrix grid length ({})",
        spectrum.len(),
        n,
    );
    let mut out = vec![0.0f64; n];
    for (i, out_i) in out.iter_mut().enumerate() {
        let start = matrix.row_starts[i] as usize;
        let end = matrix.row_starts[i + 1] as usize;
        let mut sum = 0.0f64;
        for e in start..end {
            let col = matrix.col_indices[e] as usize;
            sum += matrix.values[e] * spectrum[col];
        }
        *out_i = sum;
    }
    out
}

/// Checked variant of [`apply_r`] that validates the matrix was
/// compiled for `energies` before applying.
///
/// Returns [`ResolutionError::LengthMismatch`] when either
/// `energies` or `spectrum` has a length that disagrees with the
/// matrix grid size.  For the `spectrum` check, the `energies` field
/// of the returned error holds the matrix grid length (the required
/// length) so callers can read it as "expected vs got".  Returns
/// [`ResolutionError::MatrixGridMismatch`] when the lengths match
/// but the grid contents differ (per-element `to_bits()` compare).
///
/// Unlike [`apply_resolution_with_plan`], this entrypoint does not
/// enforce an ascending `energies` grid through the crate's internal
/// `validate_inputs` helper.  That check is redundant here: the plan
/// that produced the matrix was itself built on a sorted grid (via
/// [`TabulatedResolution::plan`], which validates sortedness), and the
/// stored `target_energies` copy is used in the `to_bits()`
/// grid-identity check above.  Any `energies` slice that is not
/// bit-identical to the matrix's stored copy — including an unsorted
/// permutation of the same values — fails with
/// [`ResolutionError::MatrixGridMismatch`].
pub fn apply_resolution_with_matrix(
    energies: &[f64],
    matrix: &ResolutionMatrix,
    spectrum: &[f64],
) -> Result<Vec<f64>, ResolutionError> {
    if energies.len() != matrix.len() {
        return Err(ResolutionError::LengthMismatch {
            energies: energies.len(),
            data: matrix.len(),
        });
    }
    if spectrum.len() != matrix.len() {
        // Reuse the `LengthMismatch` variant for the spectrum branch:
        // `energies` = expected length (matrix grid size), `data` =
        // actual spectrum length.  See docstring above.
        return Err(ResolutionError::LengthMismatch {
            energies: matrix.len(),
            data: spectrum.len(),
        });
    }
    for (i, (e_cur, e_ref)) in energies.iter().zip(matrix.target_energies()).enumerate() {
        // `to_bits()` equality catches `-0.0 vs +0.0` and NaN-bit
        // differences that float `==` silently accepts or rejects.
        if e_cur.to_bits() != e_ref.to_bits() {
            return Err(ResolutionError::MatrixGridMismatch {
                first_diff_index: i,
            });
        }
    }
    Ok(apply_r(matrix, spectrum))
}

impl TabulatedResolution {
    /// Parse a VENUS/FTS resolution file.
    ///
    /// # Arguments
    /// * `text` — File contents as a string.
    /// * `flight_path_m` — Flight path length in meters.
    pub fn from_text(text: &str, flight_path_m: f64) -> Result<Self, ResolutionParseError> {
        if !flight_path_m.is_finite() || flight_path_m <= 0.0 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Flight path must be positive and finite, got {flight_path_m}"
            )));
        }
        let mut lines = text.lines();

        // Skip header and separator
        let _header = lines
            .next()
            .ok_or(ResolutionParseError::InvalidFormat("Empty file".into()))?;
        let _sep = lines.next().ok_or(ResolutionParseError::InvalidFormat(
            "Missing separator".into(),
        ))?;

        let mut ref_energies = Vec::new();
        let mut kernels = Vec::new();
        let mut current_energy: Option<f64> = None;
        let mut current_offsets: Vec<f64> = Vec::new();
        let mut current_weights: Vec<f64> = Vec::new();

        for line in lines {
            let trimmed = line.trim();
            if trimmed.is_empty() {
                // End of current block
                if let Some(e) = current_energy.take() {
                    ref_energies.push(e);
                    kernels.push((
                        std::mem::take(&mut current_offsets),
                        std::mem::take(&mut current_weights),
                    ));
                }
                continue;
            }

            let parts: Vec<&str> = trimmed.split_whitespace().collect();
            if parts.len() != 2 {
                if current_energy.is_some() {
                    return Err(ResolutionParseError::InvalidFormat(format!(
                        "Expected 2 columns inside energy block, got {}: '{}'",
                        parts.len(),
                        trimmed
                    )));
                }
                // Outside a data block (e.g. extra header lines) — skip
                continue;
            }

            let x: f64 = parts[0].parse().map_err(|_| {
                ResolutionParseError::InvalidFormat(format!("Cannot parse float: '{}'", parts[0]))
            })?;
            let y: f64 = parts[1].parse().map_err(|_| {
                ResolutionParseError::InvalidFormat(format!("Cannot parse float: '{}'", parts[1]))
            })?;

            if current_energy.is_none() {
                // First line of block: energy + 0.0 marker
                current_energy = Some(x);
            } else {
                current_offsets.push(x);
                current_weights.push(y);
            }
        }

        // Flush last block
        if let Some(e) = current_energy.take() {
            ref_energies.push(e);
            kernels.push((current_offsets, current_weights));
        }

        if ref_energies.is_empty() {
            return Err(ResolutionParseError::InvalidFormat(
                "No energy blocks found".into(),
            ));
        }

        // Validate finite, POSITIVE, strictly ascending reference
        // energies.  The finiteness check must come first: NaN compares
        // false against everything, so a NaN energy would slip through
        // the ascending check below and then poison the bracketing
        // binary search.  Positivity is load-bearing twice over: the
        // TOF map t = TOF_FACTOR·L/√E needs E > 0, and the
        // between-reference width interpolation takes ln(E_ref) — a
        // non-positive reference would turn every blended weight into
        // NaN, which bypasses the broadener's norm guard and silently
        // disables broadening (NaN comparisons are false).
        if let Some(bad) = ref_energies.iter().find(|e| !(e.is_finite() && **e > 0.0)) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Reference energies must be finite and positive, got {bad}"
            )));
        }
        for i in 1..ref_energies.len() {
            if ref_energies[i] <= ref_energies[i - 1] {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Reference energies must be strictly ascending, but E[{}]={} <= E[{}]={}",
                    i,
                    ref_energies[i],
                    i - 1,
                    ref_energies[i - 1],
                )));
            }
        }

        // Per-block kernel invariants — the same set `from_kernels`
        // enforces, because both constructors feed the same broadener:
        //
        // * non-empty: an empty block would make `broaden_presorted`
        //   accumulate `norm == 0` and silently pass the spectrum
        //   through;
        // * all-finite offsets and weights (checked BEFORE sortedness:
        //   a NaN offset fails the ascending comparison too, and the
        //   "must be strictly ascending" message would mislead): a NaN
        //   weight poisons the kernel norm and bypasses the division
        //   guard;
        // * strictly ascending offsets: the monotonic two-pointer
        //   bracket walk and the trapezoidal `dt` widths
        //   (`offsets[k+1] − offsets[k−1]`) both assume sorted offsets
        //   — an unsorted block would broaden silently-wrong.
        for (i, (offsets, weights)) in kernels.iter().enumerate() {
            if offsets.is_empty() {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} (E = {} eV) has no (offset, weight) points",
                    ref_energies[i],
                )));
            }
            if offsets.iter().any(|v| !v.is_finite()) || weights.iter().any(|v| !v.is_finite()) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} (E = {} eV) contains non-finite offsets or weights",
                    ref_energies[i],
                )));
            }
            // Negative weights have no physical meaning (a resolution
            // kernel is an emission-time density); the broadener skips
            // w <= 0 entries, and the width machinery's trapezoidal
            // moments must integrate over the same domain — reject at
            // the door rather than let the two disagree.
            if weights.iter().any(|&v| v < 0.0) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} (E = {} eV) contains negative weights",
                    ref_energies[i],
                )));
            }
            if !offsets.windows(2).all(|w| w[0] < w[1]) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} (E = {} eV) TOF offsets must be strictly ascending",
                    ref_energies[i],
                )));
            }
        }

        Ok(TabulatedResolution {
            ref_energies: Arc::new(ref_energies),
            kernels: Arc::new(kernels),
            flight_path_m,
        })
    }

    /// Build a tabulated resolution directly from synthesized kernels.
    ///
    /// Used by the analytical [`crate::ikeda_carpenter::IkedaCarpenter`] model,
    /// which generates `(tof_offset_µs, weight)` kernels at a set of reference
    /// energies and then rides the exact same broadening machinery as a
    /// Monte-Carlo file. Validates the same invariants `from_text` enforces:
    /// non-empty + strictly ascending reference energies, one kernel per energy,
    /// each kernel non-empty with matching offset/weight lengths, and all-finite
    /// offsets/weights.
    ///
    /// # Errors
    /// Returns [`ResolutionParseError::InvalidFormat`] if the flight path is
    /// not positive and finite, if the reference energies are empty / not
    /// strictly ascending, if the energy and kernel counts differ, if any
    /// kernel is empty, if a kernel's offset and weight vectors differ in
    /// length, or if any offset/weight is non-finite.
    pub fn from_kernels(
        ref_energies: Vec<f64>,
        kernels: Vec<(Vec<f64>, Vec<f64>)>,
        flight_path_m: f64,
    ) -> Result<Self, ResolutionParseError> {
        if !flight_path_m.is_finite() || flight_path_m <= 0.0 {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Flight path must be positive and finite, got {flight_path_m}"
            )));
        }
        if ref_energies.is_empty() {
            return Err(ResolutionParseError::InvalidFormat(
                "No reference energies provided".into(),
            ));
        }
        if ref_energies.len() != kernels.len() {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Reference-energy count {} != kernel count {}",
                ref_energies.len(),
                kernels.len(),
            )));
        }
        // Finiteness first: NaN compares false against everything, so a
        // NaN energy would slip through the ascending check below and
        // then poison the bracketing binary search.  Positivity is
        // load-bearing twice over: the TOF map needs E > 0, and the
        // between-reference width interpolation takes ln(E_ref) — a
        // non-positive reference would turn every blended weight into
        // NaN and silently disable broadening.
        if let Some(bad) = ref_energies.iter().find(|e| !(e.is_finite() && **e > 0.0)) {
            return Err(ResolutionParseError::InvalidFormat(format!(
                "Reference energies must be finite and positive, got {bad}"
            )));
        }
        for i in 1..ref_energies.len() {
            if ref_energies[i] <= ref_energies[i - 1] {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Reference energies must be strictly ascending, but E[{}]={} <= E[{}]={}",
                    i,
                    ref_energies[i],
                    i - 1,
                    ref_energies[i - 1],
                )));
            }
        }
        for (i, (offsets, weights)) in kernels.iter().enumerate() {
            // Reject empty kernels: an empty `(offsets, weights)` passes the
            // length-match check (0 == 0) but later makes `broaden_presorted`
            // accumulate `norm == 0` and silently fall back to pass-through.
            if offsets.is_empty() {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} is empty; each kernel needs at least one (offset, weight) point"
                )));
            }
            if offsets.len() != weights.len() {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} offset length {} != weight length {}",
                    offsets.len(),
                    weights.len(),
                )));
            }
            // Reject non-finite synthesized kernels (e.g. a poisoned analytic
            // pulse) so a NaN table fails loudly here rather than silently
            // degrading the convolution to pass-through (a NaN norm bypasses the
            // division guard in `broaden_presorted`).
            if offsets.iter().chain(weights.iter()).any(|v| !v.is_finite()) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} contains non-finite offset/weight values"
                )));
            }
            // Negative weights have no physical meaning; the broadener
            // skips w <= 0 entries and the trapezoidal width moments
            // must integrate over the same domain (see `from_text`).
            if weights.iter().any(|&v| v < 0.0) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} contains negative weights"
                )));
            }
            // Offsets must be strictly ascending: `broaden_presorted`/`plan` walk a
            // monotonic two-pointer bracket and derive trapezoidal `dt` widths from
            // `offsets[k+1] − offsets[k−1]`, both of which assume sorted offsets.
            // An unsorted kernel would otherwise broaden silently-wrong, not error.
            if !offsets.windows(2).all(|w| w[0] < w[1]) {
                return Err(ResolutionParseError::InvalidFormat(format!(
                    "Kernel {i} TOF offsets must be strictly ascending"
                )));
            }
        }
        Ok(TabulatedResolution {
            ref_energies: Arc::new(ref_energies),
            kernels: Arc::new(kernels),
            flight_path_m,
        })
    }

    /// Parse a VENUS/FTS resolution file from disk.
    pub fn from_file(path: &str, flight_path_m: f64) -> Result<Self, ResolutionParseError> {
        let text = std::fs::read_to_string(path)
            .map_err(|e| ResolutionParseError::IoError(format!("Cannot read '{}': {}", path, e)))?;
        Self::from_text(&text, flight_path_m)
    }

    /// Apply tabulated resolution broadening to a spectrum.
    ///
    /// For each energy point:
    /// 1. Find bracketing reference energies and interpolate kernel (log-space)
    /// 2. Convert TOF offsets to energy offsets using exact TOF↔energy relation
    /// 3. Convolve spectrum with interpolated kernel (trapezoidal integration)
    ///
    /// Kernel points whose delayed-emission offset reaches the nominal
    /// flight time at the target energy (`dt ≥ TOF(E)`) gather from
    /// past infinite energy; they are dropped and the kernel
    /// renormalized over the surviving points, mirroring the grid-edge
    /// handling — see the tail-truncation note on `broaden_presorted`.
    ///
    /// # Errors
    /// Returns [`ResolutionError::LengthMismatch`] if the arrays differ in
    /// length, or [`ResolutionError::UnsortedEnergies`] if the energy grid is
    /// not sorted in non-descending order.
    pub fn broaden(&self, energies: &[f64], spectrum: &[f64]) -> Result<Vec<f64>, ResolutionError> {
        validate_inputs(energies, spectrum)?;
        Ok(self.broaden_presorted(energies, spectrum))
    }

    /// Tabulated resolution broadening assuming the energy grid is already
    /// validated (sorted ascending, same length as spectrum).
    ///
    /// ## Convolution orientation
    ///
    /// The broadened value at measured TOF `t` gathers theory at
    /// `t − dt`: a neutron *measured* at `t` whose emission was delayed
    /// by `dt` really flew for `t − dt`, i.e. it is faster than nominal,
    /// so the kernel's positive-offset (delayed-emission) tail pulls
    /// theory from earlier TOF = **higher** energy, and a resonance dip
    /// acquires its tail toward lower apparent energy. This is the
    /// convolution ∫R(τ)·S(t−τ)dτ, matching SAMMY's user-defined
    /// resolution: `sammy/src/udr/mudr4.f90` `Ud_Convolute` (line 288)
    /// computes `Cc(Tc) = ∫Bb(τ)·Aa(Tc−τ)dτ` (the theory-segment search
    /// binds `Tc−Ta` to the kernel grid), and `Ud_Mesh_Time` (line 7)
    /// sizes the needed theory window as `[T0−UdT_last, T0−UdT_first]`.
    ///
    /// ## Delayed-tail truncation at short nominal flight times
    ///
    /// A kernel point with `dt ≥ tof_center` would gather at
    /// `tof_prime = tof_center − dt ≤ 0`, i.e. from past infinite
    /// energy where no theory value exists.  Such points are silently
    /// dropped and the trapezoidal normalisation runs over the
    /// surviving points — the same truncate-and-renormalize treatment
    /// applied when `e_prime` falls outside the target grid
    /// (`[e_min, e_max]`).  This engages when the nominal TOF at the
    /// target energy is shorter than the kernel's delayed-emission
    /// reach — very high energy and/or a short flight path — and is
    /// reachable at *any* target energy because `interpolated_kernel`
    /// clamps to the nearest reference kernel outside the tabulated
    /// range.
    ///
    /// ## Inner-loop optimization
    ///
    /// The per-kernel-point spectrum interpolation uses a **two-pointer
    /// walk** instead of a binary search: `e_prime` is monotonically
    /// increasing in `k` (since `dt = offsets[k]` is non-decreasing,
    /// `TOF' = tof_center − dt` is non-increasing, and `E' = (L/TOF')²`
    /// is non-decreasing).  We maintain `bracket_hi` as the smallest
    /// index into `energies[]` whose value is `>= e_prime`; the walk is
    /// a bidirectional fixed-point search (first kernel point descends
    /// from `n−1`, subsequent points walk upward).  Amortized O(1) per
    /// kernel point within a target.
    ///
    /// Math is identical to the reference implementation pinned by
    /// `broaden_presorted_reference` in the test module.
    ///
    /// For callers that broaden many spectra on the same target grid —
    /// LM iterations with fixed TZERO, spatial maps with a pre-calibrated
    /// energy axis — [`TabulatedResolution::plan`] +
    /// [`ResolutionPlan::apply`] produce bit-exact output while
    /// hoisting the per-target invariants (TOF conversion, kernel
    /// interpolation, bracket lookup, trapezoidal widths) out of the
    /// broadening hot loop.  This `broaden_presorted` entry is the
    /// single-broadening path and keeps the original inline
    /// implementation to avoid plan-construction overhead on one-shot
    /// callers.
    pub(crate) fn broaden_presorted(&self, energies: &[f64], spectrum: &[f64]) -> Vec<f64> {
        let n = energies.len();
        if n == 0 {
            return vec![];
        }
        if n == 1 {
            return spectrum.to_vec();
        }

        let e_min = energies[0];
        let e_max = energies[n - 1];

        let mut result = vec![0.0f64; n];

        for i in 0..n {
            let e = energies[i];
            if e <= 0.0 {
                result[i] = spectrum[i];
                continue;
            }

            let tof_center = TOF_FACTOR * self.flight_path_m / e.sqrt();
            let (offsets, weights) = self.interpolated_kernel(e);
            let n_k = offsets.len();
            let mut bracket_hi: usize = n - 1;

            let mut sum = 0.0;
            let mut norm = 0.0;

            for k in 0..n_k {
                let dt = offsets[k];
                let w = weights[k];
                if w <= 0.0 {
                    continue;
                }

                // Convolution gather: theory at t − dt (see docstring;
                // SAMMY mudr4.f90 Ud_Convolute).
                let tof_prime = tof_center - dt;
                if tof_prime <= 0.0 {
                    continue;
                }

                let e_prime = (TOF_FACTOR * self.flight_path_m / tof_prime).powi(2);

                if e_prime < e_min || e_prime > e_max {
                    continue;
                }

                while bracket_hi > 1 && energies[bracket_hi - 1] > e_prime {
                    bracket_hi -= 1;
                }
                while bracket_hi < n - 1 && energies[bracket_hi] <= e_prime {
                    bracket_hi += 1;
                }

                let lo = bracket_hi - 1;
                let hi = bracket_hi;
                let span = energies[hi] - energies[lo];
                let s = if span.abs() < NEAR_ZERO_FLOOR {
                    spectrum[lo]
                } else {
                    let frac = (e_prime - energies[lo]) / span;
                    spectrum[lo] + frac * (spectrum[hi] - spectrum[lo])
                };

                let dt_width = if k > 0 && k < n_k - 1 {
                    (offsets[k + 1] - offsets[k - 1]) * 0.5
                } else if k == 0 && n_k > 1 {
                    offsets[1] - offsets[0]
                } else if k == n_k - 1 && n_k > 1 {
                    offsets[k] - offsets[k - 1]
                } else {
                    1.0
                };

                let weight = w * dt_width.abs();
                sum += weight * s;
                norm += weight;
            }

            result[i] = if norm > DIVISION_FLOOR {
                sum / norm
            } else {
                spectrum[i]
            };
        }

        result
    }

    /// Build a reusable broadening plan for a specific target energy grid.
    ///
    /// Validates that `energies` is non-descending — the same sorted-grid
    /// precondition enforced by [`TabulatedResolution::broaden`] via
    /// `validate_inputs`.  An
    /// unsorted grid would produce a silently-wrong plan (misbracketed
    /// `e_prime` lookups against `e_min` / `e_max`), so it must be
    /// caught at build time rather than returning garbage from
    /// [`ResolutionPlan::apply`].
    ///
    /// The plan hoists every quantity that depends only on
    /// `(target_energies, self.ref_energies, self.flight_path_m)` —
    /// namely the TOF conversion, the log-space kernel interpolation,
    /// the per-kernel-point `e_prime` and spectrum-bracket lookup, and
    /// the trapezoidal integration widths.  Applying the plan to a
    /// spectrum becomes a pure gather + multiply-add loop.
    ///
    /// Build cost: same as one call to the private `broaden_presorted`
    /// helper (O(N_target × N_kernel) TOF / bracket / interp work, plus
    /// ~2 × N_kernel log-interp ops per target energy for
    /// `interpolated_kernel`).  Apply cost per target: 1 branch +
    /// ~3 loads + 3 flops per retained entry, plus the final divide —
    /// typically < 10 % of the build cost.  The payoff comes from
    /// reusing one plan across many spectra.
    ///
    /// Bit-exact with `broaden_presorted`: pre-computes the same
    /// floating-point sequences (TOF, `e_prime`, `dt_width`, `frac`,
    /// `weight`, `norm`) in the same order.
    ///
    /// # Errors
    /// Returns [`ResolutionError::UnsortedEnergies`] if `energies` is
    /// not non-descending.
    pub fn plan(&self, energies: &[f64]) -> Result<ResolutionPlan, ResolutionError> {
        if !energies.windows(2).all(|w| w[0] <= w[1]) {
            return Err(ResolutionError::UnsortedEnergies);
        }
        Ok(self.plan_presorted(energies))
    }

    /// Build a plan assuming `energies` is already validated as
    /// non-descending.  Used internally by `broaden_presorted` (whose
    /// caller already validated the grid) and by `plan()` after its
    /// validation succeeded.
    fn plan_presorted(&self, energies: &[f64]) -> ResolutionPlan {
        let n = energies.len();
        if n == 0 {
            return ResolutionPlan {
                target_energies: Vec::new(),
                starts: vec![0],
                lo_idx: Vec::new(),
                frac: Vec::new(),
                weight: Vec::new(),
                norm: Vec::new(),
            };
        }
        if n == 1 {
            // No bracket available; passthrough. Represent as n=1 with
            // zero entries and norm=0, which triggers the passthrough
            // branch in `ResolutionPlan::apply`.
            return ResolutionPlan {
                target_energies: energies.to_vec(),
                starts: vec![0, 0],
                lo_idx: Vec::new(),
                frac: Vec::new(),
                weight: Vec::new(),
                norm: vec![0.0],
            };
        }

        let e_min = energies[0];
        let e_max = energies[n - 1];

        // Preallocate the entry Vecs to ~n × 2·kernel_len: the
        // width-normalized shape blend merges the two bracketing
        // blocks, so between-reference targets emit up to
        // n_lo + n_hi points (~2× a single block; real VENUS grids
        // push ~n × 998 entries). Over-allocating for at-reference
        // targets is cheap vs. repeated grow-and-memcpy during the
        // plan build.
        let estimated_kernel_len = self.kernels.first().map_or(0, |(off, _)| off.len());
        let estimated_entries = n.saturating_mul(estimated_kernel_len.saturating_mul(2));

        let mut starts: Vec<u32> = Vec::with_capacity(n + 1);
        let mut lo_idx: Vec<u32> = Vec::with_capacity(estimated_entries);
        let mut frac: Vec<f64> = Vec::with_capacity(estimated_entries);
        let mut weight: Vec<f64> = Vec::with_capacity(estimated_entries);
        let mut norm: Vec<f64> = Vec::with_capacity(n);

        starts.push(0);

        for i in 0..n {
            let e = energies[i];
            if e <= 0.0 {
                // Passthrough: no entries contribute, norm=0.
                norm.push(0.0);
                // Guard the u32 invariant for diagnostic callers; the
                // headroom is enormous for any realistic grid (VENUS
                // 3471 × 499 ≈ 1.7M entries, u32::MAX ≈ 4.29B), but
                // the debug-only assert documents the contract.
                debug_assert!(
                    lo_idx.len() <= u32::MAX as usize,
                    "plan entry count overflows u32"
                );
                starts.push(lo_idx.len() as u32);
                continue;
            }

            // TOF at this energy: t = TOF_FACTOR * L / sqrt(E).
            // Computed here in the plan build and NOT at apply time — this
            // is the main invariant we hoist.
            let tof_center = TOF_FACTOR * self.flight_path_m / e.sqrt();

            // Interpolated kernel at this target energy.  Allocates two
            // ~N_kernel Vecs; those allocations happen once per plan
            // build instead of once per broadening call.
            let (offsets, weights) = self.interpolated_kernel(e);
            let n_k = offsets.len();

            // Two-pointer walk state (same invariant as broaden_presorted).
            let mut bracket_hi: usize = n - 1;

            let mut target_norm = 0.0;

            for k in 0..n_k {
                let dt = offsets[k];
                let w = weights[k];
                if w <= 0.0 {
                    continue;
                }

                // Convolution gather: theory at t − dt (see
                // broaden_presorted; SAMMY mudr4.f90 Ud_Convolute).
                // Points with dt ≥ tof_center gather from past infinite
                // energy and are dropped, renormalizing over the
                // survivors (see the tail-truncation note there).
                let tof_prime = tof_center - dt;
                if tof_prime <= 0.0 {
                    continue;
                }

                let e_prime = (TOF_FACTOR * self.flight_path_m / tof_prime).powi(2);

                if e_prime < e_min || e_prime > e_max {
                    continue;
                }

                // Two-pointer walk — same logic + invariants as
                // broaden_presorted, in the same order, so bracket_hi
                // reaches the identical position for each kept (i, k).
                while bracket_hi > 1 && energies[bracket_hi - 1] > e_prime {
                    bracket_hi -= 1;
                }
                while bracket_hi < n - 1 && energies[bracket_hi] <= e_prime {
                    bracket_hi += 1;
                }

                let lo = bracket_hi - 1;
                let hi = bracket_hi;
                let span = energies[hi] - energies[lo];
                // Degenerate-bracket guard: if span < NEAR_ZERO_FLOOR,
                // broaden_presorted returns `spectrum[lo]` directly
                // without the interp arithmetic.  Store `frac = -0.0`
                // — the apply path short-circuits on the exact bit
                // pattern of `-0.0` and returns `spectrum[lo]` without
                // touching `spectrum[lo+1]`, so bit-exactness holds
                // even if `spectrum[lo+1]` is NaN or ±∞.
                //
                // `-0.0` (negative-signed zero) is used as the sentinel
                // because the non-degenerate path can legitimately
                // produce `frac == +0.0` when `e_prime == energies[lo]`
                // exactly — in that case `broaden_presorted` still
                // reads `spectrum[lo+1]` (and propagates NaN if present
                // there), so the apply path MUST do the same.  `+0.0`
                // and `-0.0` compare equal under `==` but differ in
                // `to_bits()`, which is what apply uses to disambiguate.
                let entry_frac = if span.abs() < NEAR_ZERO_FLOOR {
                    -0.0_f64
                } else {
                    (e_prime - energies[lo]) / span
                };

                let dt_width = if k > 0 && k < n_k - 1 {
                    (offsets[k + 1] - offsets[k - 1]) * 0.5
                } else if k == 0 && n_k > 1 {
                    offsets[1] - offsets[0]
                } else if k == n_k - 1 && n_k > 1 {
                    offsets[k] - offsets[k - 1]
                } else {
                    1.0
                };

                let entry_weight = w * dt_width.abs();

                debug_assert!(
                    lo_idx.len() < u32::MAX as usize,
                    "plan entry count overflows u32"
                );
                lo_idx.push(lo as u32);
                frac.push(entry_frac);
                weight.push(entry_weight);
                target_norm += entry_weight;
            }

            norm.push(target_norm);
            starts.push(lo_idx.len() as u32);
        }

        ResolutionPlan {
            target_energies: energies.to_vec(),
            starts,
            lo_idx,
            frac,
            weight,
            norm,
        }
    }

    /// Interpolate the kernel at an arbitrary energy as a
    /// **width-normalized shape blend** between the two bracketing
    /// reference kernels:
    ///
    /// 1. Exact hits (a reference energy, or outside the reference
    ///    range) return that reference kernel unchanged.
    /// 2. Each bracketing block's trapezoidal RMS width `σ_b` is
    ///    computed ([`trapezoidal_moments`]); the target width is the
    ///    **geometric** interpolation `σ_t = σ_lo·(σ_hi/σ_lo)^frac`
    ///    with `frac` linear in log E — exact for the physical
    ///    power-law width `σ_t ∝ E^p` (log σ linear in log E).
    /// 3. Both blocks' offsets are scaled about the mode (offset 0 —
    ///    the anchoring convention; see the discussion in
    ///    `ikeda_carpenter`) by `σ_t/σ_b`, merged into one
    ///    strictly-ascending grid, and the weights blended pointwise:
    ///    `w = w_lo(x) + frac·(w_hi(x) − w_lo(x))`, each block's weight
    ///    linearly interpolated (zero outside its support). Identical
    ///    blocks reduce to a **bitwise identity** (the scale ratios are
    ///    exactly 1.0 and the blend form is exact when `w_lo == w_hi`).
    ///
    /// Degenerate blocks (single-point, zero mass → `σ_b ≤ 0`) fall
    /// back to the nearer reference clone.
    ///
    /// ## INTENTIONAL DEPARTURE from SAMMY
    ///
    /// SAMMY's user-defined resolution blends both the amplitude and
    /// the time-point arrays element-wise, **linear in E**: in the
    /// active `Gen_Udr_Par` (sammy/src/udr/mudr3.f90, subroutine at
    /// line 164; blend block at lines 241–255: `UdR_E(J,Nud) =
    /// UdR(J,I−1,Nud)·a + UdR(J,I,Nud)·b` and likewise `UdT_E`), i.e.
    /// the arithmetic width chord. (The file's first routine
    /// `Gen_Udr_Par_x` holds the same blend at lines 92–112 but is
    /// marked "never called" at line 10 — cite the live twin.)
    /// Because the physical width law
    /// `σ_t ∝ ~E^{−1/2}` is convex, that chord systematically
    /// over-widens every between-reference energy: +7.8 % at the
    /// midpoint of synthetic 10/50 eV Gaussian blocks, +4.1…+7.2 %
    /// across the production VENUS 5→50 eV reference gap — a direct
    /// resolution-width systematic that biases fitted temperatures
    /// low. The geometric-width shape blend above removes it (and the
    /// nearest-reference width sawtooth that unequal point counts used
    /// to produce). SAMMY additionally re-aligns the blended kernel so
    /// its trapezoidal centroid `Ct` sits at T = 0 ("Realign so that
    /// centroid is at T=0", mudr3.f90 lines 266–292);
    /// NEREIDS does not re-align: the blend scales offsets about 0, so each
    /// table keeps its own origin — a loaded UDR file its peak, a synthesized
    /// Ikeda–Carpenter table the emission instant.
    ///
    /// Exactness caveat: the blend reproduces `σ_t` exactly when the
    /// two blocks' width-normalized shapes agree (self-similar
    /// families, e.g. any pure power-law file). Genuinely different
    /// bracketing shapes add a second-order mixture-spread term —
    /// inherent to shape blending and far below the removed chord
    /// error for real moderator files.
    ///
    /// Allocates the two output Vecs per call (≈ n_lo + n_hi points
    /// between references); scratch reuse is tracked separately as a
    /// performance follow-up.
    ///
    /// `ref_energies` is validated as strictly ascending by `from_text()` /
    /// `from_file()` at construction time, so no per-call sort check is needed.
    fn interpolated_kernel(&self, energy: f64) -> (Vec<f64>, Vec<f64>) {
        debug_assert!(
            self.ref_energies.windows(2).all(|w| w[0] < w[1]),
            "ref_energies must be strictly ascending (invariant broken)"
        );
        let n_ref = self.ref_energies.len();

        // NaN target energies never reach here through the validated
        // public paths (a NaN in a multi-point grid fails the sorted
        // check), but every comparison below is false for NaN, and the
        // width-scaled merge would then emit NaN offsets — violating
        // the strictly-ascending, all-finite invariants the broadener
        // assumes. Clamp to the lowest reference so the output kernel
        // is well-formed unconditionally (the old element-wise blend
        // degraded to its nearest-reference fallback here by accident
        // of its monotonicity guard; this keeps that graceful
        // behaviour explicit). +∞ needs no guard: the high clamp
        // below already catches it.
        if energy.is_nan() {
            return self.kernels[0].clone();
        }

        // Clamp to nearest reference if outside range
        if energy <= self.ref_energies[0] || n_ref == 1 {
            return self.kernels[0].clone();
        }
        if energy >= self.ref_energies[n_ref - 1] {
            return self.kernels[n_ref - 1].clone();
        }

        // Find bracketing indices
        let pos = self.ref_energies.partition_point(|&e| e < energy);
        // Interior exact hit: return that reference unchanged rather than
        // blending it with itself, which reproduces it only up to ULPs.
        if self.ref_energies[pos] == energy {
            return self.kernels[pos].clone();
        }
        let idx = if pos == 0 {
            0
        } else {
            (pos - 1).min(n_ref - 2)
        };

        let e_lo = self.ref_energies[idx];
        let e_hi = self.ref_energies[idx + 1];

        // Log-space interpolation fraction
        let frac = (energy.ln() - e_lo.ln()) / (e_hi.ln() - e_lo.ln());

        let (off_lo, w_lo) = &self.kernels[idx];
        let (off_hi, w_hi) = &self.kernels[idx + 1];

        let nearest = || -> (Vec<f64>, Vec<f64>) {
            let k = if frac < 0.5 {
                &self.kernels[idx]
            } else {
                &self.kernels[idx + 1]
            };
            (k.0.clone(), k.1.clone())
        };

        let (_, s_lo) = trapezoidal_moments(off_lo, w_lo);
        let (_, s_hi) = trapezoidal_moments(off_hi, w_hi);
        // Degenerate blocks (σ ≤ 0) cannot be width-scaled, and a
        // non-finite fraction (constructors enforce positive reference
        // energies, so defense-in-depth only) would poison every
        // blended weight with NaN — both take the nearest-clone
        // fallback so the output is well-formed unconditionally.
        if !(s_lo.is_finite() && s_lo > 0.0 && s_hi.is_finite() && s_hi > 0.0 && frac.is_finite()) {
            return nearest();
        }

        // Geometric width interpolation. This float form (ratio +
        // powf) is a bitwise no-op when σ_lo == σ_hi: the ratio is
        // exactly 1.0, powf(1.0, f) == 1.0, so both scale factors are
        // exactly 1.0 and scaled offsets are the originals.
        let s_t = s_lo * (s_hi / s_lo).powf(frac);
        let r_lo = s_t / s_lo;
        let r_hi = s_t / s_hi;

        // Single-pass sorted merge of the two scaled grids. Each
        // block's weight at a merged point is its own tabulated value
        // when the point came from that block, else the linear
        // interpolation of its shape (zero outside its support).
        let n_lo = off_lo.len();
        let n_hi = off_hi.len();
        let mut out_off: Vec<f64> = Vec::with_capacity(n_lo + n_hi);
        let mut out_w: Vec<f64> = Vec::with_capacity(n_lo + n_hi);

        // Piecewise-linear sample of one block's shape at `x`, with a
        // monotone cursor (merged points arrive in ascending order).
        let sample = |offs: &[f64], ws: &[f64], r: f64, cursor: &mut usize, x: f64| -> f64 {
            let n = offs.len();
            if x < offs[0] * r || x > offs[n - 1] * r {
                return 0.0;
            }
            while *cursor + 1 < n && offs[*cursor + 1] * r <= x {
                *cursor += 1;
            }
            if *cursor + 1 >= n {
                return ws[n - 1];
            }
            let x0 = offs[*cursor] * r;
            let x1 = offs[*cursor + 1] * r;
            let span = x1 - x0;
            if x <= x0 || span <= 0.0 {
                return ws[*cursor];
            }
            ws[*cursor] + (x - x0) / span * (ws[*cursor + 1] - ws[*cursor])
        };

        let (mut i, mut j) = (0usize, 0usize);
        let (mut ci, mut cj) = (0usize, 0usize);
        while i < n_lo || j < n_hi {
            let xa = if i < n_lo {
                off_lo[i] * r_lo
            } else {
                f64::INFINITY
            };
            let xb = if j < n_hi {
                off_hi[j] * r_hi
            } else {
                f64::INFINITY
            };
            // Near-duplicate merge (covers the exact 0 == 0 mode point):
            // emit one point carrying both blocks' exact tabulated
            // weights, so identical blocks blend to their exact values.
            let near_dup = i < n_lo
                && j < n_hi
                && (xa - xb).abs() <= 4.0 * f64::EPSILON * xa.abs().max(xb.abs());
            let (x, wl, wh) = if near_dup {
                let v = (xa, w_lo[i], w_hi[j]);
                i += 1;
                j += 1;
                v
            } else if xa < xb {
                let v = (xa, w_lo[i], sample(off_hi, w_hi, r_hi, &mut cj, xa));
                i += 1;
                v
            } else {
                let v = (xb, sample(off_lo, w_lo, r_lo, &mut ci, xb), w_hi[j]);
                j += 1;
                v
            };
            // Defensive strict-ascension guard: the broadener's
            // trapezoidal quadrature and two-pointer walk require it
            // unconditionally, regardless of the dedup epsilon.
            if let Some(&last) = out_off.last()
                && x <= last
            {
                continue;
            }
            out_off.push(x);
            // Exact when `wl == wh` (identical blocks), and exactly the
            // endpoint values at frac → 0/1.
            out_w.push(wl + frac * (wh - wl));
        }

        // Blended weights inherit the blocks' scale (peak-normalized by
        // convention); the broadener renormalizes at apply time, so no
        // re-normalization is done here — preserving the bitwise
        // identity for identical blocks unconditionally.
        (out_off, out_w)
    }
}

/// Apply resolution broadening using either Gaussian or tabulated kernel.
///
/// # Errors
/// Returns [`ResolutionError`] if the energy grid is unsorted or array
/// lengths do not match.
pub fn apply_resolution(
    energies: &[f64],
    spectrum: &[f64],
    resolution: &ResolutionFunction,
) -> Result<Vec<f64>, ResolutionError> {
    match resolution {
        ResolutionFunction::Gaussian(params) => resolution_broaden(energies, spectrum, params),
        ResolutionFunction::Tabulated(tab) => tab.broaden(energies, spectrum),
        ResolutionFunction::IkedaCarpenter(ic) => ic.tabulated().broaden(energies, spectrum),
    }
}

/// Apply resolution broadening assuming the energy grid is already validated
/// (sorted ascending, same length as spectrum).
///
/// Used by `transmission.rs` to avoid redundant O(N) sort checks when
/// broadening multiple isotopes on the same pre-validated energy grid.
pub(crate) fn apply_resolution_presorted(
    energies: &[f64],
    spectrum: &[f64],
    resolution: &ResolutionFunction,
) -> Vec<f64> {
    match resolution {
        ResolutionFunction::Gaussian(params) => {
            resolution_broaden_presorted(energies, spectrum, params)
        }
        ResolutionFunction::Tabulated(tab) => tab.broaden_presorted(energies, spectrum),
        ResolutionFunction::IkedaCarpenter(ic) => {
            ic.tabulated().broaden_presorted(energies, spectrum)
        }
    }
}

/// Build a broadening plan for `(energies, resolution)`.
///
/// Returns `Some(plan)` for [`ResolutionFunction::Tabulated`] and
/// [`ResolutionFunction::IkedaCarpenter`] (which rides its synthesized tabulated
/// kernel) — the plan hoists the per-target TOF / kernel-interpolation / bracket
/// / trap-weight work that would otherwise run on every call to
/// [`apply_resolution`].  Returns `None` for
/// [`ResolutionFunction::Gaussian`] — the Gaussian path has no
/// meaningful pixel-invariant kernel structure to cache at this
/// level, so callers fall back to the per-call broadening path with
/// no loss.
///
/// Callers that want a single-branch API can unconditionally call
/// [`apply_resolution_with_plan`] passing `plan.as_ref()`; when the
/// plan is `None` it transparently forwards to the non-plan path and
/// returns byte-identical output.
///
/// # Errors
/// Returns [`ResolutionError::UnsortedEnergies`] if `energies` is not
/// non-descending — the same precondition that [`apply_resolution`]
/// enforces per-call.
pub fn build_resolution_plan(
    energies: &[f64],
    resolution: &ResolutionFunction,
) -> Result<Option<ResolutionPlan>, ResolutionError> {
    match resolution {
        ResolutionFunction::Gaussian(_) => {
            if !energies.windows(2).all(|w| w[0] <= w[1]) {
                return Err(ResolutionError::UnsortedEnergies);
            }
            Ok(None)
        }
        ResolutionFunction::Tabulated(tab) => tab.plan(energies).map(Some),
        ResolutionFunction::IkedaCarpenter(ic) => ic.tabulated().plan(energies).map(Some),
    }
}

/// Apply resolution broadening, optionally via a pre-built
/// [`ResolutionPlan`].
///
/// When `plan` is `Some(p)` and `resolution` is a tabulated kernel,
/// `p.apply(spectrum)` runs the cached per-target broadening inner
/// loop — the expensive TOF / kernel-interpolation / bracket work
/// was already captured at plan build time.
///
/// When `plan` is `None`, or when `resolution` is Gaussian, the call
/// forwards to [`apply_resolution`] and is byte-identical to the
/// un-planned path.
///
/// # Errors
/// * Returns the same errors as [`apply_resolution`] on the non-plan
///   path.
/// * Returns [`ResolutionError::LengthMismatch`] if the plan was built
///   for a different-length grid than `energies`, or if
///   `energies.len() != spectrum.len()`.
/// * Returns [`ResolutionError::PlanGridMismatch`] if the plan was
///   built for a different grid of the same length — the cached
///   `(lo_idx, frac, weight)` entries encode brackets into the old
///   grid and would silently produce a wrong broadened spectrum if
///   applied.
pub fn apply_resolution_with_plan(
    plan: Option<&ResolutionPlan>,
    energies: &[f64],
    spectrum: &[f64],
    resolution: &ResolutionFunction,
) -> Result<Vec<f64>, ResolutionError> {
    if let Some(p) = plan
        && matches!(
            resolution,
            ResolutionFunction::Tabulated(_) | ResolutionFunction::IkedaCarpenter(_)
        )
    {
        validate_inputs(energies, spectrum)?;
        if p.len() != energies.len() {
            return Err(ResolutionError::LengthMismatch {
                energies: energies.len(),
                data: p.len(),
            });
        }
        // Grid-identity check.  A plan built for a different grid of
        // the same length would still pass the length check and then
        // gather spectrum values at brackets that belong to the old
        // grid — silently corrupt output.  Pointer identity is not
        // enough here because callers legitimately hold the plan and
        // the target grid in separate `Arc`s whose storage may or
        // may not alias; bit-exact content equality is the only
        // robust invariant.  The cost is one full grid scan per
        // broadening call (O(n), ~27 KB of f64 values for the VENUS
        // 3471-point grid) — orders of magnitude cheaper than the
        // broadening itself and cheap vs the silent-staleness
        // failure mode.
        let plan_grid = p.target_energies();
        for i in 0..plan_grid.len() {
            if plan_grid[i].to_bits() != energies[i].to_bits() {
                return Err(ResolutionError::PlanGridMismatch {
                    first_diff_index: i,
                });
            }
        }
        return Ok(p.apply(spectrum));
    }
    apply_resolution(energies, spectrum, resolution)
}

/// Errors from resolution file parsing.
#[derive(Debug)]
pub enum ResolutionParseError {
    InvalidFormat(String),
    IoError(String),
}

impl fmt::Display for ResolutionParseError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidFormat(msg) => write!(f, "Invalid resolution file format: {}", msg),
            Self::IoError(msg) => write!(f, "I/O error: {}", msg),
        }
    }
}

impl std::error::Error for ResolutionParseError {}

/// Test-only helpers for building synthetic [`ResolutionPlan`] /
/// [`TabulatedResolution`] instances without going through the full
/// parse-from-text path.  Gated behind `#[cfg(test)]` for in-crate
/// tests and the `test-support` feature flag for downstream-crate
/// tests.  Never ships in a release build with `test-support =
/// false` (the default), so the raw constructors remain out of the
/// production API surface.
#[cfg(any(test, feature = "test-support"))]
pub mod test_support {
    use super::{DIVISION_FLOOR, NEAR_ZERO_FLOOR, ResolutionPlan, TabulatedResolution};
    use std::sync::Arc;

    /// Build a [`ResolutionPlan`] directly from its SoA fields.
    ///
    /// The caller is responsible for maintaining the invariants that
    /// `plan_presorted` normally enforces (`starts.last() ==
    /// lo_idx.len()`, lo_idx in [0, n-2] for regular entries, etc.).
    /// Used by surrogate-module tests + downstream crate tests to
    /// construct hand-designed plans that exercise specific CSR
    /// patterns.
    pub fn plan_from_raw_parts(
        target_energies: Vec<f64>,
        starts: Vec<u32>,
        lo_idx: Vec<u32>,
        frac: Vec<f64>,
        weight: Vec<f64>,
        norm: Vec<f64>,
    ) -> ResolutionPlan {
        ResolutionPlan {
            target_energies,
            starts,
            lo_idx,
            frac,
            weight,
            norm,
        }
    }

    /// Build a minimal [`TabulatedResolution`] with a single
    /// reference energy and a trivial delta-like kernel — just
    /// enough for tests that need an `InstrumentParams` with a
    /// tabulated resolution (e.g., to exercise cubature dispatch
    /// guards that refuse Gaussian resolution).  The broadening
    /// would be effectively identity if anyone ever called it, but
    /// typical consumers (cubature dispatch tests) never invoke the
    /// kernel.
    pub fn trivial_tabulated_resolution(flight_path_m: f64) -> TabulatedResolution {
        TabulatedResolution {
            ref_energies: Arc::new(vec![100.0]),
            kernels: Arc::new(vec![(vec![-1e-6, 0.0, 1e-6], vec![0.0, 1.0, 0.0])]),
            flight_path_m,
        }
    }

    /// Thin shim exposing the crate-internal
    /// [`TabulatedResolution::broaden_presorted`] to integration
    /// tests living under `crates/nereids-physics/tests/`.  The
    /// internal method stays `pub(crate)` so the broader public
    /// API surface (the operator-style `apply_resolution`,
    /// `plan` / `apply` / `compile_to_matrix`) remains the
    /// recommended entry point; this shim exists solely so the
    /// fixture-gated bit-exact regression and microbenchmark
    /// tests can call the optimized two-pointer walk directly.
    pub fn broaden_presorted(
        tab: &TabulatedResolution,
        energies: &[f64],
        spectrum: &[f64],
    ) -> Vec<f64> {
        tab.broaden_presorted(energies, spectrum)
    }

    /// Thin shim exposing the crate-internal
    /// `TabulatedResolution::interpolated_kernel` to integration
    /// tests.  Needed by the bit-exact equivalence oracle that
    /// the fixture-gated regression test runs against the
    /// optimized `broaden_presorted` path.
    pub fn interpolated_kernel(tab: &TabulatedResolution, energy: f64) -> (Vec<f64>, Vec<f64>) {
        tab.interpolated_kernel(energy)
    }

    /// The TOF↔energy conversion factor used by
    /// `broaden_presorted` and its oracle.  Exposed so the
    /// integration-test oracle uses the exact same constant as
    /// the SUT, preserving bit-exact equivalence.
    pub const TOF_FACTOR: f64 = super::TOF_FACTOR;

    /// Bit-exact regression oracle: binary-search piecewise-linear
    /// interpolation of `spectrum` onto target energy `e`.  Mirror of
    /// the pre-optimization in-src reference; called transitively by
    /// [`broaden_presorted_reference`] inside its inner convolution
    /// loop.
    ///
    /// **Do not "clean up" this function.** Consumers are bit-exact
    /// equivalence tests that pin the optimized two-pointer path in
    /// `TabulatedResolution::broaden_presorted` against this byte-
    /// identical reference.  A rewrite that shifts edge cases by even
    /// one bit (e.g. swapping the upper-bound binary search for
    /// `partition_point`, or changing the `<=` to `<` in the midpoint
    /// comparison) would flip the comparison and invalidate the
    /// regression suite.
    pub fn interp_spectrum(energies: &[f64], spectrum: &[f64], e: f64) -> Option<f64> {
        let n = energies.len();
        if n == 0 {
            return None;
        }
        if e < energies[0] || e > energies[n - 1] {
            return None;
        }
        let mut lo = 0;
        let mut hi = n - 1;
        while hi - lo > 1 {
            let mid = (lo + hi) / 2;
            if energies[mid] <= e {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        let span = energies[hi] - energies[lo];
        if span.abs() < NEAR_ZERO_FLOOR {
            return Some(spectrum[lo]);
        }
        let frac = (e - energies[lo]) / span;
        Some(spectrum[lo] + frac * (spectrum[hi] - spectrum[lo]))
    }

    /// Bit-exact regression oracle: pre-optimization reference
    /// implementation of [`TabulatedResolution::broaden_presorted`].
    /// Used by the in-src + integration + microbench bit-exact test
    /// suites that pin the optimized two-pointer path against this
    /// reference.
    ///
    /// Same "do not refactor" caveat as [`interp_spectrum`].  Reads
    /// `tab.flight_path_m()` via the public getter so the oracle stays
    /// callable from integration tests (the underlying field is
    /// private; the getter is a no-op wrapper, so byte-equivalence
    /// against the original in-src field access is preserved).
    pub fn broaden_presorted_reference(
        tab: &TabulatedResolution,
        energies: &[f64],
        spectrum: &[f64],
    ) -> Vec<f64> {
        let n = energies.len();
        if n == 0 {
            return vec![];
        }

        let mut result = vec![0.0f64; n];

        for i in 0..n {
            let e = energies[i];
            if e <= 0.0 {
                result[i] = spectrum[i];
                continue;
            }

            let tof_center = TOF_FACTOR * tab.flight_path_m() / e.sqrt();
            let (offsets, weights) = tab.interpolated_kernel(e);

            let mut sum = 0.0;
            let mut norm = 0.0;

            for k in 0..offsets.len() {
                let dt = offsets[k];
                let w = weights[k];
                if w <= 0.0 {
                    continue;
                }

                // Convolution gather: theory at t − dt (see
                // broaden_presorted; SAMMY mudr4.f90 Ud_Convolute).
                // Points with dt ≥ tof_center gather from past infinite
                // energy and are dropped, renormalizing over the
                // survivors (see the tail-truncation note there).
                let tof_prime = tof_center - dt;
                if tof_prime <= 0.0 {
                    continue;
                }

                let e_prime = (TOF_FACTOR * tab.flight_path_m() / tof_prime).powi(2);

                let s = match interp_spectrum(energies, spectrum, e_prime) {
                    Some(v) => v,
                    None => continue,
                };

                let dt_width = if k > 0 && k < offsets.len() - 1 {
                    (offsets[k + 1] - offsets[k - 1]) * 0.5
                } else if k == 0 && offsets.len() > 1 {
                    offsets[1] - offsets[0]
                } else if k == offsets.len() - 1 && offsets.len() > 1 {
                    offsets[k] - offsets[k - 1]
                } else {
                    1.0
                };

                let weight = w * dt_width.abs();
                sum += weight * s;
                norm += weight;
            }

            result[i] = if norm > DIVISION_FLOOR {
                sum / norm
            } else {
                spectrum[i]
            };
        }

        result
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use nereids_core::constants;

    fn kernel_centroid_std(offs: &[f64], wts: &[f64]) -> (f64, f64) {
        let wsum: f64 = wts.iter().sum();
        let c = offs.iter().zip(wts).map(|(o, w)| o * w).sum::<f64>() / wsum;
        let var = offs
            .iter()
            .zip(wts)
            .map(|(o, w)| w * (o - c).powi(2))
            .sum::<f64>()
            / wsum;
        (c, var.sqrt())
    }

    const SQRT_PI: f64 = 1.772_453_850_905_516;

    /// Largest relative difference from the independent erfc of `libm`.
    fn worst(points: impl Iterator<Item = (f64, f64, f64)>) -> (f64, f64) {
        points
            .map(|(x, ours, want)| (x, (ours - want).abs() / want.abs()))
            .fold((f64::NAN, 0.0), |w, p| if p.1 > w.1 { p } else { w })
    }

    #[test]
    fn exerfc_and_erfc_match_an_independent_erfc_on_both_sides_of_zero() {
        let xs = (-2600..=2600).map(|i| f64::from(i) * 0.01);
        let (x, relative) = worst(xs.map(|x| {
            let want = SQRT_PI * (x * x).exp() * libm::erfc(x);
            (x, exerfc(x), want)
        }));
        assert!(relative < 2e-7, "exerfc({x}) is {relative:.2e} off");

        let xs = (-2600..=2600).map(|i| f64::from(i) * 0.01);
        let (x, relative) = worst(xs.map(|x| (x, erfc(x), libm::erfc(x))));
        assert!(relative < 2e-7, "erfc({x}) is {relative:.2e} off");
    }

    /// Z(A, B) = √π·exp(B² − A²)·erfc(B) with B = A + C, so
    /// B² − A² = 2BC − C². Below B = −5.01 the kernel is the exponential
    /// tail 2√π·exp(2BC − C²).
    #[test]
    fn the_combined_kernel_matches_an_independent_erfc_down_its_tail() {
        for c in [0.05, 0.1, 0.25, 1.0] {
            let bs = (-6000..=2600).map(|i| f64::from(i) * 0.01);
            let (b, relative) = worst(bs.map(|b| {
                let want = SQRT_PI * (2.0 * b * c - c * c).exp() * libm::erfc(b);
                (b, gauss_exp_kernel(b - c, b), want)
            }));
            assert!(
                relative < 2e-7,
                "C = {c}: Z at B = {b} is {relative:.2e} off"
            );
        }
    }

    /// ∫ Z dA = √π/C: integrate by parts with u = A + C,
    /// ∫ e^{2Cu} erfc(u) du = e^{C²}/C.
    #[test]
    fn the_combined_kernel_keeps_its_whole_area() {
        for c in [0.05, 0.1, 0.25] {
            let step = 0.01;
            let low = -50.0 / c;
            let count = ((10.0 - low) / step) as usize;
            let area: f64 = (0..=count)
                .map(|i| {
                    let a = low + step * i as f64;
                    let end = if i == 0 || i == count { 0.5 } else { 1.0 };
                    end * gauss_exp_kernel(a, a + c)
                })
                .sum::<f64>()
                * step;
            let want = SQRT_PI / c;
            let relative = (area - want).abs() / want;
            assert!(
                relative < 1e-9,
                "C = {c}: area {area:.9} against √π/C = {want:.9}, {relative:.2e} off"
            );
        }
    }

    #[test]
    fn width_corrected_preserves_centroid_scales_width_and_energy_dependence() {
        // asymmetric kernel straddling 0 (peak off-centre), two ref energies.
        let offs = vec![-2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0];
        let wts = vec![0.1, 0.3, 1.0, 0.8, 0.5, 0.3, 0.1];
        let tab = TabulatedResolution::from_kernels(
            vec![5.0, 50.0],
            vec![(offs.clone(), wts.clone()), (offs.clone(), wts.clone())],
            25.0,
        )
        .unwrap();

        // Uniform 2.5x width scale (p=0): centroid fixed, std scales by s0, weights unchanged.
        let s0 = 2.5;
        let wc = tab.width_corrected(s0, 0.0, 10.0).unwrap();
        for (orig, scaled) in tab.kernels().iter().zip(wc.kernels()) {
            let (c0, std0) = kernel_centroid_std(&orig.0, &orig.1);
            let (c1, std1) = kernel_centroid_std(&scaled.0, &scaled.1);
            assert!((c0 - c1).abs() < 1e-12, "centroid moved {c0} -> {c1}");
            assert!(
                (std1 / std0 - s0).abs() < 1e-12,
                "width ratio {} != {s0}",
                std1 / std0
            );
        }
        assert_eq!(
            tab.kernels()[0].1,
            wc.kernels()[0].1,
            "weights must be unchanged"
        );

        // s0=1,p=0 is an exact width-identical copy.
        let id = tab.width_corrected(1.0, 0.0, 10.0).unwrap();
        assert_eq!(id.kernels()[0].0, tab.kernels()[0].0);

        // p<0 -> higher energy is narrower (energy-dependent width).
        let wc2 = tab.width_corrected(1.0, -0.5, 10.0).unwrap();
        let (_, std_lo) = kernel_centroid_std(&wc2.kernels()[0].0, &wc2.kernels()[0].1); // 5 eV
        let (_, std_hi) = kernel_centroid_std(&wc2.kernels()[1].0, &wc2.kernels()[1].1); // 50 eV
        assert!(
            std_hi < std_lo,
            "p<0 should narrow higher E: {std_hi} !< {std_lo}"
        );
    }

    #[test]
    fn width_corrected_preserves_trapezoidal_centroid_on_nonuniform_grid() {
        // On a NON-uniform offset grid the trapezoidal-weighted centroid (what the
        // broadening integral actually integrates against) differs from the plain
        // centroid. The width scale must pivot about the former, else the fitted
        // width leaks into position. The uniform-grid test above cannot see this —
        // there the two centroids coincide.
        let offs = vec![-2.0, -1.5, 0.0, 0.5, 3.0]; // deliberately non-uniform
        let wts = vec![0.2, 0.6, 1.0, 0.7, 0.2];
        let tab =
            TabulatedResolution::from_kernels(vec![10.0], vec![(offs.clone(), wts.clone())], 25.0)
                .unwrap();

        // Trapezoidal-weighted centroid, mirroring broaden_presorted's dt weights.
        let trap_centroid = |o: &[f64], w: &[f64]| -> f64 {
            let n = o.len();
            let dt = |k: usize| -> f64 {
                if n <= 1 {
                    1.0
                } else if k == 0 {
                    o[1] - o[0]
                } else if k == n - 1 {
                    o[k] - o[k - 1]
                } else {
                    (o[k + 1] - o[k - 1]) * 0.5
                }
            };
            let (mut num, mut den) = (0.0, 0.0);
            for (k, (&oi, &wi)) in o.iter().zip(w).enumerate() {
                let tw = wi * dt(k).abs();
                num += oi * tw;
                den += tw;
            }
            num / den
        };
        let plain_centroid = |o: &[f64], w: &[f64]| -> f64 {
            o.iter().zip(w).map(|(a, b)| a * b).sum::<f64>() / w.iter().sum::<f64>()
        };

        let c_trap_before = trap_centroid(&offs, &wts);
        // Sanity: on this grid the trapezoidal and plain centroids genuinely differ,
        // so the test would fail under the old plain-centroid pivot.
        assert!(
            (c_trap_before - plain_centroid(&offs, &wts)).abs() > 1e-3,
            "test grid not non-uniform enough"
        );

        let wc = tab.width_corrected(2.0, 0.0, 10.0).unwrap();
        let c_trap_after = trap_centroid(&wc.kernels()[0].0, &wc.kernels()[0].1);
        assert!(
            (c_trap_after - c_trap_before).abs() < 1e-12,
            "integrated centroid leaked under width scale: {c_trap_before} -> {c_trap_after}"
        );
    }

    #[test]
    fn width_corrected_zero_weight_block_falls_back_to_zero_pivot() {
        // A degenerate all-zero-weight kernel has no centroid; the scale pivots
        // about 0 rather than dividing by a zero weight sum.
        let tab = TabulatedResolution::from_kernels(
            vec![10.0],
            vec![(vec![-1.0, 0.0, 2.0], vec![0.0, 0.0, 0.0])],
            25.0,
        )
        .unwrap();
        let wc = tab.width_corrected(2.0, 0.0, 10.0).unwrap();
        assert_eq!(wc.kernels()[0].0, vec![-2.0, 0.0, 4.0]); // scaled about 0
    }

    #[test]
    fn width_corrected_rejects_invalid_params() {
        // Public API: invalid whole-configuration inputs must hard-error up front
        // (not silently clamp), since a non-positive s0 reverses the offset order.
        let tab = TabulatedResolution::from_kernels(
            vec![10.0],
            vec![(vec![-1.0, 0.0, 2.0], vec![0.1, 1.0, 0.1])],
            25.0,
        )
        .unwrap();
        for (s0, p, e_ref) in [
            (0.0, 0.0, 10.0),          // s0 == 0
            (-1.0, 0.0, 10.0),         // s0 < 0
            (f64::NAN, 0.0, 10.0),     // s0 non-finite
            (1.0, f64::NAN, 10.0),     // p non-finite
            (1.0, 0.0, 0.0),           // e_ref == 0
            (1.0, 0.0, -5.0),          // e_ref < 0
            (1.0, 0.0, f64::INFINITY), // e_ref non-finite
        ] {
            assert!(
                matches!(
                    tab.width_corrected(s0, p, e_ref),
                    Err(ResolutionError::InvalidWidthCorrection { .. })
                ),
                "expected InvalidWidthCorrection for s0={s0}, p={p}, e_ref={e_ref}"
            );
        }
    }

    #[test]
    fn from_kernels_rejects_empty_kernel() {
        // An empty kernel passes the length-match check (0 == 0) but makes
        // broadening silently fall back to pass-through; reject it up front.
        let err = TabulatedResolution::from_kernels(vec![10.0], vec![(vec![], vec![])], 25.0);
        assert!(
            matches!(err, Err(ResolutionParseError::InvalidFormat(_))),
            "empty kernel should be rejected, got {err:?}"
        );
        // A non-empty kernel still constructs fine.
        assert!(
            TabulatedResolution::from_kernels(vec![10.0], vec![(vec![0.0], vec![1.0])], 25.0)
                .is_ok()
        );
        // Unsorted offsets are rejected — the broadener walks a monotonic two-
        // pointer bracket and derives trapezoidal widths assuming sorted offsets.
        let unsorted = TabulatedResolution::from_kernels(
            vec![10.0],
            vec![(vec![0.0, -1.0, 2.0], vec![0.2, 1.0, 0.2])],
            25.0,
        );
        assert!(
            matches!(unsorted, Err(ResolutionParseError::InvalidFormat(_))),
            "unsorted offsets should be rejected, got {unsorted:?}"
        );
        // Duplicate offsets (non-strict) are also rejected.
        let dup = TabulatedResolution::from_kernels(
            vec![10.0],
            vec![(vec![-1.0, 0.0, 0.0, 2.0], vec![0.1, 1.0, 1.0, 0.1])],
            25.0,
        );
        assert!(matches!(dup, Err(ResolutionParseError::InvalidFormat(_))));
    }

    #[test]
    fn from_text_rejects_unsorted_offsets_within_block() {
        // The second energy block has an out-of-order offset pair
        // (1.0 followed by 0.0): the trapezoidal `dt_width` quadrature
        // and the two-pointer bracket walk both assume sorted offsets,
        // so the parser must error rather than construct a
        // silently-corrupt kernel (same invariant `from_kernels`
        // enforces).
        let text = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
0.0 1.0
1.0 0.5

10.0 0.0
-1.0 0.1
1.0 0.5
0.0 1.0
";
        let err = TabulatedResolution::from_text(text, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("unsorted offsets must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("Kernel 1") && msg.contains("E = 10 eV"),
            "error must name the offending block and reference energy: {msg}"
        );
        // The same file with sorted offsets parses fine.
        let sorted = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
0.0 1.0
1.0 0.5

10.0 0.0
-1.0 0.1
0.0 1.0
1.0 0.5
";
        assert!(TabulatedResolution::from_text(sorted, 25.0).is_ok());
    }

    /// The remaining `from_kernels` invariants hold for parsed files too:
    /// empty energy blocks (silent pass-through at broaden time), non-finite
    /// kernel values (NaN norm bypasses the division guard), and non-finite
    /// reference energies (NaN compares false, slipping through the
    /// ascending check into the bracketing binary search) must all error.
    #[test]
    fn from_text_rejects_empty_block_and_non_finite_values() {
        // Empty block: header line for E=10 with no data lines.
        let empty_block = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
0.0 1.0

10.0 0.0

";
        let err = TabulatedResolution::from_text(empty_block, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("empty energy block must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("E = 10 eV") && msg.contains("no (offset, weight) points"),
            "error must name the empty block: {msg}"
        );

        // Non-finite kernel weight.
        let nan_weight = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
0.0 NaN
1.0 0.5
";
        let err = TabulatedResolution::from_text(nan_weight, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("non-finite weight must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("non-finite"),
            "error must name the non-finite value class: {msg}"
        );

        // Non-finite kernel OFFSET: must report the finiteness problem,
        // not the misleading "must be strictly ascending" (a NaN offset
        // fails the ascending comparison too — finiteness is checked
        // first so the message stays accurate).
        let nan_offset = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
NaN 1.0
1.0 0.5
";
        let err = TabulatedResolution::from_text(nan_offset, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("non-finite offset must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("non-finite") && !msg.contains("ascending"),
            "NaN offset must report finiteness, not sortedness: {msg}"
        );

        // Non-finite reference energy.
        let nan_energy = "\
Resolution file
---------------
NaN 0.0
-1.0 0.1
0.0 1.0
";
        let err = TabulatedResolution::from_text(nan_energy, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("non-finite reference energy must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("finite"),
            "error must name the finiteness requirement: {msg}"
        );
    }

    /// Non-positive reference energies must be rejected by BOTH
    /// constructors: the TOF map needs E > 0, and the width
    /// interpolation takes ln(E_ref) — a zero/negative reference would
    /// turn every blended weight into NaN, which bypasses the
    /// broadener's norm guard and silently disables broadening (a
    /// stray `0.0 0.0` line after a blank line in a VENUS/FTS file
    /// parses as an energy-block header).
    #[test]
    fn constructors_reject_non_positive_reference_energies() {
        let zero_energy = "\
Resolution file
---------------
0.0 0.0
-1.0 0.1
0.0 1.0
1.0 0.5

10.0 0.0
-1.0 0.1
0.0 1.0
1.0 0.5
";
        let err = TabulatedResolution::from_text(zero_energy, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("zero reference energy must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("positive"),
            "error must name the positivity requirement: {msg}"
        );

        let err = TabulatedResolution::from_kernels(
            vec![-10.0, 10.0],
            vec![
                (vec![-1.0, 0.0, 1.0], vec![0.1, 1.0, 0.1]),
                (vec![-1.0, 0.0, 1.0], vec![0.1, 1.0, 0.1]),
            ],
            25.0,
        );
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("negative reference energy must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("positive"),
            "error must name the positivity requirement: {msg}"
        );
    }

    /// Negative kernel weights (no physical meaning; the broadener
    /// skips them while the width moments would otherwise fold them
    /// in) must be rejected by BOTH constructors.
    #[test]
    fn constructors_reject_negative_weights() {
        let neg_weight = "\
Resolution file
---------------
5.0 0.0
-1.0 0.1
0.0 1.0
1.0 -0.2
";
        let err = TabulatedResolution::from_text(neg_weight, 25.0);
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("negative weight must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("negative weights"),
            "error must name the negative-weight problem: {msg}"
        );

        let err = TabulatedResolution::from_kernels(
            vec![10.0],
            vec![(vec![-1.0, 0.0, 1.0], vec![-1.0, 0.2, -1.0])],
            25.0,
        );
        let Err(ResolutionParseError::InvalidFormat(msg)) = err else {
            panic!("negative weights must be rejected, got {err:?}");
        };
        assert!(
            msg.contains("negative weights"),
            "error must name the negative-weight problem: {msg}"
        );
    }

    #[test]
    fn interpolated_kernel_blend_stays_ascending_and_mode_anchored() {
        // Two equal-length kernels with their peak at offset 0 — a loaded
        // UDR file's anchoring — through the width-normalized shape blend
        // that every between-reference energy takes. The blend must stay
        // strictly ascending (the sorted invariant the broadener relies on)
        // and leave whatever sits at 0 there.
        let off_lo = vec![-1.0, -0.4, 0.0, 0.6, 1.5, 3.0];
        let off_hi = vec![-0.5, -0.2, 0.0, 0.3, 0.8, 1.6]; // narrower, same length
        let wts = vec![0.1, 0.5, 1.0, 0.6, 0.3, 0.1]; // mode at index 2 (offset 0)
        let tab = TabulatedResolution::from_kernels(
            vec![5.0, 50.0],
            vec![(off_lo, wts.clone()), (off_hi, wts.clone())],
            25.0,
        )
        .unwrap();
        let (blended, weights) = tab.interpolated_kernel(15.0); // between 5 and 50 eV
        assert!(
            blended.windows(2).all(|w| w[0] < w[1]),
            "blended offsets not strictly ascending: {blended:?}"
        );
        let kmax = (0..weights.len())
            .max_by(|&a, &b| weights[a].total_cmp(&weights[b]))
            .unwrap();
        assert!(
            blended[kmax].abs() < 1e-9,
            "mode not anchored at offset 0: {}",
            blended[kmax]
        );
    }

    /// For a self-similar family (each block a width-scaled copy of one
    /// asymmetric template), the width-normalized shapes agree, so the
    /// blended kernel's trapezoidal width must equal the geometric
    /// interpolation `σ_lo·(σ_hi/σ_lo)^frac` — near-exactly (the scaled
    /// grids coincide point-for-point in z, so the merge degenerates to
    /// the template shape at the target width). An interior reference
    /// energy must return that block bitwise (exact-hit arm).
    ///
    /// SCOPE NOTE: this measures the blend with the implementation's
    /// own `trapezoidal_moments` and target formula, so it pins the
    /// merge machinery against its coded target (plus the chord
    /// separation below), not the physical width law independently —
    /// that independent, end-to-end anchor lives in
    /// `tests/kernel_width_interpolation.rs` and the pytest port,
    /// which measure the APPLIED broadening of parsed kernels against
    /// the analytic power law.
    #[test]
    fn interpolated_kernel_width_follows_power_law_for_self_similar_blocks() {
        let template_off = [-1.0, -0.4, 0.0, 0.8, 2.0, 4.0];
        let template_w = vec![0.05, 0.5, 1.0, 0.6, 0.2, 0.02];
        // σ ∝ E^{-1/2} scaling across refs 10 / 100 / 1000 eV.
        let scale = |e: f64| (e / 10.0f64).powf(-0.5);
        let block = |e: f64| -> (Vec<f64>, Vec<f64>) {
            (
                template_off.iter().map(|&o| o * scale(e)).collect(),
                template_w.clone(),
            )
        };
        let tab = TabulatedResolution::from_kernels(
            vec![10.0, 100.0, 1000.0],
            vec![block(10.0), block(100.0), block(1000.0)],
            25.0,
        )
        .unwrap();

        // Interior exact hit: the middle reference comes back bitwise.
        let (off_ref, w_ref) = test_support::interpolated_kernel(&tab, 100.0);
        let (exp_off, exp_w) = block(100.0);
        assert_eq!(off_ref.len(), exp_off.len());
        for (a, b) in off_ref.iter().zip(&exp_off) {
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "exact-hit offsets must be bitwise"
            );
        }
        for (a, b) in w_ref.iter().zip(&exp_w) {
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "exact-hit weights must be bitwise"
            );
        }

        // Between references: measured width equals the geometric law.
        let (off_lo, w_lo) = block(10.0);
        let (off_hi, w_hi) = block(100.0);
        let (_, s_lo) = trapezoidal_moments(&off_lo, &w_lo);
        let (_, s_hi) = trapezoidal_moments(&off_hi, &w_hi);
        for e in [16.0f64, 25.0, 40.0, 70.0] {
            let frac = (e.ln() - 10.0f64.ln()) / (100.0f64.ln() - 10.0f64.ln());
            let expected = s_lo * (s_hi / s_lo).powf(frac);
            let (offs, ws) = test_support::interpolated_kernel(&tab, e);
            let (_, got) = trapezoidal_moments(&offs, &ws);
            assert!(
                (got - expected).abs() / expected < 1e-9,
                "blended width at {e} eV: got {got}, expected {expected} \
                 (the pre-fix arithmetic chord gave {})",
                s_lo + frac * (s_hi - s_lo)
            );
            // Non-vacuity: the removed chord error is resolvable at
            // this tolerance (σ_hi/σ_lo = 10^{-1/2} → several % apart).
            assert!(
                (s_lo + frac * (s_hi - s_lo) - expected).abs() / expected > 1e-2,
                "fixture must separate chord from geometric law at {e} eV"
            );
        }
    }

    /// Identical bracketing blocks must come back bitwise-identical
    /// between the references (scale ratios exactly 1.0; the blend form
    /// `a + frac·(b − a)` is exact when a == b).
    #[test]
    fn interpolated_kernel_is_bitwise_identity_for_identical_blocks() {
        let off = vec![-1.0, -0.4, 0.0, 0.8, 2.0, 4.0];
        let w = vec![0.05, 0.5, 1.0, 0.6, 0.2, 0.02];
        let tab = TabulatedResolution::from_kernels(
            vec![5.0, 500.0],
            vec![(off.clone(), w.clone()), (off.clone(), w.clone())],
            25.0,
        )
        .unwrap();
        let (offs, ws) = test_support::interpolated_kernel(&tab, 42.0);
        assert_eq!(offs.len(), off.len());
        for (a, b) in offs.iter().zip(&off) {
            assert_eq!(a.to_bits(), b.to_bits(), "identity offsets must be bitwise");
        }
        for (a, b) in ws.iter().zip(&w) {
            assert_eq!(a.to_bits(), b.to_bits(), "identity weights must be bitwise");
        }
    }

    /// A NaN target energy (unreachable through validated public paths,
    /// but defended anyway) must yield a well-formed reference clone —
    /// never NaN offsets that break the broadener's invariants.
    #[test]
    fn interpolated_kernel_nan_energy_clamps_to_lowest_reference() {
        let off = vec![-1.0, 0.0, 2.0];
        let w = vec![0.3, 1.0, 0.2];
        let tab = TabulatedResolution::from_kernels(
            vec![10.0, 1000.0],
            vec![(off.clone(), w.clone()), (off.clone(), w.clone())],
            25.0,
        )
        .unwrap();
        let (offs, ws) = test_support::interpolated_kernel(&tab, f64::NAN);
        assert_eq!(offs, off);
        assert_eq!(ws, w);
    }

    /// Degenerate blocks (single-point → σ = 0) cannot be width-scaled;
    /// the nearer reference is cloned instead.
    #[test]
    fn interpolated_kernel_degenerate_blocks_fall_back_to_nearest() {
        let tab = TabulatedResolution::from_kernels(
            vec![10.0, 1000.0],
            vec![(vec![0.0], vec![1.0]), (vec![0.5], vec![1.0])],
            25.0,
        )
        .unwrap();
        // frac < 0.5 → lower block; frac > 0.5 → upper block.
        let (lo_off, _) = test_support::interpolated_kernel(&tab, 15.0);
        assert_eq!(lo_off, vec![0.0]);
        let (hi_off, _) = test_support::interpolated_kernel(&tab, 700.0);
        assert_eq!(hi_off, vec![0.5]);
    }

    // ── Smoke tests for the test_support oracles (`interp_spectrum` +
    //    `broaden_presorted_reference`).  The 7+ bit-exact tests below
    //    exercise the math thoroughly; these smoke tests just pin the
    //    boundary-condition return-shape behavior of the oracles so a
    //    future refactor that breaks empty-input or out-of-range
    //    handling fails loudly rather than only via bit-exact diffs.

    #[test]
    fn test_support_interp_spectrum_empty_returns_none() {
        assert_eq!(test_support::interp_spectrum(&[], &[], 1.0), None);
    }

    #[test]
    fn test_support_interp_spectrum_out_of_range_returns_none() {
        let energies = [1.0, 2.0, 3.0];
        let spectrum = [10.0, 20.0, 30.0];
        assert_eq!(
            test_support::interp_spectrum(&energies, &spectrum, 0.5),
            None
        );
        assert_eq!(
            test_support::interp_spectrum(&energies, &spectrum, 3.5),
            None
        );
    }

    #[test]
    fn test_support_broaden_presorted_reference_empty_returns_empty() {
        let tab = test_support::trivial_tabulated_resolution(25.0);
        let out = test_support::broaden_presorted_reference(&tab, &[], &[]);
        assert!(out.is_empty());
    }

    #[test]
    fn test_tof_factor_consistency() {
        // Verify our TOF_FACTOR matches the constants module.
        let e = 10.0; // eV
        let l = 25.0; // meters
        let tof_constants = constants::energy_to_tof(e, l);
        let tof_ours = TOF_FACTOR * l / e.sqrt();
        let rel_diff = (tof_constants - tof_ours).abs() / tof_constants;
        assert!(
            rel_diff < 1e-10,
            "TOF mismatch: constants={}, ours={}, diff={:.4}%",
            tof_constants,
            tof_ours,
            rel_diff * 100.0
        );
    }

    #[test]
    fn test_resolution_width_scaling() {
        let params = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();

        // Resolution width should increase with energy.
        let w1 = params.gaussian_width(1.0);
        let w10 = params.gaussian_width(10.0);
        let w100 = params.gaussian_width(100.0);

        assert!(w10 > w1, "Width should increase with energy");
        assert!(w100 > w10, "Width should increase with energy");

        // At low energies, timing dominates: ΔE ∝ E^(3/2)
        // At high energies, path dominates: ΔE ∝ E
        // The ratio ΔE(10)/ΔE(1) should be between 10 and 31.6 (= 10^1.5)
        let ratio = w10 / w1;
        assert!(
            ratio > 5.0 && ratio < 40.0,
            "Width ratio = {}, expected between 10 and 31.6",
            ratio
        );
    }

    #[test]
    fn test_zero_width_passthrough() {
        // If resolution parameters are zero, output should equal input.
        let energies = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let xs = vec![10.0, 20.0, 30.0, 20.0, 10.0];
        let params = ResolutionParams::new(25.0, 0.0, 0.0, 0.0).unwrap();
        let broadened = resolution_broaden(&energies, &xs, &params).unwrap();
        assert_eq!(broadened, xs);
    }

    #[test]
    fn test_broadening_reduces_peak() {
        // Resolution broadening should reduce peak heights and fill valleys.
        let n = 1001;
        let energies: Vec<f64> = (0..n).map(|i| 5.0 + (i as f64) * 0.01).collect();
        let center = 10.0;
        let gamma: f64 = 0.1; // Resonance width
        let xs: Vec<f64> = energies
            .iter()
            .map(|&e| {
                let de = e - center;
                1000.0 * (gamma / 2.0).powi(2) / (de * de + (gamma / 2.0).powi(2))
            })
            .collect();

        let params = ResolutionParams::new(25.0, 5.0, 0.01, 0.0).unwrap();
        let broadened = resolution_broaden(&energies, &xs, &params).unwrap();

        let orig_peak = xs.iter().cloned().fold(0.0_f64, f64::max);
        let broad_peak = broadened.iter().cloned().fold(0.0_f64, f64::max);

        assert!(
            broad_peak < orig_peak,
            "Broadened peak ({}) should be < original ({})",
            broad_peak,
            orig_peak
        );
        assert!(
            broad_peak > 1.0,
            "Broadened peak ({}) should still be substantial",
            broad_peak
        );
    }

    #[test]
    fn test_broadening_conserves_area() {
        // Resolution broadening should approximately conserve the area
        // under the cross-section curve.
        let n = 2001;
        let energies: Vec<f64> = (0..n).map(|i| 5.0 + (i as f64) * 0.005).collect();
        let center = 10.0;
        let gamma: f64 = 0.5;
        let xs: Vec<f64> = energies
            .iter()
            .map(|&e| {
                let de = e - center;
                1000.0 * (gamma / 2.0).powi(2) / (de * de + (gamma / 2.0).powi(2))
            })
            .collect();

        let params = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();
        let broadened = resolution_broaden(&energies, &xs, &params).unwrap();

        // Trapezoidal area
        let area_orig: f64 = (0..n - 1)
            .map(|i| 0.5 * (xs[i] + xs[i + 1]) * (energies[i + 1] - energies[i]))
            .sum();
        let area_broad: f64 = (0..n - 1)
            .map(|i| 0.5 * (broadened[i] + broadened[i + 1]) * (energies[i + 1] - energies[i]))
            .sum();

        let rel_diff = (area_orig - area_broad).abs() / area_orig;
        assert!(
            rel_diff < 0.02,
            "Area not conserved: orig={:.2}, broad={:.2}, rel_diff={:.4}",
            area_orig,
            area_broad,
            rel_diff
        );
    }

    #[test]
    fn test_gaussian_broadening_analytical() {
        // Broadening a Gaussian with a Gaussian should give a wider Gaussian.
        //
        // Input:  exp(-x²/(2σ₁²)) with σ₁ = 0.5 eV (standard Gaussian form)
        // Kernel: exp(-x²/W²) with W = 0.3 eV → std dev σ₂ = W/√2 = 0.2121 eV
        // Output: Gaussian with σ_out = √(σ₁² + σ₂²) = √(0.25 + 0.045) = 0.543 eV
        //
        // Note: kernel width varies slightly with energy (σ_E ∝ E for the
        // path-length contribution), so we allow ~5% tolerance.
        let n = 2001;
        let center = 10.0;
        let sigma_input = 0.5; // eV (standard deviation)
        let energies: Vec<f64> = (0..n).map(|i| 5.0 + (i as f64) * 0.005).collect();
        let xs: Vec<f64> = energies
            .iter()
            .map(|&e| {
                let de = e - center;
                1000.0 * (-de * de / (2.0 * sigma_input * sigma_input)).exp()
            })
            .collect();

        // Set delta_l such that W = gaussian_width(E=10) ≈ 0.3 eV.
        // W = 2·ΔL·E/L, so ΔL = W·L/(2E) = 0.3×25/(20) = 0.375 m
        let w_kernel = 0.3; // Kernel parameter W (exp(-x²/W²))
        let params =
            ResolutionParams::new(25.0, 0.0, w_kernel * 25.0 / (2.0 * center), 0.0).unwrap();

        // Verify kernel W at center energy
        let w_at_center = params.gaussian_width(center);
        assert!(
            (w_at_center - w_kernel).abs() / w_kernel < 0.01,
            "Kernel W at center: {}, expected {}",
            w_at_center,
            w_kernel
        );

        let broadened = resolution_broaden(&energies, &xs, &params).unwrap();

        // Kernel std dev = W/√2
        let sigma_kernel = w_kernel / 2.0_f64.sqrt();
        let sigma_expected = (sigma_input * sigma_input + sigma_kernel * sigma_kernel).sqrt();
        let fwhm_expected = 2.0 * (2.0_f64.ln() * 2.0).sqrt() * sigma_expected;

        // Measure FWHM from the broadened output
        let peak_idx = broadened
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap()
            .0;
        let peak_val = broadened[peak_idx];
        let half_max = peak_val / 2.0;

        let mut left_hm = energies[0];
        for i in (0..peak_idx).rev() {
            if broadened[i] < half_max {
                let t = (half_max - broadened[i]) / (broadened[i + 1] - broadened[i]);
                left_hm = energies[i] + t * (energies[i + 1] - energies[i]);
                break;
            }
        }
        let mut right_hm = energies[n - 1];
        for i in peak_idx..n - 1 {
            if broadened[i + 1] < half_max {
                let t = (half_max - broadened[i]) / (broadened[i + 1] - broadened[i]);
                right_hm = energies[i] + t * (energies[i + 1] - energies[i]);
                break;
            }
        }

        let fwhm_measured = right_hm - left_hm;
        let rel_err = (fwhm_measured - fwhm_expected).abs() / fwhm_expected;

        assert!(
            rel_err < 0.05,
            "FWHM: measured={:.4}, expected={:.4}, rel_err={:.2}%",
            fwhm_measured,
            fwhm_expected,
            rel_err * 100.0
        );
    }

    #[test]
    fn test_venus_typical_resolution() {
        // Verify resolution width for typical VENUS parameters.
        // VENUS: L ≈ 25 m, Δt ≈ 10 μs (pulsed source), ΔL ≈ 0.01 m
        let params = ResolutionParams::new(25.0, 10.0, 0.01, 0.0).unwrap();

        // At 1 eV: ΔE/E should be small (good resolution)
        let de_1 = params.gaussian_width(1.0);
        let de_over_e_1 = de_1 / 1.0;
        assert!(
            de_over_e_1 < 0.05,
            "ΔE/E at 1 eV = {:.4}, should be < 5%",
            de_over_e_1
        );

        // At 100 eV: resolution degrades
        let de_100 = params.gaussian_width(100.0);
        let de_over_e_100 = de_100 / 100.0;
        assert!(
            de_over_e_100 > de_over_e_1,
            "Resolution should degrade at higher energies"
        );
    }

    #[test]
    fn test_unsorted_energies_returns_error() {
        let energies = vec![1.0, 3.0, 2.0, 4.0]; // not sorted
        let xs = vec![10.0, 30.0, 20.0, 40.0];
        let params = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();
        let result = resolution_broaden(&energies, &xs, &params);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ResolutionError::UnsortedEnergies
        ));
    }

    #[test]
    fn test_length_mismatch_returns_error() {
        let energies = vec![1.0, 2.0, 3.0];
        let xs = vec![10.0, 20.0]; // wrong length
        let params = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();
        let result = resolution_broaden(&energies, &xs, &params);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ResolutionError::LengthMismatch {
                energies: 3,
                data: 2
            }
        ));
    }

    // --- ResolutionParams validation tests ---

    #[test]
    fn test_resolution_params_valid() {
        let p = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();
        assert!((p.flight_path_m() - 25.0).abs() < 1e-15);
        assert!((p.delta_t_us() - 1.0).abs() < 1e-15);
        assert!((p.delta_l_m() - 0.01).abs() < 1e-15);
    }

    #[test]
    fn test_resolution_params_rejects_zero_flight_path() {
        let err = ResolutionParams::new(0.0, 1.0, 0.01, 0.0).unwrap_err();
        assert_eq!(err, ResolutionParamsError::InvalidFlightPath(0.0));
    }

    #[test]
    fn test_resolution_params_rejects_negative_flight_path() {
        let err = ResolutionParams::new(-1.0, 1.0, 0.01, 0.0).unwrap_err();
        assert_eq!(err, ResolutionParamsError::InvalidFlightPath(-1.0));
    }

    #[test]
    fn test_resolution_params_rejects_nan_flight_path() {
        let err = ResolutionParams::new(f64::NAN, 1.0, 0.01, 0.0).unwrap_err();
        assert!(matches!(err, ResolutionParamsError::InvalidFlightPath(_)));
    }

    #[test]
    fn test_resolution_params_rejects_negative_delta_t() {
        let err = ResolutionParams::new(25.0, -1.0, 0.01, 0.0).unwrap_err();
        assert_eq!(err, ResolutionParamsError::InvalidDeltaT(-1.0));
    }

    #[test]
    fn test_resolution_params_rejects_nan_delta_t() {
        let err = ResolutionParams::new(25.0, f64::NAN, 0.01, 0.0).unwrap_err();
        assert!(matches!(err, ResolutionParamsError::InvalidDeltaT(_)));
    }

    #[test]
    fn test_resolution_params_rejects_negative_delta_l() {
        let err = ResolutionParams::new(25.0, 1.0, -0.01, 0.0).unwrap_err();
        assert_eq!(err, ResolutionParamsError::InvalidDeltaL(-0.01));
    }

    #[test]
    fn test_resolution_params_rejects_inf_delta_l() {
        let err = ResolutionParams::new(25.0, 1.0, f64::INFINITY, 0.0).unwrap_err();
        assert!(matches!(err, ResolutionParamsError::InvalidDeltaL(_)));
    }

    #[test]
    fn test_resolution_params_rejects_negative_delta_e() {
        let err = ResolutionParams::new(25.0, 1.0, 0.01, -0.05).unwrap_err();
        assert_eq!(err, ResolutionParamsError::InvalidDeltaE(-0.05));
    }

    #[test]
    fn test_resolution_params_rejects_nan_delta_e() {
        let err = ResolutionParams::new(25.0, 1.0, 0.01, f64::NAN).unwrap_err();
        assert!(matches!(err, ResolutionParamsError::InvalidDeltaE(_)));
    }

    #[test]
    fn test_resolution_params_accepts_zero_delta_e() {
        let p = ResolutionParams::new(25.0, 1.0, 0.01, 0.0).unwrap();
        assert!((p.delta_e_us() - 0.0).abs() < 1e-15);
        assert!(!p.has_exponential_tail());
    }

    // ─── broaden_presorted bit-exact equivalence harness ─────────────────────
    //
    // The optimized broaden_presorted uses a two-pointer walk instead of
    // binary search inside the inner convolution loop.  These tests pin
    // the math: the same formula, in the same order, must yield bit-exact
    // output against a canonical reference implementation that preserves
    // the pre-optimization code path.

    // The `interp_spectrum` + `broaden_presorted_reference` oracles
    // were promoted to `test_support` so the integration tests
    // (`tests/venus_usr_resolution{,_microbench}.rs`) share the same
    // byte-identical reference.  Imported below.
    use super::test_support::broaden_presorted_reference;

    /// Synthetic TabulatedResolution with 3 reference energies and a
    /// triangular kernel of varying widths.  Deterministic, no I/O.
    fn synthetic_tab_resolution() -> TabulatedResolution {
        fn triangle(width_us: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
            let half = width_us;
            let dt_step = 2.0 * half / (n - 1) as f64;
            let offsets: Vec<f64> = (0..n).map(|i| -half + i as f64 * dt_step).collect();
            let weights: Vec<f64> = offsets
                .iter()
                .map(|&dt| (1.0 - dt.abs() / half).max(0.0))
                .collect();
            (offsets, weights)
        }
        TabulatedResolution {
            ref_energies: Arc::new(vec![5.0, 50.0, 500.0]),
            kernels: Arc::new(vec![
                triangle(0.5, 31),
                triangle(1.0, 41),
                triangle(2.0, 51),
            ]),
            flight_path_m: 25.0,
        }
    }

    fn assert_bit_exact(reference: &[f64], actual: &[f64], label: &str) {
        assert_eq!(reference.len(), actual.len(), "{label}: length mismatch");
        for (i, (&a, &b)) in reference.iter().zip(actual.iter()).enumerate() {
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "{label}: element {i} mismatch: reference={a:.17e} actual={b:.17e}"
            );
        }
    }

    #[test]
    fn test_broaden_presorted_bit_exact_synthetic_uniform() {
        let tab = synthetic_tab_resolution();
        // Uniform log-spaced grid typical of VENUS analysis.
        let energies: Vec<f64> = (0..401).map(|i| 7.0 + i as f64 * 0.4825).collect();
        // Triangular dip + smooth background (resonance-like spectrum).
        let spectrum: Vec<f64> = energies
            .iter()
            .enumerate()
            .map(|(i, &e)| 0.9 - 0.7 * (-((e - 50.0).powi(2) / 4.0)).exp() + 0.001 * (i as f64))
            .collect();

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "synthetic_uniform");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_synthetic_nonuniform() {
        let tab = synthetic_tab_resolution();
        // Non-uniform: denser near 6.674 eV (resonance-like), sparser far away.
        let energies: Vec<f64> = {
            let mut e = Vec::new();
            for i in 0..200 {
                e.push(5.0 + (i as f64) * 0.05);
            }
            for i in 0..100 {
                e.push(15.0 + (i as f64) * 0.5);
            }
            for i in 0..50 {
                e.push(65.0 + (i as f64) * 2.0);
            }
            e
        };
        let spectrum: Vec<f64> = energies
            .iter()
            .map(|&e| 1.0 - 0.5 * (-((e - 6.674).powi(2) / 0.1)).exp())
            .collect();

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "synthetic_nonuniform");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_constant_spectrum() {
        // Constant spectrum must pass through unchanged (within trapezoid
        // normalization) — preserves integral exactly.
        let tab = synthetic_tab_resolution();
        let energies: Vec<f64> = (0..501).map(|i| 1.0 + i as f64 * 0.5).collect();
        let spectrum = vec![0.42f64; energies.len()];

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "constant_spectrum");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_short_grid() {
        // Edge case: 2-point grid.  Exercises the smallest grid that has
        // a valid (lo, hi) bracket — tests bracket_hi bounds handling.
        let tab = synthetic_tab_resolution();
        let energies = vec![10.0, 12.0];
        let spectrum = vec![0.5, 0.8];

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "short_grid");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_single_point_grid() {
        // Edge case: 1-point grid.  Exercises the n == 1 early-return
        // pass-through guard that the optimized path adds (no bracket
        // available for interpolation).
        let tab = synthetic_tab_resolution();
        let energies = vec![10.0];
        let spectrum = vec![0.5];

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "single_point_grid");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_exact_equality_target() {
        // Regression: exercise the tie-break case where `e_prime` lands
        // exactly on a grid point.  The kernel has a point at dt=0, so
        // `e_prime == energies[i]` exactly at the center kernel offset
        // for every target `i`.  The optimized path must match the
        // reference's upper-bound binary-search semantics bit-exactly.
        let tab = synthetic_tab_resolution();
        // Irregular-spacing grid so the spectrum interp at the equality
        // point isn't trivially reducible to the input value.
        let mut energies: Vec<f64> = Vec::new();
        let mut e = 3.0f64;
        for k in 0..800 {
            energies.push(e);
            e += 0.05 + 0.01 * (k as f64).sin();
        }
        // Spectrum with large local gradient so `a + (b - a)` vs `b`
        // would diverge at 1 ULP if the tie-break were wrong.
        let spectrum: Vec<f64> = energies
            .iter()
            .map(|&e| 1.0e10 * (-(((e - 6.0) / 0.2).powi(2))).exp() + 1.0e-10 * e)
            .collect();

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "exact_equality_target");
    }

    #[test]
    fn test_broaden_presorted_bit_exact_random_spectrum() {
        // Random spectrum with varied magnitudes exercises the interpolation
        // arithmetic across sign changes and scales.
        let tab = synthetic_tab_resolution();
        let energies: Vec<f64> = (0..1001).map(|i| 2.0 + i as f64 * 0.2).collect();
        // Deterministic pseudo-random via a simple LCG (no external dep).
        let mut state: u64 = 0xDEAD_BEEF_CAFE_BABE;
        let spectrum: Vec<f64> = energies
            .iter()
            .map(|_| {
                state = state
                    .wrapping_mul(6364136223846793005)
                    .wrapping_add(1442695040888963407);
                let f = ((state >> 33) as f64) / (u32::MAX as f64);
                f * 2.0 - 1.0
            })
            .collect();

        let reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        let actual = tab.broaden_presorted(&energies, &spectrum);
        assert_bit_exact(&reference, &actual, "random_spectrum");
    }

    // ---------------------------------------------------------------
    // `crates/nereids-physics/tests/venus_usr_resolution.rs`
    // (`test_broaden_presorted_bit_exact_on_venus_usr`) — see issue
    //. The integration test parses a synthetic SAMMY USR-format
    // kernel via `common::synthetic_venus_usr_tab()` (the real VENUS
    // BL10 fixture is not approved for public release; ).
    // ---------------------------------------------------------------

    #[test]
    fn test_plan_reuse_bit_exact_across_multiple_spectra() {
        // Core promise of ResolutionPlan: building the plan once and
        // applying it to K different spectra must yield the same output
        // as K independent `broaden_presorted` calls.
        let tab = synthetic_tab_resolution();
        let energies: Vec<f64> = (0..401).map(|i| 7.0 + i as f64 * 0.4825).collect();

        // Build plan ONCE.
        let plan = tab.plan(&energies).expect("sorted grid must validate");
        assert_eq!(plan.len(), energies.len());
        assert_eq!(plan.target_energies(), &energies[..]);

        // Apply across 5 varied spectra.
        let mut state: u64 = 0xCAFE_BABE_DEAD_BEEF;
        for spec_idx in 0..5 {
            let spectrum: Vec<f64> = energies
                .iter()
                .enumerate()
                .map(|(i, &e)| {
                    state = state
                        .wrapping_mul(6364136223846793005)
                        .wrapping_add(1442695040888963407);
                    let noise = ((state >> 33) as f64) / (u32::MAX as f64);
                    // Varied magnitudes and shapes per spectrum to catch
                    // spectrum-dependent arithmetic drift.
                    (10.0f64).powi(spec_idx - 2) * (1.0 - 0.5 * noise)
                        + 0.3 * (-((e - 50.0).powi(2) / 4.0)).exp()
                        + (spec_idx as f64) * 1e-8 * (i as f64)
                })
                .collect();

            let via_plan = plan.apply(&spectrum);
            let via_reference = broaden_presorted_reference(&tab, &energies, &spectrum);
            assert_bit_exact(
                &via_reference,
                &via_plan,
                &format!("plan_reuse[spec_idx={spec_idx}]"),
            );
        }
    }

    #[test]
    fn test_plan_passthrough_cases() {
        // n == 0, n == 1, and e <= 0.0 must all produce the same
        // passthrough behaviour via plan as via broaden_presorted.
        let tab = synthetic_tab_resolution();

        // n == 0: empty plan, empty result.
        let plan = tab.plan(&[]).unwrap();
        assert_eq!(plan.len(), 0);
        assert!(plan.is_empty());
        let out: Vec<f64> = plan.apply(&[]);
        assert!(out.is_empty());

        // n == 1: passthrough for any spectrum value.
        let plan1 = tab.plan(&[5.0]).unwrap();
        assert_eq!(plan1.len(), 1);
        let out1 = plan1.apply(&[0.42]);
        assert_eq!(out1, vec![0.42]);

        // e <= 0.0 in the middle of a grid: passthrough at that index.
        // Mixed positive / non-positive energies are pathological but
        // the current implementation handles them, and the plan must
        // match.  Grid is still non-descending (0.0 ≤ 10.0 etc.) so
        // plan() accepts it.
        let energies = vec![1.0, 1.0, 10.0, 100.0];
        let spectrum = vec![0.1, 0.5, 0.9, 0.3];
        let via_plan = {
            let plan = tab.plan(&energies).unwrap();
            plan.apply(&spectrum)
        };
        let via_reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        assert_bit_exact(
            &via_reference,
            &via_plan,
            "mixed_positive_and_zero_energies",
        );
    }

    #[test]
    fn test_plan_rejects_unsorted_energies() {
        // `broaden()` rejects unsorted grids via validate_inputs; `plan()`
        // must do the same so a caller doesn't silently build a plan with
        // misbracketed e_prime lookups and then produce wrong σ output
        // from `ResolutionPlan::apply`.
        let tab = synthetic_tab_resolution();
        let result = tab.plan(&[10.0, 1.0, 100.0]);
        assert!(matches!(result, Err(ResolutionError::UnsortedEnergies)));
    }

    #[test]
    fn test_plan_apply_is_nan_safe_at_degenerate_bracket() {
        // When two adjacent target energies are equal (span = 0), the
        // plan encodes `frac = 0.0` and the apply path must short-
        // circuit to `spectrum[lo]` without reading `spectrum[lo+1]`.
        // A NaN at the upper bracket would propagate through
        // `0.0 * NaN = NaN` and corrupt the result otherwise.
        let tab = synthetic_tab_resolution();
        // Grid has a degenerate duplicate at indices 1 and 2.
        let energies = vec![8.0, 10.0, 10.0, 12.0, 50.0, 100.0];
        // Spectrum with NaN exactly at the upper-bracket index (2) that
        // the degenerate pair maps to; any retained (target, kernel-
        // point) entry whose `e_prime` lands inside that duplicate
        // bracket MUST NOT pull the NaN into the output.
        let spectrum = vec![0.1, 0.5, f64::NAN, 0.9, 0.2, 0.05];
        let plan = tab.plan(&energies).unwrap();
        let via_plan = plan.apply(&spectrum);
        let via_reference = broaden_presorted_reference(&tab, &energies, &spectrum);
        // Both paths must agree on the non-pathological targets.  The
        // reference path returns `spectrum[lo]` directly in the
        // degenerate case (no touch of `spectrum[lo+1]`) and the plan
        // path's `frac == 0.0` short-circuit matches bit-exactly.
        // Targets whose kernel legitimately interpolates across index 2
        // will pull the NaN in BOTH paths equally — that's physics, not
        // a bug — so we compare bit-pattern with a NaN-aware helper.
        assert_eq!(via_plan.len(), via_reference.len());
        for (i, (&a, &b)) in via_reference.iter().zip(via_plan.iter()).enumerate() {
            // Both NaN or both finite and bit-exact.
            if a.is_nan() {
                assert!(b.is_nan(), "plan[{i}]={b} but reference is NaN");
            } else {
                assert_eq!(
                    a.to_bits(),
                    b.to_bits(),
                    "nan_safe mismatch at {i}: reference={a} plan={b}"
                );
            }
        }
    }

    #[test]
    fn test_plan_apply_exact_match_frac_plus_zero_propagates_nan() {
        // Regression gate for a subtle sign-of-zero short-circuit bug.
        //
        // When `e_prime` aligns EXACTLY with a grid point `energies[lo]`,
        // `plan_presorted`'s interp fraction computes to `+0.0`, yet the
        // bracket is NOT degenerate (span is a normal positive float).
        // In that case `broaden_presorted` still evaluates
        //   s = spectrum[lo] + (+0.0) * (spectrum[lo+1] - spectrum[lo])
        // which, for `spectrum[lo+1] = NaN`, reads `0.0 * NaN = NaN` and
        // propagates `NaN` into `s`.  The earlier `frac == 0.0` branch
        // in `ResolutionPlan::apply` incorrectly treated this case as
        // degenerate (since `+0.0 == -0.0` under `==`) and short-circuited
        // to `spectrum[lo]`, producing a finite output where the scalar
        // reference produced NaN.
        //
        // Fix: `plan_presorted` now stores `-0.0` (negative-signed zero)
        // for the degenerate sentinel, and `apply` disambiguates via
        // `to_bits()`, so the non-degenerate `+0.0` path correctly reads
        // `spectrum[lo+1]` and propagates NaN.
        let tab = synthetic_tab_resolution();

        // Grid has a point at energy = 10.0.  We engineer a target grid
        // where the broadened kernel at one of the targets produces an
        // `e_prime` that aligns exactly with `energies[lo]` of one of its
        // retained entries.  Achieved by building a coarse target grid
        // and letting the two-pointer walk land on an exact match on at
        // least one (target, kernel-point) pair.
        let energies: Vec<f64> = (0..32).map(|i| 1.0 + i as f64).collect();
        // Spectrum with NaN scattered at multiple lo+1 indices.  At
        // least one retained plan entry in this synthetic configuration
        // will have `frac == +0.0` from an exact-match case, which must
        // propagate NaN in apply.
        let mut spectrum = vec![0.5_f64; energies.len()];
        for v in &mut spectrum[3..] {
            *v = f64::NAN;
        }
        let plan = tab.plan(&energies).unwrap();
        let via_plan = plan.apply(&spectrum);
        let via_reference = broaden_presorted_reference(&tab, &energies, &spectrum);

        // Bit-exact equivalence on all targets, including NaN-propagated
        // ones.  This test would FAIL pre-fix (plan returns finite where
        // reference returns NaN for any exact-match plan entry with a
        // NaN at `lo+1`).
        assert_eq!(via_plan.len(), via_reference.len());
        for (i, (&a, &b)) in via_reference.iter().zip(via_plan.iter()).enumerate() {
            if a.is_nan() {
                assert!(
                    b.is_nan(),
                    "target {i}: reference produced NaN (NaN propagated through \
                     exact-match `frac = +0.0` path) but plan returned finite {b}",
                );
            } else {
                assert_eq!(
                    a.to_bits(),
                    b.to_bits(),
                    "target {i}: reference={a} plan={b}"
                );
            }
        }
    }

    #[test]
    #[should_panic(expected = "must match plan target-grid length")]
    fn test_plan_apply_spectrum_length_mismatch_panics() {
        let tab = synthetic_tab_resolution();
        let plan = tab.plan(&[1.0, 2.0, 3.0]).unwrap();
        // Wrong spectrum length — caller error should panic with a
        // clear message rather than silently producing garbage.
        let _ = plan.apply(&[0.1, 0.2]);
    }

    // ─── apply_resolution_with_plan / _presorted_with_plan dispatch harness ───
    //
    // These gates cover the public/crate-visible wrappers added for
    // production plan-caching.  Every production caller (fit-model
    // layer, spatial dispatch) goes through one of these two entries;
    // their dispatch choices must be byte-identical to the non-plan
    // paths they replace.

    #[test]
    fn test_apply_resolution_with_plan_tabulated_matches_non_plan_path() {
        let tab = synthetic_tab_resolution();
        let resolution = ResolutionFunction::Tabulated(Arc::new(tab.clone()));
        let energies: Vec<f64> = (0..128).map(|i| 1.0 + i as f64 * (200.0 / 128.0)).collect();
        let spectrum: Vec<f64> = energies
            .iter()
            .map(|&e| 1.0 - 0.3 * (-((e - 20.0).powi(2) / 4.0)).exp())
            .collect();

        let baseline = apply_resolution(&energies, &spectrum, &resolution).unwrap();

        let plan = build_resolution_plan(&energies, &resolution).unwrap();
        assert!(
            plan.is_some(),
            "tabulated resolution must produce Some(plan)"
        );
        let planned =
            apply_resolution_with_plan(plan.as_ref(), &energies, &spectrum, &resolution).unwrap();
        assert_eq!(planned.len(), baseline.len());
        for (i, (&a, &b)) in baseline.iter().zip(planned.iter()).enumerate() {
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "apply_resolution_with_plan mismatch at {i}: baseline={a} planned={b}"
            );
        }
    }

    #[test]
    fn test_apply_resolution_with_plan_gaussian_returns_none_plan_and_matches() {
        let resolution =
            ResolutionFunction::Gaussian(ResolutionParams::new(25.0, 1.0e-3, 0.02, 0.01).unwrap());
        let energies: Vec<f64> = (0..64).map(|i| 1.0 + i as f64 * 3.0).collect();
        let spectrum: Vec<f64> = energies.iter().map(|&e| 1.0 / e).collect();

        let plan = build_resolution_plan(&energies, &resolution).unwrap();
        assert!(
            plan.is_none(),
            "Gaussian resolution must not produce a plan"
        );

        let baseline = apply_resolution(&energies, &spectrum, &resolution).unwrap();
        let planned =
            apply_resolution_with_plan(plan.as_ref(), &energies, &spectrum, &resolution).unwrap();
        assert_eq!(planned.len(), baseline.len());
        for (i, (&a, &b)) in baseline.iter().zip(planned.iter()).enumerate() {
            assert_eq!(
                a.to_bits(),
                b.to_bits(),
                "gaussian fallback mismatch at {i}: baseline={a} planned={b}"
            );
        }
    }

    #[test]
    fn test_apply_resolution_with_plan_rejects_same_length_different_grid() {
        // `p.len() == energies.len()` is necessary
        // but not sufficient.  A plan built for one grid and applied
        // to a different same-length grid would silently gather
        // spectrum values at brackets belonging to the original grid
        // — wrong σ output without any error surfaced.  The grid-
        // identity check in `apply_resolution_with_plan` guards this
        // failure mode and reports the first differing index.
        let tab = synthetic_tab_resolution();
        let resolution = ResolutionFunction::Tabulated(Arc::new(tab.clone()));
        let energies_plan: Vec<f64> = (0..32).map(|i| 1.0 + i as f64).collect();
        let mut energies_apply = energies_plan.clone();
        // Perturb a single interior point so lengths still match.
        energies_apply[5] += 0.25;
        let spectrum = vec![0.7; energies_apply.len()];

        let plan = tab.plan(&energies_plan).unwrap();
        let result =
            apply_resolution_with_plan(Some(&plan), &energies_apply, &spectrum, &resolution);
        match result {
            Err(ResolutionError::PlanGridMismatch { first_diff_index }) => {
                assert_eq!(first_diff_index, 5);
            }
            other => panic!("expected PlanGridMismatch, got {:?}", other),
        }
    }

    #[test]
    fn test_apply_resolution_with_plan_rejects_length_mismatch() {
        let tab = synthetic_tab_resolution();
        let resolution = ResolutionFunction::Tabulated(Arc::new(tab.clone()));
        let energies_plan: Vec<f64> = (0..32).map(|i| 1.0 + i as f64).collect();
        let energies_apply: Vec<f64> = (0..48).map(|i| 1.0 + i as f64).collect();
        let spectrum = vec![0.5; energies_apply.len()];

        let plan = tab.plan(&energies_plan).unwrap();
        let result =
            apply_resolution_with_plan(Some(&plan), &energies_apply, &spectrum, &resolution);
        match result {
            Err(ResolutionError::LengthMismatch { energies, data }) => {
                assert_eq!(energies, 48);
                assert_eq!(data, 32);
            }
            other => panic!("expected LengthMismatch, got {:?}", other),
        }
    }

    #[test]
    fn test_build_resolution_plan_rejects_unsorted_energies_for_gaussian() {
        // Gaussian returns None on success, but must still reject an
        // unsorted grid — callers use `build_resolution_plan` to
        // centralise the sort-check so the downstream `apply` path can
        // skip it.
        let resolution =
            ResolutionFunction::Gaussian(ResolutionParams::new(25.0, 1.0e-3, 0.02, 0.01).unwrap());
        let result = build_resolution_plan(&[3.0, 1.0, 2.0], &resolution);
        assert!(matches!(result, Err(ResolutionError::UnsortedEnergies)));
    }

    // ---------------------------------------------------------------
    // `crates/nereids-physics/tests/venus_usr_resolution_microbench.rs`
    // (`test_broaden_presorted_bench`, `test_plan_reuse_bench`,
    // `resolution_matrix_apply_microbench`) — see. They
    // parse a synthetic SAMMY USR-format kernel via
    // `common::synthetic_venus_usr_tab()` (the real VENUS BL10
    // fixture is not approved for public release; ).
    // ---------------------------------------------------------------

    // ---------- ResolutionMatrix (CSR compile) tests ----------
    //
    // CI-hermetic synthetic tests — use hand-constructed
    // `ResolutionPlan`s via `make_synthetic_plan`; no fixture
    // dependency, run on every `cargo test` invocation.  Cover
    // passthrough rows, `-0.0` sentinel rows, regular linear-interp
    // rows, CSR invariants, and the non-finite contract exclusion.
    //
    // End-to-end equivalence tests against the VENUS-like USR
    // operator (synthetic SAMMY-format kernel) at realistic grid
    // sizes (512, 3471) live in
    // `crates/nereids-physics/tests/venus_usr_resolution.rs` — see
    // issues and.

    /// Hybrid abs+rel tolerance used across equivalence tests.  Guards
    /// against the `a ≈ 0` trap where `a.abs().max(1e-300)` produces
    /// meaningless relative errors for genuinely-zero reference values.
    fn max_hybrid_err(a: &[f64], b: &[f64]) -> f64 {
        a.iter()
            .zip(b)
            .map(|(x, y)| {
                let denom = x.abs().max(y.abs()).max(1e-12);
                (x - y).abs() / denom
            })
            .fold(0.0_f64, f64::max)
    }

    /// Build a synthetic multi-row plan with realistic overlap
    /// patterns — used as a CI-hermetic stand-in for the VENUS
    /// kernel.  Each target row `i` draws weights from a triangular
    /// kernel around column `i`, normalized so the row is
    /// row-stochastic.  `half_kernel` controls the spread.
    fn make_synthetic_overlap_plan(n_grid: usize, half_kernel: usize) -> ResolutionPlan {
        assert!(n_grid > 2 * half_kernel, "grid too small for kernel");
        let energies: Vec<f64> = (0..n_grid).map(|i| 10.0 + i as f64).collect();
        let mut rows: Vec<SyntheticRow> = Vec::with_capacity(n_grid);
        for i in 0..n_grid {
            let lo_min = i.saturating_sub(half_kernel);
            // Clamp so `lo ∈ [0, n_grid - 2]` — the linear-interp
            // branch reads `spec[lo + 1]`, and the `-0.0` sentinel is
            // the only way to safely go up to `lo = n_grid - 1`.  We
            // keep all synthetic entries on the regular branch here.
            let lo_max = (i + half_kernel).min(n_grid - 2);
            let entries: Vec<SyntheticEntry> = (lo_min..=lo_max)
                .map(|lo| {
                    let d = (lo as i64 - i as i64).abs() as f64;
                    let w = 1.0 - d / (half_kernel as f64 + 1.0);
                    // A uniform `frac = 0.5` distributes each entry's
                    // weight evenly across `lo` and `lo + 1`, which
                    // exercises the regular linear-interp branch of
                    // `compile_to_matrix`.
                    SyntheticEntry {
                        lo: lo as u32,
                        frac: 0.5,
                        weight: w,
                    }
                })
                .collect();
            let norm: f64 = entries.iter().map(|e| e.weight).sum();
            rows.push(SyntheticRow { entries, norm });
        }
        make_synthetic_plan(energies, rows)
    }

    /// CI-hermetic: row-stochasticity on a synthetic multi-row plan.
    #[test]
    fn resolution_matrix_is_row_stochastic_synthetic() {
        let plan = make_synthetic_overlap_plan(40, 5);
        let matrix = plan.compile_to_matrix();
        for i in 0..matrix.len() {
            let start = matrix.row_starts()[i] as usize;
            let end = matrix.row_starts()[i + 1] as usize;
            let row_sum: f64 = matrix.values()[start..end].iter().sum();
            assert!(
                (row_sum - 1.0).abs() < 1e-13,
                "row {} sum = {} (expected 1.0 within 1e-13)",
                i,
                row_sum,
            );
        }
    }

    /// CI-hermetic: equivalence of `apply_r` and `plan.apply` on a
    /// synthetic multi-row plan, 40-point grid, half-kernel 5.
    #[test]
    fn resolution_matrix_apply_equivalent_to_plan_apply_synthetic() {
        let plan = make_synthetic_overlap_plan(40, 5);
        let matrix = plan.compile_to_matrix();
        // Beer-Lambert-shaped synthetic spectrum, bounded [0, 1].
        let spec: Vec<f64> = (0..matrix.len())
            .map(|i| {
                let x = i as f64 / 39.0;
                1.0 - 0.7 * (-((x - 0.5).powi(2)) / 0.01).exp()
            })
            .collect();
        let plan_out = plan.apply(&spec);
        let matrix_out = apply_r(&matrix, &spec);
        let max_err = max_hybrid_err(&plan_out, &matrix_out);
        assert!(
            max_err < 1e-12,
            "synthetic apply_r vs plan.apply max hybrid err = {:.3e} (expected < 1e-12)",
            max_err,
        );
    }

    /// CI-hermetic: CSR column indices strictly ascending per row on
    /// a synthetic multi-row plan.
    #[test]
    fn resolution_matrix_csr_column_indices_sorted_per_row_synthetic() {
        let plan = make_synthetic_overlap_plan(30, 4);
        let matrix = plan.compile_to_matrix();
        for i in 0..matrix.len() {
            let start = matrix.row_starts()[i] as usize;
            let end = matrix.row_starts()[i + 1] as usize;
            let row_cols = &matrix.col_indices()[start..end];
            for w in row_cols.windows(2) {
                assert!(
                    w[0] < w[1],
                    "row {} col_indices not strictly ascending: {:?}",
                    i,
                    row_cols,
                );
            }
        }
    }

    /// CI-hermetic: grid-mismatch / length-mismatch detection via
    /// `apply_resolution_with_matrix` on a synthetic plan.
    #[test]
    fn resolution_matrix_grid_and_length_mismatch_synthetic() {
        let plan = make_synthetic_overlap_plan(16, 3);
        let matrix = plan.compile_to_matrix();
        let n = matrix.len();
        let energies: Vec<f64> = (0..n).map(|i| 10.0 + i as f64).collect();
        let spec = vec![1.0_f64; n];

        // Same grid + length → passes.
        assert!(apply_resolution_with_matrix(&energies, &matrix, &spec).is_ok());

        // Perturb one energy → MatrixGridMismatch with offending
        // index.
        let mut mutated = energies.clone();
        mutated[7] += 1e-12;
        let err = apply_resolution_with_matrix(&mutated, &matrix, &spec)
            .expect_err("grid mismatch must error");
        assert_eq!(
            err,
            ResolutionError::MatrixGridMismatch {
                first_diff_index: 7,
            }
        );

        // Short spectrum → LengthMismatch.
        let short = vec![1.0_f64; n - 1];
        let err = apply_resolution_with_matrix(&energies, &matrix, &short)
            .expect_err("length mismatch must error");
        assert!(matches!(err, ResolutionError::LengthMismatch { .. }));
    }

    // ---------------------------------------------------------------
    // `crates/nereids-physics/tests/venus_usr_resolution.rs`
    // (`resolution_matrix_is_row_stochastic_on_venus_kernel`,
    //  `resolution_matrix_apply_equivalent_to_plan_apply_on_venus_kernel`,
    //  `resolution_matrix_apply_equivalent_at_production_grid`,
    //  `resolution_matrix_apply_equivalent_across_densities`,
    //  `resolution_matrix_csr_column_indices_sorted_per_row`,
    //  `resolution_matrix_grid_mismatch_detected`,
    //  `resolution_matrix_length_mismatch_detected`) — see issues
    // and. They parse a synthetic SAMMY USR-format kernel
    // via `common::synthetic_venus_usr_tab()`.
    // ---------------------------------------------------------------

    #[test]
    fn resolution_matrix_empty_plan() {
        // Compile must not panic and must produce a valid empty
        // matrix when the plan itself is empty.  Build the empty
        // plan synthetically (no fixture needed) — an empty
        // `target_energies` plus empty `norm` / `starts = [0]`
        // yields the same zero-row plan that
        // `TabulatedResolution::plan(&[])` would produce.
        let plan = make_synthetic_plan(Vec::new(), Vec::new());
        let matrix = plan.compile_to_matrix();
        assert_eq!(matrix.len(), 0);
        assert!(matrix.is_empty());
        assert_eq!(matrix.nnz(), 0);
    }

    /// Hand-construct a `ResolutionPlan` that deliberately exercises
    /// both the passthrough branch (`norm ≤ DIVISION_FLOOR`) and the
    /// `-0.0` degenerate-bracket sentinel — neither of which is
    /// reached on the VENUS fixture at the tested grid sizes, which
    /// made the earlier fixture-based passthrough test vacuous.  This
    /// replacement verifies the two unreached branches with direct
    /// assertions on the resulting CSR.
    fn make_synthetic_plan(target_energies: Vec<f64>, rows: Vec<SyntheticRow>) -> ResolutionPlan {
        let n = target_energies.len();
        assert_eq!(rows.len(), n);
        let mut starts: Vec<u32> = Vec::with_capacity(n + 1);
        starts.push(0);
        let mut lo_idx: Vec<u32> = Vec::new();
        let mut frac: Vec<f64> = Vec::new();
        let mut weight: Vec<f64> = Vec::new();
        let mut norm: Vec<f64> = Vec::with_capacity(n);
        for row in &rows {
            norm.push(row.norm);
            for entry in &row.entries {
                lo_idx.push(entry.lo);
                frac.push(entry.frac);
                weight.push(entry.weight);
            }
            starts.push(lo_idx.len() as u32);
        }
        ResolutionPlan {
            target_energies,
            starts,
            lo_idx,
            frac,
            weight,
            norm,
        }
    }

    struct SyntheticRow {
        entries: Vec<SyntheticEntry>,
        norm: f64,
    }

    struct SyntheticEntry {
        lo: u32,
        frac: f64,
        weight: f64,
    }

    #[test]
    fn resolution_matrix_passthrough_row_compiles_to_identity_entry() {
        // Row 0: passthrough via norm ≤ DIVISION_FLOOR.
        // Row 1: regular linear-interp entry (lo=1 → reads cols 1, 2).
        // Row 2: degenerate `-0.0` sentinel entry (lo=2 → reads col 2 only).
        //
        // Grid has 4 cells so `lo ∈ [0, n-2] = [0, 2]` holds for all
        // entries — this preserves the `ResolutionPlan::apply` SAFETY
        // invariant that `lo + 1 < n` even if a future refactor
        // weakens the `-0.0` sentinel short-circuit.
        let plan = make_synthetic_plan(
            vec![10.0, 20.0, 30.0, 40.0],
            vec![
                SyntheticRow {
                    entries: vec![],
                    // 0.0 is <= DIVISION_FLOOR, so row 0 goes through
                    // the passthrough branch.
                    norm: 0.0,
                },
                SyntheticRow {
                    entries: vec![SyntheticEntry {
                        lo: 1,
                        frac: 0.25,
                        weight: 1.0,
                    }],
                    norm: 1.0,
                },
                SyntheticRow {
                    entries: vec![SyntheticEntry {
                        lo: 2,
                        frac: -0.0,
                        weight: 1.0,
                    }],
                    norm: 1.0,
                },
                // Row 3: passthrough too, to round out the 4-cell grid.
                SyntheticRow {
                    entries: vec![],
                    norm: 0.0,
                },
            ],
        );
        let matrix = plan.compile_to_matrix();

        // Row 0 — single (0, 0, 1.0).
        let r0_start = matrix.row_starts()[0] as usize;
        let r0_end = matrix.row_starts()[1] as usize;
        assert_eq!(r0_end - r0_start, 1, "passthrough row must have 1 entry");
        assert_eq!(matrix.col_indices()[r0_start], 0);
        assert_eq!(matrix.values()[r0_start].to_bits(), 1.0_f64.to_bits());

        // Row 1 — linear-interp: contributes at col 1 and col 2.
        let r1_start = matrix.row_starts()[1] as usize;
        let r1_end = matrix.row_starts()[2] as usize;
        assert_eq!(
            r1_end - r1_start,
            2,
            "linear-interp row must have 2 entries"
        );
        assert_eq!(matrix.col_indices()[r1_start], 1);
        assert_eq!(matrix.col_indices()[r1_start + 1], 2);
        assert!((matrix.values()[r1_start] - 0.75).abs() < 1e-14);
        assert!((matrix.values()[r1_start + 1] - 0.25).abs() < 1e-14);

        // Row 2 — `-0.0` sentinel: single entry at col 2 (no col 3).
        let r2_start = matrix.row_starts()[2] as usize;
        let r2_end = matrix.row_starts()[3] as usize;
        assert_eq!(
            r2_end - r2_start,
            1,
            "-0.0 sentinel row must have exactly 1 entry (not 2)",
        );
        assert_eq!(matrix.col_indices()[r2_start], 2);
        assert_eq!(matrix.values()[r2_start].to_bits(), 1.0_f64.to_bits());

        // Cross-check with apply semantics: spec[3] is chosen so the
        // sentinel row, if buggy, would contaminate the output.
        // Both `plan.apply` and `apply_r` must ignore spec[3] at
        // row 2.
        let spec = vec![7.0, 11.0, 13.0, 999.0];
        let plan_out = plan.apply(&spec);
        let matrix_out = apply_r(&matrix, &spec);
        // Row 0 passthrough: out[0] = spec[0] = 7.
        assert!((matrix_out[0] - 7.0).abs() < 1e-14);
        assert!((plan_out[0] - 7.0).abs() < 1e-14);
        // Row 1: 0.75 * spec[1] + 0.25 * spec[2] = 0.75*11 + 0.25*13 = 11.5.
        assert!((matrix_out[1] - 11.5).abs() < 1e-14);
        assert!((plan_out[1] - 11.5).abs() < 1e-14);
        // Row 2 sentinel: 1.0 * spec[2] = 13 — NOT 999 (would indicate
        // spec[lo+1] was read).
        assert!((matrix_out[2] - 13.0).abs() < 1e-14);
        assert!((plan_out[2] - 13.0).abs() < 1e-14);
        // Row 3 passthrough: out[3] = spec[3] = 999.
        assert!((matrix_out[3] - 999.0).abs() < 1e-14);
        assert!((plan_out[3] - 999.0).abs() < 1e-14);
    }

    /// Documents (and guards) the explicit contract exclusion on
    /// non-finite spectra between `ResolutionPlan::apply` and
    /// `apply_r`.  See [`ResolutionPlan::compile_to_matrix`] docstring
    /// for the full reasoning; this test simply pins the divergence
    /// so a future unification attempt fails loudly.
    #[test]
    fn resolution_matrix_nonfinite_contract() {
        // 3-cell grid so `lo = 0` for the regular row reads cols 0, 1
        // and the sentinel row at `lo = 1` reads col 1 only — `lo ∈
        // [0, n-2] = [0, 1]` satisfied.
        let plan = make_synthetic_plan(
            vec![10.0, 20.0, 30.0],
            vec![
                SyntheticRow {
                    entries: vec![SyntheticEntry {
                        lo: 0,
                        frac: 0.5,
                        weight: 1.0,
                    }],
                    norm: 1.0,
                },
                SyntheticRow {
                    entries: vec![SyntheticEntry {
                        lo: 1,
                        frac: -0.0, // sentinel: short-circuit to spec[lo]
                        weight: 1.0,
                    }],
                    norm: 1.0,
                },
                SyntheticRow {
                    entries: vec![],
                    norm: 0.0, // passthrough
                },
            ],
        );
        let matrix = plan.compile_to_matrix();

        // Spectrum with same-sign infinities in both bins of row 0's
        // non-degenerate bracket.
        let inf_spec = vec![f64::INFINITY, f64::INFINITY, 0.0];
        let plan_out = plan.apply(&inf_spec);
        let matrix_out = apply_r(&matrix, &inf_spec);

        // Row 0: plan.apply evaluates `s_lo + frac * (s_hi - s_lo)`
        // = `+∞ + 0.5 * (+∞ - +∞)` = `+∞ + 0.5 * NaN` = NaN.
        // apply_r evaluates `0.5 * +∞ + 0.5 * +∞` = `+∞`.
        assert!(plan_out[0].is_nan(), "plan.apply must produce NaN on ∞+∞");
        assert!(matrix_out[0].is_infinite(), "apply_r collapses ∞+∞ to ∞");

        // Row 1 (sentinel): both paths short-circuit to spec[lo] = ∞,
        // so there is no divergence here.
        assert!(plan_out[1].is_infinite());
        assert!(matrix_out[1].is_infinite());
    }

    /// Documents (and guards) the analogous
    /// divergence on **finite spectra near f64 overflow**.  With
    /// opposite-sign neighboring bins at f64::MAX, `plan.apply`'s
    /// `s_lo + frac * (s_hi - s_lo)` overflows in the subtraction
    /// and returns `±∞`, while `apply_r`'s `(1 - frac) * s_lo +
    /// frac * s_hi` stays finite because the overflow is avoided by
    /// scaling before summation.  This is why the equivalence
    /// contract on [`ResolutionPlan::compile_to_matrix`] is scoped
    /// to bounded finite spectra (Beer-Lambert `T ∈ [0, 1]`) — no
    /// production forward model can hit this case.
    #[test]
    fn resolution_matrix_large_finite_contract() {
        let plan = make_synthetic_plan(
            vec![10.0, 20.0, 30.0],
            vec![
                SyntheticRow {
                    entries: vec![SyntheticEntry {
                        lo: 0,
                        frac: 0.5,
                        weight: 1.0,
                    }],
                    norm: 1.0,
                },
                SyntheticRow {
                    entries: vec![],
                    norm: 0.0, // passthrough
                },
                SyntheticRow {
                    entries: vec![],
                    norm: 0.0,
                },
            ],
        );
        let matrix = plan.compile_to_matrix();

        // Opposite-sign large finite bins at row 0's non-degenerate
        // bracket.  `s_hi - s_lo = -f64::MAX - f64::MAX = -∞`.
        let big_spec = vec![f64::MAX, -f64::MAX, 0.0];
        let plan_out = plan.apply(&big_spec);
        let matrix_out = apply_r(&matrix, &big_spec);

        // plan.apply: s_lo + frac * (s_hi - s_lo) = MAX + 0.5 * (-∞)
        // = MAX + -∞ = -∞.
        assert!(
            plan_out[0].is_infinite() && plan_out[0] < 0.0,
            "plan.apply must overflow to -∞ on opposite-sign MAX bins; got {}",
            plan_out[0],
        );
        // apply_r: 0.5 * MAX + 0.5 * -MAX = 0.
        assert!(
            matrix_out[0].is_finite(),
            "apply_r must stay finite (scaled before summation); got {}",
            matrix_out[0],
        );
        assert!(matrix_out[0].abs() < 1e-280);
    }

    // ------------------------------------------------------------------
    // TabulatedResolution::kernel_support_ev — used by SAMMY EMIN/EMAX
    // -equivalent fit-energy-range margin computation.
    // ------------------------------------------------------------------

    /// `synthetic_tab_resolution` uses triangular kernels whose
    /// outermost entries (`weight = 1 - |dt|/half`) are exactly zero
    /// at `dt = ±half`; the actual non-zero support is the next-
    /// outermost entry at `±half · (1 - 1/(n-1))`.
    fn triangle_dt_max(half: f64, n: usize) -> f64 {
        // dt_step = 2*half / (n-1); next-outermost = half - dt_step
        half - 2.0 * half / (n - 1) as f64
    }

    /// Exact-map expected support: the larger of the up-side excursion
    /// `E·((t/(t−dt⁺))²−1)` and the down-side `E·(1−(t/(t+dt⁻))²)`
    /// with `t = TOF_FACTOR·L/√E` — the same map the broadener applies
    /// per kernel point, so this oracle is exact by construction.
    fn exact_support(e: f64, dt_pos: f64, dt_neg: f64, l: f64) -> f64 {
        let t = TOF_FACTOR * l / e.sqrt();
        let up = e * ((t / (t - dt_pos)).powi(2) - 1.0);
        let down = e * (1.0 - (t / (t + dt_neg)).powi(2));
        up.max(down)
    }

    /// At a reference energy with a known kernel half-width in TOF, the
    /// support must equal the exact TOF→E excursion of the outermost
    /// non-zero offsets.  At E = 50 eV the triangle kernel has
    /// half = 1.0 μs, n = 41, so the largest non-zero offset is
    /// `1.0 · (1 − 1/40) = 0.975` on both sides.
    #[test]
    fn test_tabulated_kernel_support_at_ref_energy_matches_exact_map() {
        let r = synthetic_tab_resolution();
        let e: f64 = 50.0;
        let dt_max = triangle_dt_max(1.0, 41);
        let expected = exact_support(e, dt_max, dt_max, 25.0);
        let got = r.kernel_support_ev(e);
        assert!(
            (got - expected).abs() / expected < 1e-12,
            "support at ref energy: got {got}, expected {expected}"
        );
    }

    /// Between two reference energies the support tracks the ACTUAL
    /// width-interpolated kernel: it must cover that kernel's non-zero
    /// offsets (non-circular — the blend comes from
    /// `interpolated_kernel` itself), while sitting strictly BELOW the
    /// old take-the-wider-bracket bound (proving the interior arm
    /// engaged rather than falling back to per-kernel extremes).
    #[test]
    fn test_tabulated_kernel_support_covers_actual_blend_between_refs() {
        let r = synthetic_tab_resolution();
        let e: f64 = 100.0; // between the 50 eV and 500 eV references
        let got = r.kernel_support_ev(e);

        // Cover: exact excursion of the blended kernel's w>0 extremes.
        let (offs, ws) = test_support::interpolated_kernel(&r, e);
        let dt_pos = offs
            .iter()
            .zip(&ws)
            .filter(|&(_, &w)| w > 0.0)
            .map(|(&o, _)| o)
            .fold(0.0f64, f64::max);
        let dt_neg = offs
            .iter()
            .zip(&ws)
            .filter(|&(_, &w)| w > 0.0)
            .map(|(&o, _)| -o)
            .fold(0.0f64, f64::max);
        let actual_excursion = exact_support(e, dt_pos, dt_neg, 25.0);
        assert!(
            got >= actual_excursion * (1.0 - 1e-12),
            "support must cover the actual blended kernel: got {got}, \
             actual excursion {actual_excursion}"
        );

        // Tightness + non-vacuity: strictly below the pre-blend bound
        // built from the wider 500 eV bracket's extreme (1.96 µs) —
        // the between-ref kernel is genuinely narrower.
        let old_bound = exact_support(e, triangle_dt_max(2.0, 51), triangle_dt_max(2.0, 51), 25.0);
        assert!(
            got < old_bound,
            "interior support must track the narrower interpolated \
             kernel: got {got}, old wider-bracket bound {old_bound}"
        );
    }

    /// Below the lowest ref energy: use the lowest ref kernel.
    /// Above the highest ref energy: use the highest ref kernel.
    #[test]
    fn test_tabulated_kernel_support_uses_nearest_outside_grid() {
        let r = synthetic_tab_resolution();
        // Below grid (ref_min = 5 eV; triangle(half=0.5, n=31)).
        let e_low: f64 = 1.0;
        let dt_low = triangle_dt_max(0.5, 31);
        let exp_low = exact_support(e_low, dt_low, dt_low, 25.0);
        assert!((r.kernel_support_ev(e_low) - exp_low).abs() / exp_low < 1e-12);
        // Above grid (ref_max = 500 eV; triangle(half=2.0, n=51)).
        let e_hi: f64 = 1000.0;
        let dt_hi = triangle_dt_max(2.0, 51);
        let exp_hi = exact_support(e_hi, dt_hi, dt_hi, 25.0);
        assert!((r.kernel_support_ev(e_hi) - exp_hi).abs() / exp_hi < 1e-12);
    }

    /// The exact up-side excursion strictly exceeds the linear
    /// chain-rule estimate for a wide positive (delayed-emission) tail
    /// at high energy — the case where the old linear margin
    /// under-covered exactly the side the convolution gather loads.
    #[test]
    fn test_tabulated_kernel_support_exceeds_linear_estimate_for_wide_tail() {
        let offsets = vec![-1.0, 0.0, 15.0];
        let weights = vec![0.3, 1.0, 0.2];
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![100.0]),
            kernels: Arc::new(vec![(offsets, weights)]),
            flight_path_m: 25.0,
        };
        let e: f64 = 100.0;
        let linear = 2.0 * e.powf(1.5) / (TOF_FACTOR * 25.0) * 15.0;
        let got = r.kernel_support_ev(e);
        assert!(
            got > linear,
            "exact support must exceed the linear estimate on the \
             high-E side: got {got}, linear {linear}"
        );
        let expected = exact_support(e, 15.0, 1.0, 25.0);
        assert!(
            (got - expected).abs() / expected < 1e-12,
            "exact support: got {got}, expected {expected}"
        );
    }

    /// An offset at or past the nominal flight time contributes no reach: the
    /// support is that of the surviving points.
    #[test]
    fn test_tabulated_kernel_support_counts_only_points_the_broadening_keeps() {
        let offsets = vec![0.0, 2.0, 10.0];
        let weights = vec![1.0, 0.5, 0.5];
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![100.0]),
            kernels: Arc::new(vec![(offsets, weights)]),
            flight_path_m: 25.0,
        };
        // Choose E high enough that 2 μs < t = K·L/√E ≤ 10 μs, so the
        // 10 μs point is dropped and the 2 μs point is the last kept.
        let t_at = |e: f64| TOF_FACTOR * 25.0 / e.sqrt();
        let mut e: f64 = 100.0;
        while t_at(e) > 10.0 {
            e *= 10.0;
        }
        let t = t_at(e);
        assert!(t > 2.0, "fixture must keep the 2 μs point (t = {t})");
        let expected_above = e * ((t / (t - 2.0)).powi(2) - 1.0);
        let (_, high) = r.gather_bounds_ev(e);
        let above = high - e;
        assert!(above.is_finite(), "reach must be finite, got {above}");
        assert!(
            (above - expected_above).abs() <= 1e-9 * expected_above,
            "reach {above} should be that of the last surviving offset, {expected_above}"
        );
    }

    /// Between two references an offset past the flight time drops only
    /// itself: the reach is that of the last surviving offset of the blended
    /// kernel.
    #[test]
    fn test_tabulated_kernel_support_between_refs_keeps_surviving_offsets() {
        let offsets = vec![0.0, 10.0, 20.0, 80.0];
        let weights = vec![1.0; 4];
        let r = TabulatedResolution::from_kernels(
            vec![1.0, 4.0],
            vec![(offsets.clone(), weights.clone()), (offsets, weights)],
            1.0,
        )
        .expect("valid two-block table");
        let e: f64 = 2.0;
        let t = TOF_FACTOR / e.sqrt();
        assert!(
            t > 20.0 && t < 80.0,
            "fixture must drop the 80 μs point and keep the 20 μs point (t = {t})"
        );
        let expected_above = e * ((t / (t - 20.0)).powi(2) - 1.0);
        let (_, high) = r.gather_bounds_ev(e);
        let above = high - e;
        assert!(
            (above - expected_above).abs() <= 1e-9 * expected_above,
            "reach {above} should be that of the last surviving offset, {expected_above}"
        );
    }

    /// Non-positive / non-finite energy → 0.0 (no broadening footprint).
    #[test]
    fn test_tabulated_kernel_support_returns_zero_for_invalid_energy() {
        let r = synthetic_tab_resolution();
        assert_eq!(r.kernel_support_ev(0.0), 0.0);
        assert_eq!(r.kernel_support_ev(-1.0), 0.0);
        assert_eq!(r.kernel_support_ev(f64::NAN), 0.0);
        assert_eq!(r.kernel_support_ev(f64::INFINITY), 0.0);
    }

    /// Zero-weight tail entries must not inflate the support; only
    /// `weights[i] > 0` entries count.
    #[test]
    fn test_tabulated_kernel_support_ignores_zero_weight_entries() {
        // Build a kernel where the outermost entries have weight 0.
        let offsets = vec![-10.0, -1.0, 0.0, 1.0, 10.0];
        let weights = vec![0.0, 0.5, 1.0, 0.5, 0.0];
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![100.0]),
            kernels: Arc::new(vec![(offsets, weights)]),
            flight_path_m: 25.0,
        };
        let e: f64 = 100.0;
        // Expected support uses dt = ±1.0 (the outermost zero-weight
        // entries at ±10 are ignored), not ±10.0.
        let expected = exact_support(e, 1.0, 1.0, 25.0);
        let got = r.kernel_support_ev(e);
        assert!(
            (got - expected).abs() / expected < 1e-12,
            "support should ignore zero-weight entries: got {got}, expected {expected}"
        );
    }

    /// Between reference kernels, the blended shape is positive on the
    /// FRINGE between a block's outermost `w > 0` entry and its
    /// adjacent `w == 0` entry (linear interpolation), and a merged
    /// point from the other block can land there — so the support's
    /// closure extremes (outermost positive weight extended to the
    /// adjacent zero-weight entry) must cover the actual blended
    /// kernel, which reaches beyond both blocks' bare `w > 0` maxima.
    #[test]
    fn test_tabulated_kernel_support_covers_blend_activated_offsets() {
        // Kernel A's positive support ends at z ≈ 3.5 (offset 5,
        // σ_A ≈ 1.44) with a zero-weight fringe out to z ≈ 7.0
        // (offset 10). Kernel B is compact (σ_B ≈ 0.97) with a
        // low-weight point at z ≈ 5.2 (offset 5) — inside A's fringe
        // after width normalization — so the blended kernel is
        // positive beyond A's bare w>0 extreme.
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![10.0, 1000.0]),
            kernels: Arc::new(vec![
                (vec![0.0, 5.0, 10.0, 20.0], vec![1.0, 0.1, 0.0, 0.0]),
                (vec![0.0, 1.0, 5.0, 6.0], vec![1.0, 0.8, 0.05, 0.0]),
            ]),
            flight_path_m: 25.0,
        };
        let e: f64 = 100.0; // strictly between the reference energies
        let got = r.kernel_support_ev(e);

        // Non-circular cover: the actual blended kernel's w>0 extremes.
        let (offs, ws) = test_support::interpolated_kernel(&r, e);
        let dt_pos = offs
            .iter()
            .zip(&ws)
            .filter(|&(_, &w)| w > 0.0)
            .map(|(&o, _)| o)
            .fold(0.0f64, f64::max);
        let dt_neg = offs
            .iter()
            .zip(&ws)
            .filter(|&(_, &w)| w > 0.0)
            .map(|(&o, _)| -o)
            .fold(0.0f64, f64::max);
        let actual_excursion = exact_support(e, dt_pos, dt_neg, 25.0);
        assert!(
            got >= actual_excursion * (1.0 - 1e-12),
            "support must cover the actual blended kernel (incl. the \
             zero-weight fringe): got {got}, actual {actual_excursion}"
        );

        // The fringe matters: the actual blend reaches beyond block
        // A's bare positive maximum (5 µs) scaled to the target width
        // — assert the blend truly is wider than a no-fringe reading
        // of block A would suggest, keeping this case load-bearing.
        let (_, s_lo) = trapezoidal_moments(&r.kernels[0].0, &r.kernels[0].1);
        let (_, s_hi) = trapezoidal_moments(&r.kernels[1].0, &r.kernels[1].1);
        let frac = (e.ln() - 10.0f64.ln()) / (1000.0f64.ln() - 10.0f64.ln());
        let s_t = s_lo * (s_hi / s_lo).powf(frac);
        let bare_positive_bound = 5.0 / s_lo * s_t;
        assert!(
            dt_pos > bare_positive_bound,
            "blend must extend into the zero-weight fringe: dt_pos {dt_pos}, \
             bare-positive bound {bare_positive_bound}"
        );
    }

    /// The support must contain the footprint of the ACTUAL broadener:
    /// broadening a flat baseline with a single narrow dip must leave
    /// every target untouched whose window
    /// `[e − support(e), e + support(e)]` excludes the dip.
    /// Non-circular by construction — the oracle here is `broaden`
    /// itself, not the `exact_support` formula mirror.
    #[test]
    fn test_tabulated_kernel_support_contains_broadener_footprint() {
        // Asymmetric kernel with a dominant delayed (positive) tail.
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![100.0]),
            kernels: Arc::new(vec![(vec![-1.0, 0.0, 6.0], vec![0.2, 1.0, 0.7])]),
            flight_path_m: 25.0,
        };
        // Dense uniform grid; flat baseline with a single-point dip.
        let de = 0.05;
        let energies: Vec<f64> = (0..2001).map(|j| 50.0 + de * j as f64).collect();
        let f = 1000; // dip mid-grid, at ~100 eV
        let e_f = energies[f];
        let mut spectrum = vec![1.0; energies.len()];
        spectrum[f] = 0.0;
        let out = r.broaden(&energies, &spectrum).unwrap();

        // Piecewise-linear reads touch spectrum[f] only for
        // e′ ∈ (energies[f−1], energies[f+1]); require one extra grid
        // step of margin so bracketing/interpolation edge effects
        // cannot straddle the window boundary.
        let margin = 2.0 * de;
        let (mut excluded_below, mut excluded_above) = (0usize, 0usize);
        let mut included_differs = false;
        let mut upside_differs = false;
        for (&e, &o) in energies.iter().zip(out.iter()) {
            let s = r.kernel_support_ev(e);
            if e + s + margin < e_f || e - s - margin > e_f {
                assert!(
                    (o - 1.0).abs() < 1e-12,
                    "target {e} eV (support {s}) must be untouched by a \
                     dip at {e_f} eV outside its window; got {o}"
                );
                if e < e_f {
                    excluded_below += 1;
                } else {
                    excluded_above += 1;
                }
            } else if (o - 1.0).abs() > 1e-3 {
                included_differs = true;
                if e < e_f - margin {
                    upside_differs = true;
                }
            }
        }
        // Non-vacuity: both exclusion regions were exercised, the dip
        // measurably alters at least one in-window target, and the
        // delayed tail reaches the dip from BELOW — the direction the
        // convolution gather loads.
        assert!(
            excluded_below > 0 && excluded_above > 0,
            "grid must exercise both exclusion regions \
             (below: {excluded_below}, above: {excluded_above})"
        );
        assert!(
            included_differs,
            "the dip must measurably alter at least one in-window target"
        );
        assert!(
            upside_differs,
            "the delayed tail must reach the dip from a target below it"
        );
    }

    /// Same containment property, exercised through the
    /// BETWEEN-REFERENCES support arm: two width-scaled reference
    /// blocks bracket the grid, so every target energy uses the
    /// width-interpolated kernel and the closure-extent support bound.
    /// The oracle is `broaden` itself — fully independent of both the
    /// support formula and the interpolation implementation.
    #[test]
    fn test_tabulated_kernel_support_contains_broadener_footprint_between_refs() {
        let r = TabulatedResolution {
            ref_energies: Arc::new(vec![10.0, 1000.0]),
            kernels: Arc::new(vec![
                (vec![-2.0, 0.0, 12.0], vec![0.2, 1.0, 0.7]),
                (vec![-0.5, 0.0, 3.0], vec![0.2, 1.0, 0.7]),
            ]),
            flight_path_m: 25.0,
        };
        let de = 0.05;
        let energies: Vec<f64> = (0..2001).map(|j| 50.0 + de * j as f64).collect();
        let f = 1000; // dip mid-grid, at ~100 eV — between the refs
        let e_f = energies[f];
        let mut spectrum = vec![1.0; energies.len()];
        spectrum[f] = 0.0;
        let out = r.broaden(&energies, &spectrum).unwrap();

        let margin = 2.0 * de;
        let (mut excluded, mut included_differs) = (0usize, false);
        for (&e, &o) in energies.iter().zip(out.iter()) {
            let s = r.kernel_support_ev(e);
            if e + s + margin < e_f || e - s - margin > e_f {
                assert!(
                    (o - 1.0).abs() < 1e-12,
                    "between-refs target {e} eV (support {s}) must be \
                     untouched by a dip at {e_f} eV outside its window; got {o}"
                );
                excluded += 1;
            } else if (o - 1.0).abs() > 1e-3 {
                included_differs = true;
            }
        }
        assert!(
            excluded > 0,
            "grid must exercise the exclusion region between references"
        );
        assert!(
            included_differs,
            "the dip must measurably alter at least one in-window target"
        );
    }

    #[test]
    fn piecewise_linear_bin_probabilities_pins_delta_and_partial_window() {
        let delta = |t: f64, edges: &[f64]| {
            piecewise_linear_bin_probabilities(&[t], &[3.5], edges)
                .expect("one-point kernel is a unit delta under normalization")
        };
        assert_eq!(delta(2.0, &[0.0, 1.0, 3.0, 5.0]), vec![0.0, 1.0, 0.0]);
        // The last bin is right-closed: a delta exactly on the final edge
        // belongs to it.
        assert_eq!(delta(5.0, &[0.0, 1.0, 3.0, 5.0]), vec![0.0, 0.0, 1.0]);
        // A delta outside every bin contributes nothing.
        assert_eq!(delta(9.0, &[0.0, 1.0, 3.0, 5.0]), vec![0.0, 0.0, 0.0]);
        let w3 = [0.0, 1.0, 0.0];
        assert!(
            piecewise_linear_bin_probabilities(
                &[0.0, 1.0, 1.0, 2.0],
                &[0.0, 1.0, 1.0, 0.0],
                &[0.5, 1.5]
            )
            .is_none()
        );
        for bad_times in [[0.0, f64::NAN, 2.0], [0.0, 1.0, f64::INFINITY]] {
            assert!(piecewise_linear_bin_probabilities(&bad_times, &w3, &[0.5, 1.5]).is_none());
        }
        for bad_edges in [[0.5, f64::NAN], [1.5, 0.5]] {
            assert!(
                piecewise_linear_bin_probabilities(&[0.0, 1.0, 2.0], &w3, &bad_edges).is_none()
            );
        }

        // Support normalization: a triangle on [0, 2] integrates to one over
        // its full support, and a partial window reports the true fraction.
        let times = [0.0, 1.0, 2.0];
        let weights = [0.0, 4.0, 0.0]; // arbitrary scale — normalization removes it
        let full = piecewise_linear_bin_probabilities(&times, &weights, &[0.0, 2.0])
            .expect("triangle integrates");
        assert!((full[0] - 1.0).abs() < 1e-15);
        let halves = piecewise_linear_bin_probabilities(&times, &weights, &[0.0, 0.5, 1.0])
            .expect("triangle integrates");
        assert!((halves[0] - 0.125).abs() < 1e-15);
        assert!((halves[1] - 0.375).abs() < 1e-15);
    }
}

#[cfg(test)]
mod width_convention_tests {
    use super::*;

    /// Grid half-span, in units of the nominal width. Beyond 6 the remaining
    /// Gaussian mass is below 1e-16, so 12 is already generous.
    const SPAN_IN_WIDTHS: usize = 12;

    /// Grid points per nominal width.
    ///
    /// Two errors scale as `(step/W)²` here: the trapezoidal second-moment
    /// quadrature, and the variance the one-bin-wide impulse carries in its own
    /// right (`step²/12`). At 50 samples per width both are ~1e-4 of the
    /// measured σ, two orders below the tolerances these tests assert, while
    /// the convolution cost scales as the square of the point count.
    const SAMPLES_PER_WIDTH: usize = 50;

    /// The kernel's second moment, measured numerically rather than taken from
    /// the width parameter the kernel was built from.
    ///
    /// Broadens a unit impulse and reads the standard deviation back off the
    /// result. `nominal_width_ev` only sizes the grid — the measurement itself
    /// never uses it, so the result cannot agree with `gaussian_width` by
    /// construction. That independence is the point of the test.
    fn measured_sigma_ev(params: &ResolutionParams, center_ev: f64, nominal_width_ev: f64) -> f64 {
        let n: usize = 2 * SPAN_IN_WIDTHS * SAMPLES_PER_WIDTH + 1;
        let half_span_ev = SPAN_IN_WIDTHS as f64 * nominal_width_ev;
        let step = 2.0 * half_span_ev / (n - 1) as f64;
        let energies: Vec<f64> = (0..n)
            .map(|i| center_ev - half_span_ev + i as f64 * step)
            .collect();
        // A unit impulse at the centre bin: broadening it returns the kernel.
        let mut impulse = vec![0.0; n];
        impulse[(n - 1) / 2] = 1.0 / step;

        let kernel = resolution_broaden(&energies, &impulse, params).expect("broadening runs");

        let mass: f64 = kernel.iter().sum::<f64>() * step;
        assert!(
            (mass - 1.0).abs() < 1.0e-3,
            "kernel is not normalised on this span: mass {mass}"
        );
        let mean: f64 = kernel
            .iter()
            .zip(&energies)
            .map(|(k, e)| k * e * step)
            .sum::<f64>()
            / mass;
        let variance: f64 = kernel
            .iter()
            .zip(&energies)
            .map(|(k, e)| k * (e - mean) * (e - mean) * step)
            .sum::<f64>()
            / mass;
        variance.sqrt()
    }

    /// `gaussian_width` returns W, and W/√2 is the standard deviation.
    ///
    /// This is the assertion that fixes the convention. If the width were a
    /// standard deviation instead, the measured σ would come back a factor √2
    /// larger than W/√2 and this fails by 41 %.
    #[test]
    fn gaussian_width_is_a_w_parameter_not_a_standard_deviation() {
        let center = 10.0;
        let params = ResolutionParams::new(25.0, 1.0, 0.0, 0.0).expect("valid params");
        let w = params.gaussian_width(center);
        assert!(w > 0.0, "test is vacuous without a width");

        let measured = measured_sigma_ev(&params, center, w);
        let expected = w / SQRT_2;
        assert!(
            (measured - expected).abs() / expected < 2.0e-3,
            "measured σ {measured:.9} but W/√2 is {expected:.9} (W = {w:.9})"
        );
        // And it is NOT the standard deviation itself, by a clear margin.
        assert!(
            (measured - w).abs() / w > 0.25,
            "measured σ {measured:.9} is indistinguishable from W {w:.9}"
        );
    }

    /// `fwhm()` is the full width at half maximum of the kernel it describes.
    #[test]
    fn fwhm_matches_the_kernel_measured_half_maximum() {
        let center = 10.0;
        let params = ResolutionParams::new(25.0, 1.0, 0.0, 0.0).expect("valid params");
        let w = params.gaussian_width(center);
        let measured_sigma = measured_sigma_ev(&params, center, w);
        // FWHM of a Gaussian in terms of its own standard deviation.
        let fwhm_from_measured = 2.0 * (2.0 * 2.0_f64.ln()).sqrt() * measured_sigma;
        let reported = params.fwhm(center);
        assert!(
            (reported - fwhm_from_measured).abs() / fwhm_from_measured < 2.0e-3,
            "fwhm() reports {reported:.9} but the kernel measures {fwhm_from_measured:.9}"
        );
    }

    /// `FWHM_PER_W` is exactly what the previous open-coded expression produced.
    #[test]
    fn fwhm_per_w_constant_is_bit_exact() {
        assert_eq!(
            FWHM_PER_W.to_bits(),
            (2.0 * (2.0_f64.ln()).sqrt()).to_bits(),
            "the named constant changed the value it replaced"
        );
    }

    /// `from_sigma` accepts a standard deviation and produces that σ.
    #[test]
    fn from_sigma_takes_a_standard_deviation() {
        let sigma_t_us = 1.0;
        let from_sigma = ResolutionParams::from_sigma(25.0, sigma_t_us, 0.0, 0.0).expect("valid");
        let direct = ResolutionParams::new(25.0, sigma_t_us, 0.0, 0.0).expect("valid");

        // The converted parameters are √2 wider than the raw ones.
        assert!(
            (from_sigma.delta_t_us() - sigma_t_us * SQRT_2).abs() < 1.0e-15,
            "from_sigma did not apply W = σ·√2"
        );

        // And the kernel it builds has the σ the caller asked for, in energy:
        // σ_E = 2·σ_t·E^{3/2}/(TOF_FACTOR·L), the ordinary propagation.
        let center = 10.0_f64;
        let expected_sigma_e =
            2.0 * sigma_t_us * center.powf(1.5) / (TOF_FACTOR * from_sigma.flight_path_m());
        let measured = measured_sigma_ev(&from_sigma, center, from_sigma.gaussian_width(center));
        assert!(
            (measured - expected_sigma_e).abs() / expected_sigma_e < 2.0e-3,
            "from_sigma kernel measures σ {measured:.9}, asked for {expected_sigma_e:.9}"
        );
        // The direct constructor, given the same number, is √2 narrower —
        // which is exactly the error the documentation used to invite.
        let measured_direct = measured_sigma_ev(&direct, center, direct.gaussian_width(center));
        assert!(
            (measured_direct * SQRT_2 - measured).abs() / measured < 5.0e-3,
            "the two constructors do not differ by √2"
        );
    }

    /// `from_fwhm` accepts a full width at half maximum.
    #[test]
    fn from_fwhm_takes_a_full_width_at_half_maximum() {
        let fwhm_t_us = 1.0;
        let params = ResolutionParams::from_fwhm(25.0, fwhm_t_us, 0.0, 0.0).expect("valid");
        assert!(
            (params.delta_t_us() - fwhm_t_us / FWHM_PER_W).abs() < 1.0e-15,
            "from_fwhm did not apply W = FWHM/(2√ln2)"
        );
    }

    /// `from_fwhm` is the conversion SAMMY's `Deltag` needs.
    ///
    /// `nereids_endf::sammy::sammy_to_nereids_resolution` divides Deltag by
    /// 2√(ln 2) to reach this module's convention; the two must agree, because
    /// the samtry baselines depend on that mapping being the right one.
    #[test]
    fn from_fwhm_agrees_with_the_sammy_deltag_conversion() {
        let delta_g = 0.022_f64; // tr007's BROADENING card value.
        let via_constructor = ResolutionParams::from_fwhm(25.0, delta_g, 0.0, 0.0).expect("valid");
        let sammy_mapping = delta_g / (2.0 * 2.0_f64.ln().sqrt());
        assert!(
            (via_constructor.delta_t_us() - sammy_mapping).abs() < 1.0e-15,
            "from_fwhm disagrees with the SAMMY Deltag conversion"
        );
    }
}

#[cfg(test)]
mod flight_path_rebinding_tests {
    use super::*;

    /// Rebinding the flight path scales the Gaussian energy width by 1/s.
    ///
    /// The oracle is the analytic law, not the code: a timing width Δt maps to
    /// an energy width `W_E = 2·Δt·E^{3/2}/(F·L)`, so the same Δt read against
    /// `L·s` gives `W_E/s`. A fit that frees `L_scale` builds its grid with
    /// `L·L_scale`; if the kernel keeps `L`, this is exactly the factor it is
    /// wrong by.
    #[test]
    fn rebinding_scales_the_gaussian_width_inversely() {
        let l_nom = 25.0;
        let delta_t = 1.0;
        let base = ResolutionParams::new(l_nom, delta_t, 0.0, 0.0).expect("valid");
        for s in [0.97_f64, 1.0, 1.03, 1.25] {
            let rebound = ResolutionFunction::Gaussian(base)
                .with_flight_path(l_nom * s)
                .expect("positive flight path");
            let ResolutionFunction::Gaussian(rebound) = rebound else {
                panic!("rebinding changed the family");
            };
            for e in [5.0_f64, 20.0, 100.0] {
                let analytic = 2.0 * delta_t * e.powf(1.5) / (TOF_FACTOR * l_nom * s);
                let got = rebound.gaussian_width(e);
                assert!(
                    (got - analytic).abs() / analytic < 1.0e-14,
                    "s={s}, E={e}: width {got:e} but the law gives {analytic:e}"
                );
                // And it is the base width divided by s, which is the
                // statement the fit depends on.
                let base_w = base.gaussian_width(e);
                assert!(
                    (got - base_w / s).abs() / (base_w / s) < 1.0e-14,
                    "s={s}, E={e}: rebound width is not base/s"
                );
            }
        }
    }

    /// Rebinding does not touch the tabulated kernel itself.
    ///
    /// The stored offsets are emission times. The flight path belongs to the
    /// map they are applied through, so a rebind must leave every offset and
    /// weight byte-identical — if it resynthesized or rescaled them it would
    /// be changing the moderator, not the geometry.
    #[test]
    fn rebinding_leaves_the_tabulated_kernel_byte_identical() {
        let offsets = vec![-2.0, -1.0, 0.0, 1.0, 3.0];
        let weights = vec![0.1, 0.6, 1.0, 0.5, 0.05];
        let base = TabulatedResolution::from_kernels(
            vec![5.0, 50.0],
            vec![
                (offsets.clone(), weights.clone()),
                (offsets.clone(), weights.clone()),
            ],
            25.0,
        )
        .expect("valid kernel table");

        let rebound = base.with_flight_path(25.0 * 1.03).expect("positive");
        assert_eq!(rebound.flight_path_m(), 25.0 * 1.03);
        assert_eq!(base.ref_energies, rebound.ref_energies);
        for (i, ((b_off, b_w), (r_off, r_w))) in
            base.kernels.iter().zip(rebound.kernels.iter()).enumerate()
        {
            assert_eq!(
                b_off.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
                r_off.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
                "rebinding moved the emission-time offsets of block {i}"
            );
            assert_eq!(
                b_w.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
                r_w.iter().map(|v| v.to_bits()).collect::<Vec<_>>(),
                "rebinding changed the weights of block {i}"
            );
        }
    }

    /// A non-positive flight path is refused rather than silently accepted.
    #[test]
    fn rebinding_refuses_a_non_physical_flight_path() {
        let params = ResolutionParams::new(25.0, 1.0, 0.0, 0.0).expect("valid");
        for bad in [0.0, -1.0, f64::NAN, f64::INFINITY] {
            assert!(
                ResolutionFunction::Gaussian(params)
                    .with_flight_path(bad)
                    .is_err(),
                "flight path {bad} was accepted"
            );
        }
    }
}