nereids_physics/

reich_moore.rs

//! Multi-formalism cross-section dispatcher.
//!
//! `cross_sections_at_energy` is the primary entry point for computing
//! energy-dependent cross-sections from ENDF resonance data.  It
//! iterates over all resonance ranges in the data and dispatches each to
//! the appropriate formalism-specific calculator:
//!
//! All formalisms route through a single dispatch pipeline —
//! [`cross_sections_on_grid`] precomputes per-range invariants once, then
//! `evaluate_precomputed_range` evaluates the correct formalism at each
//! energy.  [`cross_sections_at_energy`] is a one-call convenience
//! wrapper around the same pipeline.  There is exactly one entry point
//! and one evaluator per formalism:
//!
//! | ENDF LRF | Formalism                  | Evaluator                                         |
//! |----------|----------------------------|---------------------------------------------------|
//! | 1        | SLBW                       | `slbw::slbw_evaluate_with_cached_jgroups`         |
//! | 2        | MLBW                       | `slbw::mlbw_evaluate_with_cached_jgroups`         |
//! | 3        | Reich-Moore                | `reich_moore_spin_group_precomputed` (+ 2ch/3ch)  |
//!
//! ## Reich-Moore Approximation
//! In the full R-matrix, all channels (neutron, capture, fission) appear
//! explicitly. The Reich-Moore approximation *eliminates the capture channel*
//! from the channel space, absorbing its effect into an imaginary part of
//! the energy denominator. This makes the level matrix smaller while
//! remaining highly accurate.
//!
//! For non-fissile isotopes (like U-238 below threshold), each spin group
//! has only ONE explicit channel (neutron elastic), making the R-matrix
//! a scalar — and the calculation is very efficient.
//!
//! ## SAMMY Reference
//! - `rml/mrml07.f` Setr: R-matrix construction
//! - `rml/mrml09.f` Yinvrs: level matrix inversion
//! - `rml/mrml11.f` Setxqx: X-matrix, Sectio: cross-sections
//! - `rml/mrml03.f` Betset: ENDF widths → reduced width amplitudes
//! - SAMMY manual Section II.B.1 (Reich-Moore approximation)

use num_complex::Complex64;

use nereids_core::constants::{DIVISION_FLOOR, LOG_FLOOR, PIVOT_FLOOR, QUANTUM_NUMBER_EPS};
use nereids_endf::resonance::{ResonanceData, ResonanceFormalism, ResonanceRange, Tab1};

use crate::channel;
use crate::penetrability;
use crate::slbw;

// ─── Per-resonance precomputed invariants ─────────────────────────────────────
//
// These quantities depend only on the resonance parameters and the channel
// radius at the resonance energy — both are energy-independent constants.
// Pre-computing them once (outside the energy loop) eliminates redundant
// `penetrability(l, rho_r)` and `group_by_j()` calls per energy point.
//
// Issue #87: "Perf: Pre-cache J-groups and per-resonance quantities"

/// Per-resonance invariants for the single-channel (non-fissile) Reich-Moore path.
///
/// Pre-computed once per resonance before the energy sweep.
/// Reference: SAMMY `rml/mrml03.f` Betset (lines 240-276)
struct PrecomputedResonanceSingle {
    /// Resonance energy E_r (eV).
    energy: f64,
    /// Capture width Γ_γ (eV).
    gamma_g: f64,
    /// Reduced width amplitude squared γ²_n = |Γ_n| / (2·P_l(E_r)).
    gamma_n_reduced_sq: f64,
}

/// Per-resonance invariants for the 2-channel (one fission) Reich-Moore path.
struct PrecomputedResonance2ch {
    /// Resonance energy E_r (eV).
    energy: f64,
    /// Capture width Γ_γ (eV).
    gamma_g: f64,
    /// Reduced width amplitude β_n = sign(Γ_n) × √(|Γ_n| / (2·P_l(E_r))).
    beta_n: f64,
    /// Fission width amplitude β_f = sign(Γ_f) × √(|Γ_f| / 2).
    beta_f: f64,
}

/// Per-resonance invariants for the 3-channel (two fission) Reich-Moore path.
struct PrecomputedResonance3ch {
    /// Resonance energy E_r (eV).
    energy: f64,
    /// Capture width Γ_γ (eV).
    gamma_g: f64,
    /// Reduced width amplitude β_n.
    beta_n: f64,
    /// Fission width amplitude β_fa.
    beta_fa: f64,
    /// Fission width amplitude β_fb.
    beta_fb: f64,
}

/// Pre-computed J-group, generic over the per-resonance invariant type.
///
/// Groups resonances by total angular momentum J, with per-resonance
/// invariants already computed. The J-grouping depends only on the
/// resonance data, not on the incident energy, so it is computed once
/// and reused for every energy point.
///
/// Used by both Reich-Moore (single/2ch/3ch) and SLBW precompute paths.
pub(crate) struct PrecomputedJGroup<R> {
    /// Total angular momentum J (signed, per SAMMY convention).
    pub(crate) j: f64,
    /// Statistical weight g_J = (2J+1) / ((2I+1)(2s+1)).
    pub(crate) g_j: f64,
    /// Pre-computed per-resonance quantities for this J-group.
    pub(crate) resonances: Vec<R>,
}

/// Type aliases for each channel count (preserves readability at call sites).
type PrecomputedJGroupSingle = PrecomputedJGroup<PrecomputedResonanceSingle>;
type PrecomputedJGroup2ch = PrecomputedJGroup<PrecomputedResonance2ch>;
type PrecomputedJGroup3ch = PrecomputedJGroup<PrecomputedResonance3ch>;

/// Compute penetrability at the resonance energy P_l(ρ_r).
///
/// ENDF widths are defined as Γ_n = 2·P_l(AP(E_r), E_r)·γ²_n,
/// so the penetrability must be evaluated at the resonance energy
/// using the channel radius AP(E_r) — not the incident-energy AP(E).
///
/// This function is the core quantity that Issue #87 caches: previously
/// it was recomputed for every resonance at every energy point.
///
/// When E_r ≈ 0, the penetrability is zero (matching SLBW behavior in
/// `slbw.rs`). This ensures the result depends only on resonance
/// parameters and is independent of the incident energy, enabling the
/// precompute to be hoisted above the energy loop.
fn penetrability_at_resonance(
    e_r: f64,
    l: u32,
    awr: f64,
    channel_radius: f64,
    ap_table: Option<&Tab1>,
) -> f64 {
    if e_r.abs() > PIVOT_FLOOR {
        let radius_at_er = ap_table.map_or(channel_radius, |t| t.evaluate(e_r.abs()));
        let rho_r = channel::rho(e_r.abs(), awr, radius_at_er);
        penetrability::penetrability(l, rho_r)
    } else {
        // E_r ≈ 0 → P_l(0) = 0, so γ²_n = 0 regardless.
        // Using 0.0 keeps this function energy-independent.
        0.0
    }
}

/// Group resonances by total angular momentum J, building per-resonance
/// precomputed invariants via a caller-supplied closure.
///
/// This is the shared core of all `precompute_jgroups_*` functions (RM single,
/// 2ch, 3ch) and SLBW's `precompute_slbw_jgroups`.  The closure `build_resonance`
/// receives each ENDF resonance and returns the per-resonance struct `R` that
/// differs between formalisms and channel counts.
///
/// The grouping logic is identical across all callers:
/// 1. Extract J from each resonance.
/// 2. Find or create a J-group (matching within `QUANTUM_NUMBER_EPS`).
/// 3. Compute `g_J = (2J+1) / ((2I+1)(2s+1))` for new groups.
/// 4. Push the precomputed resonance into the matching group.
pub(crate) fn group_resonances_by_j<R>(
    resonances: &[nereids_endf::resonance::Resonance],
    target_spin: f64,
    mut build_resonance: impl FnMut(&nereids_endf::resonance::Resonance) -> R,
) -> Vec<PrecomputedJGroup<R>> {
    let mut j_values: Vec<f64> = Vec::new();
    let mut groups: Vec<PrecomputedJGroup<R>> = Vec::new();

    for res in resonances {
        let j = res.j;
        let precomp = build_resonance(res);

        if let Some(idx) = j_values
            .iter()
            .position(|&gj| (gj - j).abs() < QUANTUM_NUMBER_EPS)
        {
            groups[idx].resonances.push(precomp);
        } else {
            j_values.push(j);
            groups.push(PrecomputedJGroup {
                j,
                g_j: channel::statistical_weight(j, target_spin),
                resonances: vec![precomp],
            });
        }
    }
    groups
}

/// Build pre-computed J-groups for the single-channel (non-fissile) path.
///
/// Groups resonances by J, pre-computes γ²_n per resonance.
/// All quantities depend only on resonance parameters (not incident energy),
/// so the result can be computed once and reused across all energy points.
fn precompute_jgroups_single(
    resonances: &[nereids_endf::resonance::Resonance],
    l: u32,
    awr: f64,
    channel_radius: f64,
    ap_table: Option<&Tab1>,
    target_spin: f64,
) -> Vec<PrecomputedJGroupSingle> {
    group_resonances_by_j(resonances, target_spin, |res| {
        let p_at_er = penetrability_at_resonance(res.energy, l, awr, channel_radius, ap_table);
        let gamma_n_reduced_sq = if p_at_er > PIVOT_FLOOR {
            res.gn.abs() / (2.0 * p_at_er)
        } else {
            0.0
        };
        PrecomputedResonanceSingle {
            energy: res.energy,
            gamma_g: res.gg,
            gamma_n_reduced_sq,
        }
    })
}

/// Build pre-computed J-groups for the 2-channel fission path.
///
/// All quantities depend only on resonance parameters (not incident energy),
/// so the result can be computed once and reused across all energy points.
fn precompute_jgroups_2ch(
    resonances: &[nereids_endf::resonance::Resonance],
    l: u32,
    awr: f64,
    channel_radius: f64,
    ap_table: Option<&Tab1>,
    target_spin: f64,
) -> Vec<PrecomputedJGroup2ch> {
    group_resonances_by_j(resonances, target_spin, |res| {
        let p_at_er = penetrability_at_resonance(res.energy, l, awr, channel_radius, ap_table);

        let beta_n = if p_at_er > PIVOT_FLOOR {
            let sign = if res.gn >= 0.0 { 1.0 } else { -1.0 };
            sign * (res.gn.abs() / (2.0 * p_at_er)).sqrt()
        } else {
            0.0
        };

        let beta_f = {
            let sign = if res.gfa >= 0.0 { 1.0 } else { -1.0 };
            sign * (res.gfa.abs() / 2.0).sqrt()
        };

        PrecomputedResonance2ch {
            energy: res.energy,
            gamma_g: res.gg,
            beta_n,
            beta_f,
        }
    })
}

/// Build pre-computed J-groups for the 3-channel fission path.
///
/// All quantities depend only on resonance parameters (not incident energy),
/// so the result can be computed once and reused across all energy points.
fn precompute_jgroups_3ch(
    resonances: &[nereids_endf::resonance::Resonance],
    l: u32,
    awr: f64,
    channel_radius: f64,
    ap_table: Option<&Tab1>,
    target_spin: f64,
) -> Vec<PrecomputedJGroup3ch> {
    group_resonances_by_j(resonances, target_spin, |res| {
        let p_at_er = penetrability_at_resonance(res.energy, l, awr, channel_radius, ap_table);

        let beta_n = if p_at_er > PIVOT_FLOOR {
            let sign = if res.gn >= 0.0 { 1.0 } else { -1.0 };
            sign * (res.gn.abs() / (2.0 * p_at_er)).sqrt()
        } else {
            0.0
        };

        let beta_fa = {
            let sign = if res.gfa >= 0.0 { 1.0 } else { -1.0 };
            sign * (res.gfa.abs() / 2.0).sqrt()
        };

        let beta_fb = {
            let sign = if res.gfb >= 0.0 { 1.0 } else { -1.0 };
            sign * (res.gfb.abs() / 2.0).sqrt()
        };

        PrecomputedResonance3ch {
            energy: res.energy,
            gamma_g: res.gg,
            beta_n,
            beta_fa,
            beta_fb,
        }
    })
}

/// Cross-section results at a single energy point.
#[derive(Debug, Clone, Copy)]
pub struct CrossSections {
    /// Total cross-section (barns).
    pub total: f64,
    /// Elastic scattering cross-section (barns).
    pub elastic: f64,
    /// Capture (n,γ) cross-section (barns).
    pub capture: f64,
    /// Fission cross-section (barns).
    pub fission: f64,
}

/// Compute cross-sections at a single energy.
///
/// Dispatches each resonance range to the appropriate formalism-specific
/// calculator (SLBW, MLBW, Reich-Moore) based on the formalism stored in
/// that range; non-evaluable ranges (LRF=7, LRU=2) resolve to `Skip` and
/// contribute zero.  See the module-level table for the full dispatch map.
///
/// Adjacent ranges that share a boundary energy use half-open intervals
/// `[e_low, e_high)` so the boundary point is counted exactly once
/// (ENDF-6 §2 convention).
///
/// Shares the **same precompute+evaluate pipeline** as
/// [`cross_sections_on_grid`] — both entry points call the same
/// (private) `precompute_range_data` and `evaluate_precomputed_range`
/// helpers, so there is exactly one dispatch table per formalism in
/// the codebase.  The difference is that this entry point does not
/// store precomputed plans in an outer `Vec` (unnecessary when only
/// one energy is evaluated), and it skips the precompute entirely
/// for ranges whose energy interval excludes `energy_ev`.  Per-call
/// overhead is therefore close to — though not exactly — the
/// pre-consolidation per-point path; a small residual cost remains
/// because each matching range is wrapped in a `PrecomputedRangeData`
/// before evaluation.  Measured A.1 LM+grouped walltime on real
/// VENUS Hf 120 min data: 1.37 s (pre-consolidation) → 1.42 s (this
/// path), ≈ 3.6 % overhead.
///
/// # Arguments
/// * `data` — Parsed resonance parameters from ENDF.
/// * `energy_ev` — Neutron energy in eV (lab frame).
///
/// # Returns
/// Cross-sections in barns.
///
/// # Panics
/// Panics if `energy_ev` is non-finite or non-positive.  The SLBW leaf
/// routine already enforces this precondition in release builds; hoisting
/// the same assert to the top-level pub fn keeps the public contract
/// symmetric so direct Rust callers cannot bypass validation by hitting a
/// range that gates entry on a finite-only check (e.g. a pure-RM range with
/// no SLBW leaf would otherwise silently return zeros when handed NaN).
pub fn cross_sections_at_energy(data: &ResonanceData, energy_ev: f64) -> CrossSections {
    // Symmetric public-API guard.  Matches `slbw_cross_sections_for_range`.
    // One branch at the entry of this O(ranges × resonances) function is
    // negligible.
    assert!(
        energy_ev.is_finite() && energy_ev > 0.0,
        "expected positive finite energy_ev, got {energy_ev}"
    );

    let awr = data.awr;

    let mut total = 0.0;
    let mut elastic = 0.0;
    let mut capture = 0.0;
    let mut fission = 0.0;

    for (range_idx, range) in data.ranges.iter().enumerate() {
        // Cheap interval check FIRST — precompute is O(n_resonances) per range,
        // so building a plan for a range that doesn't cover `energy_ev` wastes
        // meaningful work on multi-range isotopes or energies outside the
        // resolved band.
        if !covers(
            range.energy_low,
            range.energy_high,
            upper_bound_is_half_open(data, range_idx),
            energy_ev,
        ) {
            continue;
        }

        let plan = precompute_range_data(range, range_idx, data, awr);
        let (t, e, c, f) = evaluate_precomputed_range(&plan, energy_ev, awr);
        total += t;
        elastic += e;
        capture += c;
        fission += f;
    }

    CrossSections {
        total,
        elastic,
        capture,
        fission,
    }
}

/// Compute cross-sections over a grid of energies.
///
/// Optimized batch evaluation: precomputes J-groups and per-resonance
/// invariants (reduced width amplitudes, penetrability at E_r) once per
/// resonance range, then evaluates each energy point using the cached data.
/// This avoids redundant `group_by_j` + `penetrability(l, rho_r)` calls
/// that the per-point API (`cross_sections_at_energy`) would repeat.
///
/// Issue #87: the precompute is hoisted above the energy loop so that
/// `precompute_jgroups_*` runs O(ranges) times total, not O(ranges × energies).
///
/// # Arguments
/// * `data` — Parsed resonance parameters from ENDF.
/// * `energies` — Slice of neutron energies in eV.
///
/// # Returns
/// Vector of cross-sections, one per energy point.
///
/// # Panics
/// Panics if any element of `energies` is non-finite or non-positive.
/// Validating the entire grid up-front (O(n) branch, one pass) means a
/// single bad energy fails fast with a clear message instead of being
/// hidden inside the inner loop, matches the symmetric contract on
/// `cross_sections_at_energy`, and protects direct Rust callers from
/// the same release-mode silent-zero footgun that the SLBW / Reich-Moore
/// leaf asserts guard against per-point.
pub fn cross_sections_on_grid(data: &ResonanceData, energies: &[f64]) -> Vec<CrossSections> {
    if energies.is_empty() {
        return Vec::new();
    }
    CrossSectionPlan::new(data).evaluate(energies)
}

/// Whether range `range_idx` excludes its upper bound.
///
/// When the next range starts exactly where this one ends and is
/// evaluable, the shared boundary energy belongs to the next range only,
/// so a point on the boundary is never evaluated with both formalisms.
/// The per-point dispatcher, the plan and the Doppler route gate all use
/// this one convention. [`crate::slbw::slbw_cross_sections`] deliberately
/// does not: it closes both ends so that a point on a resolved/unresolved
/// boundary is not dropped by both ranges, which its own comment explains.
/// Treat that as the documented exception rather than a site left behind.
pub(crate) fn upper_bound_is_half_open(data: &ResonanceData, range_idx: usize) -> bool {
    let range = &data.ranges[range_idx];
    data.ranges
        .get(range_idx + 1)
        .is_some_and(|next| next.energy_low == range.energy_high && range_is_evaluable(next))
}

/// Whether `[low, high]` covers `energy_ev` under the given upper-bound
/// convention. The one definition of range coverage: the per-point
/// dispatcher, the plan and the Doppler route gate all call it, so they
/// cannot drift apart on a boundary energy.
pub(crate) fn covers(low: f64, high: f64, half_open_upper: bool, energy_ev: f64) -> bool {
    if half_open_upper {
        energy_ev >= low && energy_ev < high
    } else {
        energy_ev >= low && energy_ev <= high
    }
}

/// Reusable evaluation plan for one immutable resonance source.
///
/// The formalism-specific resonance caches (J-groups, reduced width
/// amplitudes, penetrability at each resonance energy) are built once by
/// [`CrossSectionPlan::new`] and reused by every evaluation. This is the
/// primitive that lets the continuous Doppler integral evaluate the same
/// ENDF equation on many adaptively chosen energies without rebuilding the
/// caches per quadrature panel; [`cross_sections_on_grid`] is the plan
/// applied to one grid.
pub struct CrossSectionPlan<'a> {
    awr: f64,
    precomputed: Vec<PrecomputedRangeData<'a>>,
}

impl<'a> CrossSectionPlan<'a> {
    /// Build the per-range caches for `data`.
    pub fn new(data: &'a ResonanceData) -> Self {
        let awr = data.awr;
        let precomputed = data
            .ranges
            .iter()
            .enumerate()
            .map(|(range_idx, range)| precompute_range_data(range, range_idx, data, awr))
            .collect();
        Self { awr, precomputed }
    }

    /// Cross-sections at one energy, summed over every range covering it.
    ///
    /// # Panics
    /// Panics if `energy_ev` is non-finite or non-positive — the same
    /// public contract as [`cross_sections_at_energy`].
    pub fn evaluate_one(&self, energy_ev: f64) -> CrossSections {
        assert!(
            energy_ev.is_finite() && energy_ev > 0.0,
            "expected positive finite energy_ev, got {energy_ev}"
        );

        let mut total = 0.0;
        let mut elastic = 0.0;
        let mut capture = 0.0;
        let mut fission = 0.0;

        for pc in &self.precomputed {
            if !covers(pc.energy_low, pc.energy_high, pc.half_open_upper, energy_ev) {
                continue;
            }

            let (t, e, c, f) = evaluate_precomputed_range(pc, energy_ev, self.awr);
            total += t;
            elastic += e;
            capture += c;
            fission += f;
        }

        CrossSections {
            total,
            elastic,
            capture,
            fission,
        }
    }

    /// Cross-sections at every energy of `energies`, in order.
    ///
    /// # Panics
    /// Panics if any element of `energies` is non-finite or non-positive.
    /// Validating the entire grid up-front (O(n) branch, one pass) means a
    /// single bad energy fails fast with a clear message instead of being
    /// hidden inside the inner loop.
    pub fn evaluate(&self, energies: &[f64]) -> Vec<CrossSections> {
        for &energy_ev in energies {
            assert!(
                energy_ev.is_finite() && energy_ev > 0.0,
                "expected positive finite energy_ev, got {energy_ev}"
            );
        }
        energies
            .iter()
            .map(|&energy_ev| self.evaluate_one(energy_ev))
            .collect()
    }
}

// ─── Precomputed range data for batch grid evaluation ────────────────────────
//
// These types hold energy-independent invariants for a single resonance range,
// precomputed once by `precompute_range_data` and reused for every energy point
// in `cross_sections_on_grid`.

/// Precomputed data for a single L-group within a Reich-Moore range.
///
/// Holds the J-group cache and metadata needed to compute energy-dependent
/// channel parameters (rho, P_l, S_l, phi_l) at each energy point.
enum PrecomputedRmLGroupData {
    /// Non-fissile: single neutron channel, capture eliminated.
    Single {
        l: u32,
        awr_l: f64,
        /// L-group override radius (fm). Used only for scattering radius
        /// (phase shifts). 0.0 means use range radius.
        apl: f64,
        /// Precomputed penetrability radius (fm): APL when set, else
        /// NAPS=0 formula. 0.0 means fall back to `scatt_radius`.
        pen_radius_override: f64,
        jgroups: Vec<PrecomputedJGroupSingle>,
    },
    /// One fission channel (gfa != 0, gfb == 0).
    TwoCh {
        l: u32,
        awr_l: f64,
        apl: f64,
        pen_radius_override: f64,
        jgroups: Vec<PrecomputedJGroup2ch>,
    },
    /// Two fission channels (both gfa and gfb != 0).
    ThreeCh {
        l: u32,
        awr_l: f64,
        apl: f64,
        pen_radius_override: f64,
        jgroups: Vec<PrecomputedJGroup3ch>,
    },
}

/// Precomputed data for a single SLBW L-group.
struct PrecomputedSlbwLGroupData {
    l: u32,
    awr_l: f64,
    /// L-group override radius (fm). 0.0 means use range radius.
    apl: f64,
    /// Precomputed penetrability radius (fm): APL when set, else
    /// NAPS=0 formula. 0.0 means fall back to `scatt_radius`.
    pen_radius_override: f64,
    jgroups: Vec<slbw::PrecomputedSlbwJGroup>,
}

/// Precomputed data for a single resonance range.
///
/// Wraps the formalism-specific precomputed L-group data plus the
/// energy interval metadata needed for range dispatch.
struct PrecomputedRangeData<'a> {
    energy_low: f64,
    energy_high: f64,
    half_open_upper: bool,
    kind: PrecomputedRangeKind<'a>,
}

/// Formalism-specific precomputed data for a range.
///
/// `Slbw` and `Mlbw` share the same precomputed J-group layout but
/// dispatch to different evaluators — SLBW's incoherent per-resonance
/// elastic sum vs MLBW's coherent-sum elastic.  Keeping them as
/// distinct variants (instead of a single variant with a formalism
/// tag) is deliberate: it makes the "pick the right evaluator" step a
/// `match` arm that the compiler checks exhaustively, which is what
/// prevents the #465 class of bug (MLBW being silently routed through
/// the SLBW evaluator) from reappearing.
enum PrecomputedRangeKind<'a> {
    /// Reich-Moore (LRF=3): precomputed J-groups per L-group.
    /// The range reference is kept for `scattering_radius_at(energy_ev)`.
    ReichMoore {
        range: &'a ResonanceRange,
        l_groups: Vec<PrecomputedRmLGroupData>,
    },
    /// Single-Level Breit-Wigner (LRF=1): incoherent per-resonance sums.
    /// The range reference is kept for `scattering_radius_at(energy_ev)`.
    Slbw {
        range: &'a ResonanceRange,
        l_groups: Vec<PrecomputedSlbwLGroupData>,
    },
    /// Multi-Level Breit-Wigner (LRF=2): coherent-sum elastic, same
    /// capture/fission as SLBW.  Uses the same precomputed-J-group
    /// layout as SLBW but dispatches to `mlbw_evaluate_with_cached_jgroups`
    /// (see issue #465 for why this MUST be a distinct variant).
    Mlbw {
        range: &'a ResonanceRange,
        l_groups: Vec<PrecomputedSlbwLGroupData>,
    },
    /// Not evaluable (skip).
    Skip,
}

/// Build precomputed range data for a single resonance range.
///
/// This extracts all energy-independent quantities (J-group structure,
/// reduced width amplitudes, penetrability at resonance energies) so they
/// can be reused across all energy points without redundant computation.
fn precompute_range_data<'a>(
    range: &'a ResonanceRange,
    range_idx: usize,
    data: &'a ResonanceData,
    awr: f64,
) -> PrecomputedRangeData<'a> {
    let make = |kind| PrecomputedRangeData {
        energy_low: range.energy_low,
        energy_high: range.energy_high,
        half_open_upper: upper_bound_is_half_open(data, range_idx),
        kind,
    };

    // Literal sync with `ResonanceRange::is_evaluable` (see
    // `range_is_evaluable` below): every non-evaluable shape — parse-and-skip
    // placeholders AND accepted-but-inert resolved ranges whose L-groups are
    // all empty — takes the Skip arm. The inert-resolved shape is unreachable
    // from the parser and the Python constructor (both guard it) but can
    // arrive via legacy serialized caches or direct Rust construction.
    if !range.is_evaluable() {
        return make(PrecomputedRangeKind::Skip);
    }

    // SLBW and MLBW share the precomputed-J-group layout but evaluate
    // with different math (see `PrecomputedRangeKind` doc).  Build the
    // shared precomputed data first, then wrap in the formalism-specific
    // variant so the evaluator dispatch is exhaustive.
    if matches!(
        range.formalism,
        ResonanceFormalism::SLBW | ResonanceFormalism::MLBW
    ) {
        let l_groups: Vec<PrecomputedSlbwLGroupData> = range
            .l_groups
            .iter()
            .map(|l_group| {
                let l = l_group.l;
                let awr_l = if l_group.awr > 0.0 { l_group.awr } else { awr };
                let jgroups = slbw::precompute_slbw_jgroups(
                    &l_group.resonances,
                    l,
                    awr_l,
                    range,
                    l_group,
                    range.target_spin,
                );
                let pen_radius_override = if l_group.apl > 0.0 {
                    l_group.apl
                } else if range.naps == 0 {
                    channel::endf_channel_radius_fm(awr_l)
                } else {
                    0.0
                };
                PrecomputedSlbwLGroupData {
                    l,
                    awr_l,
                    apl: l_group.apl,
                    pen_radius_override,
                    jgroups,
                }
            })
            .collect();
        return match range.formalism {
            ResonanceFormalism::SLBW => make(PrecomputedRangeKind::Slbw { range, l_groups }),
            ResonanceFormalism::MLBW => make(PrecomputedRangeKind::Mlbw { range, l_groups }),
            // Unreachable: the outer `matches!` already restricted to SLBW|MLBW.
            _ => unreachable!("formalism guard admits only SLBW/MLBW"),
        };
    }

    // Reich-Moore ranges: precompute J-groups per L-group.
    if range.formalism == ResonanceFormalism::ReichMoore {
        let rm_l_groups: Vec<PrecomputedRmLGroupData> = range
            .l_groups
            .iter()
            .map(|l_group| {
                let l = l_group.l;
                let awr_l = if l_group.awr > 0.0 { l_group.awr } else { awr };

                // Pen-radius override: energy-independent when APL > 0 or NAPS=0.
                // Stored in the precomputed struct to avoid recomputation per energy.
                let pen_radius_override = if l_group.apl > 0.0 {
                    l_group.apl
                } else if range.naps == 0 {
                    channel::endf_channel_radius_fm(awr_l)
                } else {
                    0.0
                };
                // Channel radius for precompute: when a penetrability radius
                // override is available (APL > 0 or NAPS=0), use that
                // precomputed radius; otherwise fall back to the constant
                // scattering_radius. The ap_table (NRO=1) case is handled
                // inside penetrability_at_resonance, which evaluates the
                // table at E_r for each resonance.
                let channel_radius = if pen_radius_override > 0.0 {
                    pen_radius_override
                } else {
                    range.scattering_radius
                };
                // NAPS=0: penetrability uses the formula radius, not the AP(E) table.
                let ap_table_ref: Option<&Tab1> = if l_group.apl > 0.0 || range.naps == 0 {
                    None
                } else {
                    range.ap_table.as_ref()
                };

                let has_fission = l_group
                    .resonances
                    .iter()
                    .any(|r| r.gfa.abs() > PIVOT_FLOOR || r.gfb.abs() > PIVOT_FLOOR);
                let has_two_fission = l_group.resonances.iter().any(|r| r.gfb.abs() > PIVOT_FLOOR);

                if !has_fission {
                    let jgroups = precompute_jgroups_single(
                        &l_group.resonances,
                        l,
                        awr_l,
                        channel_radius,
                        ap_table_ref,
                        range.target_spin,
                    );
                    PrecomputedRmLGroupData::Single {
                        l,
                        awr_l,
                        apl: l_group.apl,
                        pen_radius_override,
                        jgroups,
                    }
                } else if !has_two_fission {
                    // P-7: R-external now applied in the 2ch evaluation path.
                    let jgroups = precompute_jgroups_2ch(
                        &l_group.resonances,
                        l,
                        awr_l,
                        channel_radius,
                        ap_table_ref,
                        range.target_spin,
                    );
                    PrecomputedRmLGroupData::TwoCh {
                        l,
                        awr_l,
                        apl: l_group.apl,
                        pen_radius_override,
                        jgroups,
                    }
                } else {
                    // P-7: R-external now applied in the 3ch evaluation path.
                    let jgroups = precompute_jgroups_3ch(
                        &l_group.resonances,
                        l,
                        awr_l,
                        channel_radius,
                        ap_table_ref,
                        range.target_spin,
                    );
                    PrecomputedRmLGroupData::ThreeCh {
                        l,
                        awr_l,
                        apl: l_group.apl,
                        pen_radius_override,
                        jgroups,
                    }
                }
            })
            .collect();
        return make(PrecomputedRangeKind::ReichMoore {
            range,
            l_groups: rm_l_groups,
        });
    }

    // Unrecognized formalism: skip.
    make(PrecomputedRangeKind::Skip)
}

/// Evaluate cross-sections for a precomputed range at a single energy.
///
/// Uses the cached J-groups and per-resonance invariants to avoid
/// redundant precomputation. Only energy-dependent quantities (rho,
/// P_l, S_l, phi_l, pi/k^2) are computed per call.
fn evaluate_precomputed_range(
    pc: &PrecomputedRangeData,
    energy_ev: f64,
    awr: f64,
) -> (f64, f64, f64, f64) {
    match &pc.kind {
        PrecomputedRangeKind::Skip => (0.0, 0.0, 0.0, 0.0),

        PrecomputedRangeKind::Slbw { range, l_groups } => {
            let pi_over_k2 = channel::pi_over_k_squared_barns(energy_ev, awr);
            let mut total = 0.0;
            let mut elastic = 0.0;
            let mut capture = 0.0;
            let mut fission = 0.0;

            for lg in l_groups {
                // Scattering radius for phase shift (always AP/APL).
                let scatt_radius = if lg.apl > 0.0 {
                    lg.apl
                } else {
                    range.scattering_radius_at(energy_ev)
                };
                // Penetrability radius: precomputed override (APL or NAPS=0
                // formula), falling back to scattering radius.
                let pen_radius = if lg.pen_radius_override > 0.0 {
                    lg.pen_radius_override
                } else {
                    scatt_radius
                };

                let rho_phase = channel::rho(energy_ev, lg.awr_l, scatt_radius);
                let rho_pen = channel::rho(energy_ev, lg.awr_l, pen_radius);
                let phi = penetrability::phase_shift(lg.l, rho_phase);
                let sin_phi = phi.sin();
                let cos_phi = phi.cos();
                let sin2_phi = sin_phi * sin_phi;
                let p_at_e = penetrability::penetrability(lg.l, rho_pen);

                let (t, e, c, f) = slbw::slbw_evaluate_with_cached_jgroups(
                    &lg.jgroups,
                    energy_ev,
                    pi_over_k2,
                    p_at_e,
                    sin_phi,
                    cos_phi,
                    sin2_phi,
                );
                total += t;
                elastic += e;
                capture += c;
                fission += f;
            }

            (total, elastic, capture, fission)
        }

        PrecomputedRangeKind::Mlbw { range, l_groups } => {
            // MLBW uses the SAME precomputed J-groups as SLBW but a
            // different evaluator (coherent-sum elastic).  Routing MLBW
            // through `slbw_evaluate_with_cached_jgroups` was the #465 bug.
            let pi_over_k2 = channel::pi_over_k_squared_barns(energy_ev, awr);
            let mut total = 0.0;
            let mut elastic = 0.0;
            let mut capture = 0.0;
            let mut fission = 0.0;

            for lg in l_groups {
                let scatt_radius = if lg.apl > 0.0 {
                    lg.apl
                } else {
                    range.scattering_radius_at(energy_ev)
                };
                let pen_radius = if lg.pen_radius_override > 0.0 {
                    lg.pen_radius_override
                } else {
                    scatt_radius
                };

                let rho_phase = channel::rho(energy_ev, lg.awr_l, scatt_radius);
                let rho_pen = channel::rho(energy_ev, lg.awr_l, pen_radius);
                let phi = penetrability::phase_shift(lg.l, rho_phase);
                let p_at_e = penetrability::penetrability(lg.l, rho_pen);

                let (t, e, c, f) = slbw::mlbw_evaluate_with_cached_jgroups(
                    &lg.jgroups,
                    energy_ev,
                    pi_over_k2,
                    p_at_e,
                    phi,
                );
                total += t;
                elastic += e;
                capture += c;
                fission += f;
            }

            (total, elastic, capture, fission)
        }

        PrecomputedRangeKind::ReichMoore { range, l_groups } => {
            let mut total = 0.0;
            let mut elastic = 0.0;
            let mut capture = 0.0;
            let mut fission = 0.0;

            for lg in l_groups {
                let (l, awr_l, apl, pen_ovr) = match lg {
                    PrecomputedRmLGroupData::Single {
                        l,
                        awr_l,
                        apl,
                        pen_radius_override,
                        ..
                    } => (*l, *awr_l, *apl, *pen_radius_override),
                    PrecomputedRmLGroupData::TwoCh {
                        l,
                        awr_l,
                        apl,
                        pen_radius_override,
                        ..
                    } => (*l, *awr_l, *apl, *pen_radius_override),
                    PrecomputedRmLGroupData::ThreeCh {
                        l,
                        awr_l,
                        apl,
                        pen_radius_override,
                        ..
                    } => (*l, *awr_l, *apl, *pen_radius_override),
                };

                // Scattering radius for phase shift (always AP/APL).
                let scatt_radius = if apl > 0.0 {
                    apl
                } else {
                    range.scattering_radius_at(energy_ev)
                };
                // Penetrability/shift radius: precomputed override (APL or
                // NAPS=0 formula), falling back to scattering radius.
                let pen_radius = if pen_ovr > 0.0 { pen_ovr } else { scatt_radius };

                let rho_phase = channel::rho(energy_ev, awr_l, scatt_radius);
                let rho_pen = channel::rho(energy_ev, awr_l, pen_radius);
                let p_l = penetrability::penetrability(l, rho_pen);
                let s_l = penetrability::shift_factor(l, rho_pen);
                let phi_l = penetrability::phase_shift(l, rho_phase);

                let (t, e, c, f) = match lg {
                    PrecomputedRmLGroupData::Single { l, jgroups, .. } => {
                        let mut t = 0.0;
                        let mut e = 0.0;
                        let mut c = 0.0;
                        let mut f = 0.0;
                        for jg in jgroups {
                            // SAFETY: Float J comparison is safe here because both R-matrix resonance J values
                            // and R-external J values originate from the same `compute_j_offsets()` map in
                            // sammy.rs and follow the same computation path. The possible J offsets differ by
                            // multiples of ~1e-6, which is many orders of magnitude larger than the 1e-10
                            // QUANTUM_NUMBER_EPS used here, so the comparison reliably identifies matching J values.
                            let r_ext = range
                                .r_external
                                .iter()
                                .find(|re| re.l == *l && (re.j - jg.j).abs() < QUANTUM_NUMBER_EPS)
                                .map(|re| re.evaluate(energy_ev))
                                .unwrap_or(0.0);
                            let (jt, je, jc, jf) = reich_moore_spin_group_precomputed(
                                &jg.resonances,
                                energy_ev,
                                awr_l,
                                jg.g_j,
                                p_l,
                                s_l,
                                phi_l,
                                r_ext,
                            );
                            t += jt;
                            e += je;
                            c += jc;
                            f += jf;
                        }
                        (t, e, c, f)
                    }
                    PrecomputedRmLGroupData::TwoCh { jgroups, l, .. } => {
                        let mut t = 0.0;
                        let mut e = 0.0;
                        let mut c = 0.0;
                        let mut f = 0.0;
                        for jg in jgroups {
                            // P-7: R-external for 2ch fission path.
                            let r_ext = range
                                .r_external
                                .iter()
                                .find(|re| re.l == *l && (re.j - jg.j).abs() < QUANTUM_NUMBER_EPS)
                                .map(|re| re.evaluate(energy_ev))
                                .unwrap_or(0.0);
                            let (jt, je, jc, jf) = reich_moore_2ch_precomputed(
                                &jg.resonances,
                                energy_ev,
                                awr_l,
                                jg.g_j,
                                p_l,
                                s_l,
                                phi_l,
                                r_ext,
                            );
                            t += jt;
                            e += je;
                            c += jc;
                            f += jf;
                        }
                        (t, e, c, f)
                    }
                    PrecomputedRmLGroupData::ThreeCh { jgroups, l, .. } => {
                        let mut t = 0.0;
                        let mut e = 0.0;
                        let mut c = 0.0;
                        let mut f = 0.0;
                        for jg in jgroups {
                            // P-7: R-external for 3ch fission path.
                            let r_ext = range
                                .r_external
                                .iter()
                                .find(|re| re.l == *l && (re.j - jg.j).abs() < QUANTUM_NUMBER_EPS)
                                .map(|re| re.evaluate(energy_ev))
                                .unwrap_or(0.0);
                            let (jt, je, jc, jf) = reich_moore_3ch_precomputed(
                                &jg.resonances,
                                energy_ev,
                                awr_l,
                                jg.g_j,
                                p_l,
                                s_l,
                                phi_l,
                                r_ext,
                            );
                            t += jt;
                            e += je;
                            c += jc;
                            f += jf;
                        }
                        (t, e, c, f)
                    }
                };

                total += t;
                elastic += e;
                capture += c;
                fission += f;
            }

            (total, elastic, capture, fission)
        }
    }
}

/// Can this range actually produce non-zero cross-sections?
///
/// Delegates to [`ResonanceRange::is_evaluable`]: evaluable = resolved
/// LRF=1/2/3 (SLBW, MLBW, Reich-Moore). LRF=7 and LRU=2 ranges are
/// parse-and-skip placeholders and are never evaluated.
///
/// **In literal sync with `precompute_range_data`**: both key on
/// [`ResonanceRange::is_evaluable`] (`precompute_range_data` returns its Skip
/// arm for every non-evaluable range), so the energy-boundary logic
/// (`next_starts_here`) and the evaluator dispatch cannot disagree. Whenever
/// a new formalism becomes evaluable, extend `ResonanceRange::is_evaluable`.
fn range_is_evaluable(range: &ResonanceRange) -> bool {
    range.is_evaluable()
}

/// Cross-sections for a single spin group (J, π) in the Reich-Moore formalism,
/// using pre-computed per-resonance invariants (γ²_n cached).
///
/// For non-fissile isotopes, the R-matrix has a single neutron channel
/// and the capture channel is eliminated (absorbed into the imaginary
/// part of the resonance denominator).
///
/// ## Mathematical Formulation
///
/// For a single neutron channel with eliminated capture:
///
/// R(E) = Σ_n γ²_n / (E_n - E - iΓ_γ,n/2)
///
/// where γ²_n = Γ_n,n / (2·P_l(E_n)) is the reduced width amplitude squared.
///
/// Level matrix (scalar): Y = (S - B + iP)⁻¹ - R
///
/// X-matrix (scalar): X = P · Y⁻¹ · R · (S - B + iP)⁻¹
///
/// The scattering matrix element is:
///   U = e^{-2iφ} · (1 + 2i·X)
///
/// Cross-sections:
///   σ_elastic = (π/k²) · g_J · |1 - U|²
///   σ_total   = (2π/k²) · g_J · (1 - Re(U))
///   σ_capture = σ_total - σ_elastic (unitarity deficit)
///
/// Reference: SAMMY `rml/mrml11.f` Sectio routine
#[allow(clippy::too_many_arguments)]
fn reich_moore_spin_group_precomputed(
    resonances: &[PrecomputedResonanceSingle],
    energy_ev: f64,
    awr: f64,
    g_j: f64,
    p_l: f64,
    s_l: f64,
    phi_l: f64,
    r_ext: f64,
) -> (f64, f64, f64, f64) {
    let pi_over_k2 = channel::pi_over_k_squared_barns(energy_ev, awr);

    // Single-channel case (neutron only, capture eliminated).
    // This is the common case for non-fissile isotopes.

    // Boundary condition B = S_l(E): the shift factor and boundary cancel.
    //
    // SAMMY convention (CalcShift=false / Ishift=0): the shift factor
    // is NOT computed in the level matrix — Pgh (src/xxx/mxxx8.f90)
    // leaves S-B = 0 for all L values when Ishift=0. The resonance
    // energies in the .par file are observed peak positions.
    //
    // ENDF-102 convention (LRF=3 Reich-Moore): B_l = S_l, giving
    // S_l - B_l = 0 identically. Resonance energies are "formal"
    // eigenvalues, but with B=S the observed peaks coincide.
    //
    // Net effect: l_real = S - B = 0 for all L, so the level-matrix
    // denominator is 0 + iP, regardless of orbital angular momentum.
    let boundary = s_l;

    // Build the R-matrix (scalar, complex) = Σ_n γ²_n / (E_n - E - iΓ_γ,n/2)
    //
    // Note: ENDF stores "observed" widths Γ_n. The reduced width amplitude is:
    //   γ²_n = Γ_n / (2 · P_l(ρ_n))
    // where ρ_n = k(E_n)·a, evaluated at the resonance energy.
    //
    // Issue #87: γ²_n is now pre-computed in PrecomputedResonanceSingle.
    //
    // Reference: SAMMY `rml/mrml03.f` Betset (lines 240-276)
    let mut r_real = 0.0;
    let mut r_imag = 0.0;

    for res in resonances {
        let e_r = res.energy;
        let gamma_g = res.gamma_g;
        let gamma_n_reduced_sq = res.gamma_n_reduced_sq;

        // Denominator: (E_n - E)² + (Γ_γ/2)²
        let de = e_r - energy_ev;
        let half_gg = gamma_g / 2.0;
        let denom = de * de + half_gg * half_gg;

        if denom > DIVISION_FLOOR {
            // R-matrix contribution:
            // R += γ²_n / (E_n - E - i·Γ_γ/2)
            //    = γ²_n · (E_n - E + i·Γ_γ/2) / denom
            r_real += gamma_n_reduced_sq * de / denom;
            r_imag += gamma_n_reduced_sq * half_gg / denom;
        }
    }

    // R-external: diagonal, real-valued background R-matrix correction.
    // Adds smooth energy-dependent contribution from distant resonances.
    // SAMMY Ref: mcro2.f90 Setr_Cro lines 180-193
    r_real += r_ext;

    // Level matrix Y = 1/(S - B + iP) - R  (scalar, complex)
    let l_real = s_l - boundary;
    let l_imag = p_l;
    let l_denom = l_real * l_real + l_imag * l_imag;
    if l_denom < LOG_FLOOR {
        return (0.0, 0.0, 0.0, 0.0);
    }

    // 1/(S - B + iP) = (S - B - iP) / |S - B + iP|²
    let l_inv_real = l_real / l_denom;
    let l_inv_imag = -l_imag / l_denom;

    let y_real = l_inv_real - r_real;
    let y_imag = l_inv_imag - r_imag;

    // Y⁻¹ = 1/Y
    let y_denom = y_real * y_real + y_imag * y_imag;
    if y_denom < LOG_FLOOR {
        return (0.0, 0.0, 0.0, 0.0);
    }
    let y_inv_real = y_real / y_denom;
    let y_inv_imag = -y_imag / y_denom;

    // X-matrix (scalar): X = P · Y⁻¹ · R · (1/(S-B+iP))
    // Actually: X = √P · Y⁻¹ · R · √P · (1/(S-B+iP))
    //
    // From SAMMY mrml11.f: XXXX = √P_J · (Y⁻¹·R)_JI · (√P_I / L_II)
    // For single channel: X = √P · Y⁻¹ · R · √P / L
    //                       = P · Y⁻¹ · R / (S-B+iP)
    //
    // Let's compute step by step:
    // 1. q = Y⁻¹ · R (complex multiply)
    let q_real = y_inv_real * r_real - y_inv_imag * r_imag;
    let q_imag = y_inv_real * r_imag + y_inv_imag * r_real;

    // 2. X = P · q / (S-B+iP) = P · q · (S-B-iP) / |S-B+iP|²
    let x_unscaled_real = q_real * l_real + q_imag * l_imag;
    let x_unscaled_imag = q_imag * l_real - q_real * l_imag;
    let x_real = p_l * x_unscaled_real / l_denom;
    let x_imag = p_l * x_unscaled_imag / l_denom;

    // Compute the collision matrix element U from X.
    //
    //   U = e^{-2iφ} · (1 + 2iX)
    //
    // The phase factor uses e^{-2iφ}, NOT e^{+2iφ}.  SAMMY's Cossin
    // subroutine (src/xxx/mxxx6.f90 line 8) generates cos(2φ) and sin(2φ),
    // and the Total subroutine (src/cro/mcro4.f90 line 232) combines them
    // as:  Re(U) = cos(2φ)·Wr + sin(2φ)·Wi = Re(e^{-2iφ} · W)
    //
    // This is the Ω² factor in Lane & Thomas: Ω = e^{-iφ}.
    //
    // Reference: ENDF-102 Section 2, Lane & Thomas R-matrix theory
    let x = Complex64::new(x_real, x_imag);
    let phase = Complex64::new((2.0 * phi_l).cos(), -(2.0 * phi_l).sin());
    let u = phase * (1.0 + 2.0 * Complex64::i() * x);

    // Cross-sections from the collision matrix U:
    //
    //   σ_total   = g_J · (2π/k²) · (1 - Re(U))
    //   σ_elastic = g_J · (π/k²) · |1 - U|²
    //   σ_capture = σ_total - σ_elastic  (unitarity deficit)
    //
    // Reference: standard R-matrix cross-section formulas
    let sigma_total = g_j * 2.0 * pi_over_k2 * (1.0 - u.re);
    let one_minus_u = 1.0 - u;
    let sigma_elastic = g_j * pi_over_k2 * one_minus_u.norm_sqr();
    let sigma_capture = sigma_total - sigma_elastic;

    // For non-fissile isotopes, all absorption is capture.
    (sigma_total, sigma_elastic, sigma_capture, 0.0)
}

/// Reich-Moore 2-channel (neutron + 1 fission) with pre-computed betas.
///
/// Reference: SAMMY `rml/mrml09.f` Twoch routine
#[allow(clippy::too_many_arguments)]
fn reich_moore_2ch_precomputed(
    resonances: &[PrecomputedResonance2ch],
    energy_ev: f64,
    awr: f64,
    g_j: f64,
    p_l: f64,
    s_l: f64,
    phi_l: f64,
    r_ext: f64,
) -> (f64, f64, f64, f64) {
    let pi_over_k2 = channel::pi_over_k_squared_barns(energy_ev, awr);
    // B = S_l(E) — see comment in reich_moore_spin_group_precomputed.
    let boundary = s_l;

    // 2-channel: neutron + one fission channel.
    // R-matrix is 2x2 complex.
    let mut r_mat = [[Complex64::new(0.0, 0.0); 2]; 2];

    for res in resonances {
        // Denominator: (E_n - E) - i*Gamma_g/2
        let de = res.energy - energy_ev;
        let half_gg = res.gamma_g / 2.0;
        let inv_denom = 1.0 / Complex64::new(de, -half_gg);

        // R_ij += beta_i * beta_j / denom
        let betas = [res.beta_n, res.beta_f];
        for i in 0..2 {
            for j in 0..2 {
                r_mat[i][j] += betas[i] * betas[j] * inv_denom;
            }
        }
    }

    // P-7: R-external adds to the neutron channel diagonal only.
    // This is a smooth energy-dependent background from distant resonances.
    // SAMMY ref: mcro2.f90 Setr_Cro lines 180-193.
    if r_ext != 0.0 {
        r_mat[0][0] += Complex64::new(r_ext, 0.0);
    }

    // Level matrix Y = diag(1/(S-B+iP)) - R
    // Channel 0 (neutron): L = S_l - B + i*P_l
    // Channel 1 (fission): L = 0 + i*1 (no penetrability, Pent=0)
    //   -> fission channel: P_f = 1, S_f = 0
    let l_n = Complex64::new(s_l - boundary, p_l);
    let l_f = Complex64::new(0.0, 1.0); // Fission: no barrier

    let l_inv = [1.0 / l_n, 1.0 / l_f];

    let mut y_mat = [[Complex64::new(0.0, 0.0); 2]; 2];
    for i in 0..2 {
        for j in 0..2 {
            y_mat[i][j] = -r_mat[i][j];
        }
        y_mat[i][i] += l_inv[i];
    }

    // Invert 2x2 Y-matrix.
    // Guard against singular matrix.
    let det = y_mat[0][0] * y_mat[1][1] - y_mat[0][1] * y_mat[1][0];
    if det.norm() < LOG_FLOOR {
        return (0.0, 0.0, 0.0, 0.0);
    }
    let inv_det = 1.0 / det;
    let y_inv = [
        [y_mat[1][1] * inv_det, -y_mat[0][1] * inv_det],
        [-y_mat[1][0] * inv_det, y_mat[0][0] * inv_det],
    ];

    // X-matrix (ENDF-102 Eq. 2.76):
    // W_ij = sqrt(P_i) * (L^{-1}_i) * (Y^{-1} * R)_ij * sqrt(P_j)
    //
    // The L^{-1} factor is applied to channel i (row), NOT channel j (column).
    // This matters for off-diagonal elements where L_n = iP ≠ L_f = i.
    // Ref: (I - RL)^{-1} = L^{-1} · Y^{-1}, so L^{-1} multiplies from the left.
    let mut q = [[Complex64::new(0.0, 0.0); 2]; 2];
    for i in 0..2 {
        for j in 0..2 {
            for k in 0..2 {
                q[i][j] += y_inv[i][k] * r_mat[k][j];
            }
        }
    }

    let sqrt_p = [p_l.sqrt(), 1.0]; // sqrt(P_n), sqrt(P_f)
    let l_vals = [l_n, l_f];
    let mut x_mat = [[Complex64::new(0.0, 0.0); 2]; 2];
    for i in 0..2 {
        for j in 0..2 {
            x_mat[i][j] = sqrt_p[i] * q[i][j] * sqrt_p[j] / l_vals[i];
        }
    }

    // Collision matrix U from X-matrix.
    // Phase: e^{-2iφ} for diagonal, e^{-iφ} for off-diagonal (one neutron leg).
    // See comment in reich_moore_spin_group_precomputed for SAMMY reference.
    let phase2 = Complex64::new((2.0 * phi_l).cos(), -(2.0 * phi_l).sin());
    let phase1 = Complex64::new(phi_l.cos(), -phi_l.sin());

    let u_nn = phase2 * (1.0 + 2.0 * Complex64::i() * x_mat[0][0]);
    let u_nf = phase1 * 2.0 * Complex64::i() * x_mat[0][1];

    // Cross-sections from U-matrix.
    let sigma_total = g_j * 2.0 * pi_over_k2 * (1.0 - u_nn.re);
    let sigma_elastic = g_j * pi_over_k2 * (1.0 - u_nn).norm_sqr();
    let sigma_fission = g_j * pi_over_k2 * u_nf.norm_sqr();
    let sigma_capture = sigma_total - sigma_elastic - sigma_fission;

    (sigma_total, sigma_elastic, sigma_capture, sigma_fission)
}

/// 3-channel Reich-Moore (neutron + 2 fission channels) with pre-computed betas.
#[allow(clippy::too_many_arguments)]
fn reich_moore_3ch_precomputed(
    resonances: &[PrecomputedResonance3ch],
    energy_ev: f64,
    awr: f64,
    g_j: f64,
    p_l: f64,
    s_l: f64,
    phi_l: f64,
    r_ext: f64,
) -> (f64, f64, f64, f64) {
    let pi_over_k2 = channel::pi_over_k_squared_barns(energy_ev, awr);
    // B = S_l(E) — see comment in reich_moore_spin_group_precomputed.
    let boundary = s_l;

    let mut r_mat = [[Complex64::new(0.0, 0.0); 3]; 3];

    for res in resonances {
        let de = res.energy - energy_ev;
        let half_gg = res.gamma_g / 2.0;
        let inv_denom = 1.0 / Complex64::new(de, -half_gg);

        let betas = [res.beta_n, res.beta_fa, res.beta_fb];
        for i in 0..3 {
            for j in 0..3 {
                r_mat[i][j] += betas[i] * betas[j] * inv_denom;
            }
        }
    }

    // P-7: R-external adds to the neutron channel diagonal only.
    if r_ext != 0.0 {
        r_mat[0][0] += Complex64::new(r_ext, 0.0);
    }

    // Level matrix Y.
    let l_n = Complex64::new(s_l - boundary, p_l);
    let l_f = Complex64::new(0.0, 1.0);
    let l_vals = [l_n, l_f, l_f];
    let l_inv: Vec<Complex64> = l_vals.iter().map(|&li| 1.0 / li).collect();

    let mut y_mat = [[Complex64::new(0.0, 0.0); 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            y_mat[i][j] = -r_mat[i][j];
        }
        y_mat[i][i] += l_inv[i];
    }

    // Invert 3x3 via cofactor expansion.
    let y_inv = match invert_3x3(y_mat) {
        Some(inv) => inv,
        None => return (0.0, 0.0, 0.0, 0.0),
    };

    // X-matrix.
    let sqrt_p = [p_l.sqrt(), 1.0, 1.0];
    let mut x_mat = [[Complex64::new(0.0, 0.0); 3]; 3];
    let mut q = [[Complex64::new(0.0, 0.0); 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            for k in 0..3 {
                q[i][j] += y_inv[i][k] * r_mat[k][j];
            }
        }
    }
    for i in 0..3 {
        for j in 0..3 {
            x_mat[i][j] = sqrt_p[i] * q[i][j] * sqrt_p[j] / l_vals[i];
        }
    }

    // Collision matrix U from X-matrix.
    // Phase: e^{-2iφ} for diagonal, e^{-iφ} for off-diagonal.
    let phase2 = Complex64::new((2.0 * phi_l).cos(), -(2.0 * phi_l).sin());
    let phase1 = Complex64::new(phi_l.cos(), -phi_l.sin());

    let u_nn = phase2 * (1.0 + 2.0 * Complex64::i() * x_mat[0][0]);
    let u_nf1 = phase1 * 2.0 * Complex64::i() * x_mat[0][1];
    let u_nf2 = phase1 * 2.0 * Complex64::i() * x_mat[0][2];

    // Cross-sections from U-matrix.
    let sigma_total = g_j * 2.0 * pi_over_k2 * (1.0 - u_nn.re);
    let sigma_elastic = g_j * pi_over_k2 * (1.0 - u_nn).norm_sqr();
    let sigma_fission = g_j * pi_over_k2 * (u_nf1.norm_sqr() + u_nf2.norm_sqr());
    let sigma_capture = sigma_total - sigma_elastic - sigma_fission;

    (sigma_total, sigma_elastic, sigma_capture, sigma_fission)
}

/// Invert a 3×3 complex matrix via cofactor expansion.
///
/// Returns `None` if the matrix is singular (|det| < LOG_FLOOR), preventing
/// NaN propagation from 1/det when det ≈ 0.
fn invert_3x3(m: [[Complex64; 3]; 3]) -> Option<[[Complex64; 3]; 3]> {
    let det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
        - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
        + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);

    if det.norm() < LOG_FLOOR {
        return None; // singular — caller returns zero cross-sections
    }

    let inv_det = 1.0 / det;

    let mut result = [[Complex64::new(0.0, 0.0); 3]; 3];
    result[0][0] = (m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv_det;
    result[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inv_det;
    result[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv_det;
    result[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inv_det;
    result[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv_det;
    result[1][2] = (m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inv_det;
    result[2][0] = (m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv_det;
    result[2][1] = (m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inv_det;
    result[2][2] = (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv_det;

    Some(result)
}

// J-group assembly is now done inline via the `precompute_jgroups_*` functions
// (Issue #87).  The old `group_by_j` import is no longer needed.

#[cfg(test)]
mod tests {
    use super::*;
    use nereids_endf::resonance::test_support::{
        SingleResonanceParams, single_resonance, u238_single_resonance, u238_with_formalism,
    };
    use nereids_endf::resonance::{LGroup, Resonance, ResonanceRange};

    #[test]
    fn test_capture_peak_single_resonance() {
        // U-238 6.674 eV resonance.
        // At the resonance energy, capture cross-section should peak at ~22,000 barns.
        let data = u238_single_resonance();

        let xs = cross_sections_at_energy(&data, 6.674);

        // The capture cross-section at peak should be approximately:
        // σ_c = g_J × π/k² × 4×Γ_n×Γ_γ / Γ² where Γ = Γ_n + Γ_γ
        // For the RM formalism the peak is very close to this BW estimate.
        // g_J = 1.0, π/k² ≈ 98,200 barns, Γ = 0.024493
        // σ_c ≈ 1.0 × 98200 × 4 × 1.493e-3 × 23.0e-3 / (24.493e-3)²
        //     ≈ 98200 × 0.2289 ≈ 22,478 barns
        assert!(
            xs.capture > 15000.0 && xs.capture < 30000.0,
            "Capture should be ~22000 barns, got {}",
            xs.capture
        );
        assert!(xs.total > xs.capture, "Total > capture");
        assert!(xs.elastic > 0.0, "Elastic should be positive");
        assert!(xs.fission.abs() < 1e-10, "No fission for U-238");
    }

    #[test]
    fn test_1_over_v_behavior() {
        // Far from resonances, capture cross-section should follow 1/v ∝ 1/√E.
        // The 6.674 eV resonance tail should dominate at low energies.
        let data = u238_single_resonance();

        let xs_01 = cross_sections_at_energy(&data, 0.1);
        let xs_04 = cross_sections_at_energy(&data, 0.4);

        // At low E, σ ∝ 1/√E, so σ(0.1)/σ(0.4) ≈ √(0.4/0.1) = 2.0
        let ratio = xs_01.capture / xs_04.capture;
        assert!(
            (ratio - 2.0).abs() < 0.3,
            "Expected ~2.0 for 1/v behavior, got {}",
            ratio
        );
    }

    #[test]
    fn test_cross_sections_positive() {
        // All cross-sections must be non-negative at all energies.
        let data = u238_single_resonance();

        for &e in &[0.01, 0.1, 1.0, 5.0, 6.0, 6.674, 7.0, 10.0, 100.0, 1000.0] {
            let xs = cross_sections_at_energy(&data, e);
            assert!(xs.total >= 0.0, "Total negative at E={}: {}", e, xs.total);
            assert!(
                xs.elastic >= 0.0,
                "Elastic negative at E={}: {}",
                e,
                xs.elastic
            );
            assert!(
                xs.capture >= -1e-10,
                "Capture negative at E={}: {}",
                e,
                xs.capture
            );
        }
    }

    /// Parse the full vendored U-238 ENDF and compute cross-sections.
    ///
    /// Validates against the SAMMY ex027 case (Doppler-broadened at 300 K),
    /// against which we compare unbroadened RM values that should bracket
    /// the broadened data. Fixture is shipped under this crate's
    /// `tests/data/u238_ex027.endf` (public-domain ENDF/B-VIII.0) so the
    /// gate runs even when the crate is built standalone (outside the
    /// workspace, where `examples/data/` is not packaged).  The original
    /// `examples/data/u238_ex027.endf` is kept for end-user example code.
    #[test]
    fn test_u238_full_endf_cross_sections() {
        let endf_path =
            std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/u238_ex027.endf");

        let endf_text = std::fs::read_to_string(&endf_path)
            .unwrap_or_else(|e| panic!("vendored U-238 fixture missing at {endf_path:?}: {e}"));
        let data = nereids_endf::parser::parse_endf_file2(&endf_text).unwrap();

        // Compute cross-sections at several energies near the 6.674 eV resonance.
        let energies = [1.0, 5.0, 6.0, 6.5, 6.674, 7.0, 8.0, 10.0, 20.0, 50.0, 100.0];

        for &e in &energies {
            let xs = cross_sections_at_energy(&data, e);
            // Basic sanity: all cross-sections non-negative.
            assert!(xs.total >= 0.0, "Total negative at E={}", e);
            assert!(xs.elastic >= 0.0, "Elastic negative at E={}", e);
            // Capture can be very slightly negative due to floating point.
            assert!(
                xs.capture >= -0.01,
                "Capture negative at E={}: {}",
                e,
                xs.capture
            );
        }

        // Check the 6.674 eV resonance peak.
        // With the full ENDF file (all resonances), the peak capture
        // should still be dominated by the 6.674 eV resonance.
        let xs_peak = cross_sections_at_energy(&data, 6.674);
        assert!(
            xs_peak.capture > 10000.0,
            "Capture at 6.674 eV should be >10,000 barns (got {})",
            xs_peak.capture
        );

        // The 20.87 eV resonance should also show a significant peak.
        let xs_20 = cross_sections_at_energy(&data, 20.87);
        assert!(
            xs_20.capture > 1000.0,
            "Capture at 20.87 eV should be >1,000 barns (got {})",
            xs_20.capture
        );

        // SAMMY ex027 broadened output at ~6.674 eV gives ~339 barns capture.
        // Our UNBROADENED result should be MUCH larger (since Doppler broadening
        // spreads the peak). This confirms we're computing the correct physics.
        assert!(
            xs_peak.capture > 339.0,
            "Unbroadened peak must exceed SAMMY broadened value"
        );
    }

    /// `cross_sections_at_energy` with an SLBW-formalism range must give
    /// the same result as `slbw::slbw_cross_sections`.
    #[test]
    fn test_dispatcher_slbw_matches_slbw_module() {
        let data = u238_with_formalism(ResonanceFormalism::SLBW);

        let test_energies = [0.1, 1.0, 5.0, 6.0, 6.674, 7.0, 10.0, 100.0];
        for &e in &test_energies {
            let via_dispatcher = cross_sections_at_energy(&data, e);
            let via_slbw = crate::slbw::slbw_cross_sections(&data, e);

            let eps = 1e-10;
            assert!(
                (via_dispatcher.total - via_slbw.total).abs() < eps,
                "total mismatch at {e} eV: dispatcher={} slbw={}",
                via_dispatcher.total,
                via_slbw.total
            );
            assert!(
                (via_dispatcher.capture - via_slbw.capture).abs() < eps,
                "capture mismatch at {e} eV: dispatcher={} slbw={}",
                via_dispatcher.capture,
                via_slbw.capture
            );
            assert!(
                (via_dispatcher.elastic - via_slbw.elastic).abs() < eps,
                "elastic mismatch at {e} eV: dispatcher={} slbw={}",
                via_dispatcher.elastic,
                via_slbw.elastic
            );
        }
    }

    /// For a **single isolated resonance**, MLBW and SLBW should produce
    /// identical results because the interference term (cross-resonance)
    /// vanishes when there is only one resonance per spin group.
    ///
    /// Capture and fission are always identical (interference only
    /// affects elastic).  Elastic should also match for single resonances.
    #[test]
    fn test_mlbw_single_resonance_matches_slbw() {
        let data_mlbw = u238_with_formalism(ResonanceFormalism::MLBW);
        let data_slbw = u238_with_formalism(ResonanceFormalism::SLBW);

        let test_energies = [1.0, 6.674, 10.0];
        for &e in &test_energies {
            let xs_mlbw = cross_sections_at_energy(&data_mlbw, e);
            let xs_slbw = cross_sections_at_energy(&data_slbw, e);

            // Capture and fission must be identical.
            assert!(
                (xs_mlbw.capture - xs_slbw.capture).abs() < 1e-10,
                "MLBW/SLBW capture mismatch at {e} eV: mlbw={} slbw={}",
                xs_mlbw.capture,
                xs_slbw.capture
            );

            // Elastic: for single resonance, MLBW formula should reduce
            // to SLBW because there are no cross-terms.
            let rel_diff = if xs_slbw.elastic.abs() > 1e-10 {
                (xs_mlbw.elastic - xs_slbw.elastic).abs() / xs_slbw.elastic
            } else {
                (xs_mlbw.elastic - xs_slbw.elastic).abs()
            };
            assert!(
                rel_diff < 0.01,
                "MLBW/SLBW elastic mismatch at {e} eV: mlbw={} slbw={} (rel_diff={rel_diff})",
                xs_mlbw.elastic,
                xs_slbw.elastic,
            );
        }

        // Sanity: peak capture at resonance energy should be large.
        let xs_peak = cross_sections_at_energy(&data_mlbw, 6.674);
        assert!(
            xs_peak.capture > 1000.0,
            "MLBW capture at 6.674 eV should be substantial (got {})",
            xs_peak.capture
        );
    }

    /// Singular Y-matrix guard: when the R-matrix contribution nearly
    /// cancels the L⁻¹ diagonal at a resonance energy, Y ≈ 0 and
    /// Y⁻¹ diverges.  The cross-sections must remain finite and
    /// non-negative (no NaN or Inf propagation).
    ///
    /// We construct a scenario where evaluation occurs exactly at E_r,
    /// maximizing the R-matrix contribution.  With an extremely narrow
    /// resonance (Γ_γ = 1e-15 eV), the imaginary denominator is tiny
    /// and the R-matrix peak is enormous, stressing the Y inversion.
    #[test]
    fn test_reich_moore_singular_y_matrix_guard() {
        // Extremely narrow resonance: Γ_γ = 1e-15 eV forces R-matrix
        // contribution to be enormous at E = E_r, pushing Y toward
        // singularity and exercising the y_denom < LOG_FLOOR guard.
        let data = single_resonance(SingleResonanceParams {
            energy: 10.0,     // E_r
            gamma_n: 1.0e-3,  // Γ_n (eV)
            gamma_g: 1.0e-15, // Γ_γ (extremely small → near-singular Y)
            j: 0.5,
            l: 0,
            awr: 236.006,
            target_spin: 0.0,
            scattering_radius: 9.4285,
        });

        // Evaluate exactly at E_r where R is maximized.
        let xs = cross_sections_at_energy(&data, 10.0);
        assert!(
            xs.total.is_finite() && xs.total >= 0.0,
            "Total must be finite and non-negative at resonance peak, got {}",
            xs.total
        );
        assert!(
            xs.elastic.is_finite() && xs.elastic >= 0.0,
            "Elastic must be finite and non-negative, got {}",
            xs.elastic
        );
        assert!(
            xs.capture.is_finite(),
            "Capture must be finite, got {}",
            xs.capture
        );
    }

    /// Zero capture width definitively triggers the y_denom guard:
    /// with Γ_γ = 0, the R-matrix denominator at E = E_r is zero,
    /// making R infinite.  The singularity guard must return zeros
    /// rather than NaN/Inf.
    #[test]
    fn test_reich_moore_zero_capture_width_guard() {
        let data = single_resonance(SingleResonanceParams {
            energy: 10.0,    // E_r
            gamma_n: 1.0e-3, // Γ_n (eV)
            gamma_g: 0.0,    // Γ_γ = 0 → guaranteed singularity
            j: 0.5,
            l: 0,
            awr: 236.006,
            target_spin: 0.0,
            scattering_radius: 9.4285,
        });

        // At E = E_r with Γ_γ = 0, the denominator (E_r - E)² + (Γ_γ/2)² = 0,
        // so the DIVISION_FLOOR guard on the R-matrix denom fires, but even if
        // it didn't, the y_denom guard would catch it downstream.
        let xs = cross_sections_at_energy(&data, 10.0);
        assert!(
            xs.total.is_finite() && xs.total >= 0.0,
            "Total must be finite and non-negative with zero capture width, got {}",
            xs.total
        );
        assert!(
            xs.elastic.is_finite() && xs.elastic >= 0.0,
            "Elastic must be finite and non-negative, got {}",
            xs.elastic
        );
        // With Γ_γ = 0 the true capture is exactly zero, but floating-point
        // arithmetic at this singularity can produce a tiny negative value
        // (machine-epsilon level).  Accept values > -1e-10 barns.
        assert!(
            xs.capture.is_finite() && xs.capture > -1e-10,
            "Capture must be finite and nearly non-negative, got {}",
            xs.capture
        );
    }

    /// Reich-Moore 2-channel (fission) with a singular det guard:
    /// when both fission and neutron widths are tiny, the 2x2 Y-matrix
    /// determinant can be near zero.  Results must be finite.
    #[test]
    fn test_reich_moore_fission_near_singular() {
        let data = ResonanceData {
            isotope: nereids_core::types::Isotope::new(94, 239).unwrap(),
            za: 94239,
            awr: 236.998,
            ranges: vec![ResonanceRange {
                energy_low: 1e-5,
                energy_high: 1e4,
                resolved: true,
                formalism: ResonanceFormalism::ReichMoore,
                target_spin: 0.5,
                scattering_radius: 9.41,
                naps: 1,
                l_groups: vec![LGroup {
                    l: 0,
                    awr: 236.998,
                    apl: 0.0,
                    qx: 0.0,
                    lrx: 0,
                    resonances: vec![Resonance {
                        energy: 10.0,
                        j: 1.0,
                        gn: 1.0e-8,  // very small neutron width
                        gg: 1.0e-8,  // very small capture width
                        gfa: 1.0e-8, // very small fission width
                        gfb: 0.0,
                    }],
                }],
                ap_table: None,
                r_external: vec![],
            }],
        };

        // Evaluate at the resonance energy.
        let xs = cross_sections_at_energy(&data, 10.0);
        assert!(
            xs.total.is_finite() && xs.total >= 0.0,
            "Total must be finite, got {}",
            xs.total
        );
        assert!(
            xs.fission.is_finite() && xs.fission >= 0.0,
            "Fission must be finite and non-negative, got {}",
            xs.fission
        );
    }

    /// A parsed-but-skipped range contributes EXACTLY zero cross-section over
    /// its span, through both evaluation entry points.
    ///
    /// This pins the branch's central safety property: a non-evaluable
    /// placeholder range (LRU=2 URR / LRF=7 RML / LRU=0) reached by real mixed
    /// tapes must add nothing to the four components, so the load-time
    /// warnings ("these spans contribute zero cross-section") are honoured at
    /// evaluation time — not just at parse time. The data has one evaluable
    /// resolved range [1e-5, 1e4] eV (a real U-238 6.674 eV resonance) and one
    /// disjoint placeholder range [1e4, 1e5] eV.
    ///
    /// Non-circular: the same shared primitive (`evaluate_precomputed_range`)
    /// backs both APIs, so a vacuous "everything is zero" implementation would
    /// pass a zeros-only check. The nonzero assertion inside the evaluable span
    /// rules that out; the boundary assertion pins that the skipped range does
    /// not corrupt the resolved range's inclusive upper edge.
    #[test]
    fn test_skipped_range_contributes_exact_zero() {
        // Evaluable resolved RM range [1e-5, 1e4] with one resonance.
        let mut data = u238_with_formalism(ResonanceFormalism::ReichMoore);
        // Append a disjoint non-evaluable placeholder over [1e4, 1e5] (URR-like:
        // empty l_groups, resolved=false), mirroring a real mixed tape.
        data.ranges.push(ResonanceRange {
            energy_low: 1e4,
            energy_high: 1e5,
            resolved: false,
            formalism: ResonanceFormalism::Unresolved,
            target_spin: 0.0,
            scattering_radius: 9.4285,
            naps: 1,
            ap_table: None,
            l_groups: vec![],
            r_external: vec![],
        });

        assert!(data.has_unevaluated_ranges());
        assert_eq!(data.unevaluated_ranges().len(), 1);

        // (1) Non-vacuity: the evaluable span must produce a real signal, so
        // "zero inside the skipped span" is not trivially true everywhere.
        let e_eval = 6.674; // strictly inside [1e-5, 1e4]
        let xs_eval = cross_sections_at_energy(&data, e_eval);
        assert!(
            xs_eval.total > 0.0 && xs_eval.capture > 0.0,
            "evaluable span must be nonzero, got {xs_eval:?}"
        );

        // (2) Exact zero strictly inside the placeholder span, via BOTH APIs.
        let e_skip = [2.0e4, 5.0e4, 9.0e4]; // strictly inside [1e4, 1e5]
        for &e in &e_skip {
            let pt = cross_sections_at_energy(&data, e);
            assert_eq!(pt.total, 0.0, "per-point total must be exactly 0 at E={e}");
            assert_eq!(
                pt.elastic, 0.0,
                "per-point elastic must be exactly 0 at E={e}"
            );
            assert_eq!(
                pt.capture, 0.0,
                "per-point capture must be exactly 0 at E={e}"
            );
            assert_eq!(
                pt.fission, 0.0,
                "per-point fission must be exactly 0 at E={e}"
            );
        }
        for (i, xs) in cross_sections_on_grid(&data, &e_skip).iter().enumerate() {
            let e = e_skip[i];
            assert_eq!(xs.total, 0.0, "grid total must be exactly 0 at E={e}");
            assert_eq!(xs.elastic, 0.0, "grid elastic must be exactly 0 at E={e}");
            assert_eq!(xs.capture, 0.0, "grid capture must be exactly 0 at E={e}");
            assert_eq!(xs.fission, 0.0, "grid fission must be exactly 0 at E={e}");
        }

        // (3) Scalar/grid agreement at shared points (evaluable and skipped).
        let shared = [6.674, 2.0e4, 5.0e4];
        let grid = cross_sections_on_grid(&data, &shared);
        for (i, &e) in shared.iter().enumerate() {
            let pt = cross_sections_at_energy(&data, e);
            assert_eq!(
                pt.total, grid[i].total,
                "scalar/grid total disagree at E={e}"
            );
            assert_eq!(
                pt.elastic, grid[i].elastic,
                "scalar/grid elastic disagree at E={e}"
            );
            assert_eq!(
                pt.capture, grid[i].capture,
                "scalar/grid capture disagree at E={e}"
            );
            assert_eq!(
                pt.fission, grid[i].fission,
                "scalar/grid fission disagree at E={e}"
            );
        }

        // (4) Boundary: the placeholder is non-evaluable, so the resolved range
        // owns its inclusive upper edge 1e4. At exactly 1e4 the placeholder adds
        // zero — the value must equal the resolved-range-only value there.
        let e_edge = 1e4;
        let edge = cross_sections_at_energy(&data, e_edge);
        let mut resolved_only = data.clone();
        resolved_only.ranges.truncate(1);
        let edge_resolved = cross_sections_at_energy(&resolved_only, e_edge);
        assert_eq!(
            edge.total, edge_resolved.total,
            "skipped range corrupts the boundary total"
        );
        assert_eq!(edge.elastic, edge_resolved.elastic);
        assert_eq!(edge.capture, edge_resolved.capture);
        assert_eq!(edge.fission, edge_resolved.fission);
    }

    /// `cross_sections_on_grid` (batch, precompute hoisted) must produce
    /// identical results to `cross_sections_at_energy` (per-point) for
    /// Reich-Moore data.
    #[test]
    fn test_grid_matches_per_point_reich_moore() {
        let data = u238_single_resonance();

        let energies = [0.01, 0.1, 1.0, 5.0, 6.0, 6.674, 7.0, 10.0, 100.0, 1000.0];
        let grid_results = cross_sections_on_grid(&data, &energies);

        for (i, &e) in energies.iter().enumerate() {
            let point = cross_sections_at_energy(&data, e);
            let grid = &grid_results[i];
            let eps = 1e-12;
            assert!(
                (point.total - grid.total).abs() < eps,
                "total mismatch at E={e}: per_point={} grid={}",
                point.total,
                grid.total
            );
            assert!(
                (point.elastic - grid.elastic).abs() < eps,
                "elastic mismatch at E={e}: per_point={} grid={}",
                point.elastic,
                grid.elastic
            );
            assert!(
                (point.capture - grid.capture).abs() < eps,
                "capture mismatch at E={e}: per_point={} grid={}",
                point.capture,
                grid.capture
            );
            assert!(
                (point.fission - grid.fission).abs() < eps,
                "fission mismatch at E={e}: per_point={} grid={}",
                point.fission,
                grid.fission
            );
        }
    }

    /// `cross_sections_on_grid` must match `cross_sections_at_energy` for
    /// SLBW-formalism data too (the batch path precomputes SLBW J-groups).
    #[test]
    fn test_grid_matches_per_point_slbw() {
        let data = u238_with_formalism(ResonanceFormalism::SLBW);

        let energies = [0.1, 1.0, 5.0, 6.0, 6.674, 7.0, 10.0, 100.0];
        let grid_results = cross_sections_on_grid(&data, &energies);

        for (i, &e) in energies.iter().enumerate() {
            let point = cross_sections_at_energy(&data, e);
            let grid = &grid_results[i];
            let eps = 1e-12;
            assert!(
                (point.total - grid.total).abs() < eps,
                "total mismatch at E={e}: per_point={} grid={}",
                point.total,
                grid.total
            );
            assert!(
                (point.capture - grid.capture).abs() < eps,
                "capture mismatch at E={e}: per_point={} grid={}",
                point.capture,
                grid.capture
            );
        }
    }

    /// `cross_sections_on_grid` must match per-point for the full U-238
    /// ENDF file (many resonances, L-groups, J-groups). Fixture is the
    /// crate-local `tests/data/u238_ex027.endf`, so the gate runs
    /// unconditionally on CI even when the crate is built standalone
    /// (outside the workspace, where `examples/data/` is not packaged).
    #[test]
    fn test_grid_matches_per_point_u238_full() {
        let endf_path =
            std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/u238_ex027.endf");

        let endf_text = std::fs::read_to_string(&endf_path)
            .unwrap_or_else(|e| panic!("vendored U-238 fixture missing at {endf_path:?}: {e}"));
        let data = nereids_endf::parser::parse_endf_file2(&endf_text).unwrap();

        let energies: Vec<f64> = (0..100).map(|i| 1.0 + i as f64 * 0.5).collect();
        let grid_results = cross_sections_on_grid(&data, &energies);

        for (i, &e) in energies.iter().enumerate() {
            let point = cross_sections_at_energy(&data, e);
            let grid = &grid_results[i];
            let eps = 1e-10;
            assert!(
                (point.total - grid.total).abs() < eps * point.total.abs().max(1.0),
                "total mismatch at E={e}: per_point={} grid={}",
                point.total,
                grid.total
            );
            assert!(
                (point.capture - grid.capture).abs() < eps * point.capture.abs().max(1.0),
                "capture mismatch at E={e}: per_point={} grid={}",
                point.capture,
                grid.capture
            );
        }
    }

    /// Verify that NAPS=0 uses the channel radius formula for penetrability
    /// while still using AP for phase shifts.
    ///
    /// The NAPS flag controls which radius is used for penetrability
    /// and shift factor calculations (ENDF-6 §2.2.1):
    /// - NAPS=0: channel radius = (0.123·A^(1/3) + 0.08) × 10  (fm)
    /// - NAPS=1: scattering radius AP (or AP(E))
    ///
    /// Note: for L=0, P_0(rho) = 1 regardless of radius, so NAPS only
    /// affects L>=1.  Even for L>=1, the neutron width uses the RATIO
    /// P_l(E)/P_l(E_r), so the radius effect largely cancels.  The main
    /// observable impact of NAPS is through the penetrability and
    /// shift-factor terms (P_l, S_l) for L>=1; the phase shift itself
    /// always uses the scattering radius AP regardless of NAPS.
    ///
    /// This test verifies the code path is wired correctly by checking:
    /// 1. NAPS=0 with AP = formula_radius matches NAPS=1 with AP = formula_radius
    ///    (confirming the formula gives the expected value)
    /// 2. NAPS=0 with a different AP still produces valid XS (no NaN/panic)
    #[test]
    fn test_naps_zero_uses_channel_radius_formula() {
        let awr: f64 = 55.345; // Fe-56-like
        let formula_radius = channel::endf_channel_radius_fm(awr);

        // NAPS=0: penetrability uses formula, phase shift uses AP (= formula here)
        let data_naps0 = ResonanceData {
            isotope: nereids_core::types::Isotope::new(26, 56).unwrap(),
            za: 26056,
            awr,
            ranges: vec![ResonanceRange {
                energy_low: 1e-5,
                energy_high: 1e5,
                resolved: true,
                formalism: ResonanceFormalism::ReichMoore,
                target_spin: 0.0,
                scattering_radius: formula_radius, // AP = formula → same as pen_radius
                naps: 0,
                l_groups: vec![LGroup {
                    l: 1,
                    awr,
                    apl: 0.0,
                    qx: 0.0,
                    lrx: 0,
                    resonances: vec![Resonance {
                        energy: 30000.0,
                        j: 1.5,
                        gn: 5.0,
                        gg: 1.0,
                        gfa: 0.0,
                        gfb: 0.0,
                    }],
                }],
                ap_table: None,
                r_external: vec![],
            }],
        };

        // NAPS=1 with AP = formula_radius: both penetrability and phase use AP
        let data_naps1 = ResonanceData {
            isotope: nereids_core::types::Isotope::new(26, 56).unwrap(),
            za: 26056,
            awr,
            ranges: vec![ResonanceRange {
                energy_low: 1e-5,
                energy_high: 1e5,
                resolved: true,
                formalism: ResonanceFormalism::ReichMoore,
                target_spin: 0.0,
                scattering_radius: formula_radius,
                naps: 1,
                l_groups: vec![LGroup {
                    l: 1,
                    awr,
                    apl: 0.0,
                    qx: 0.0,
                    lrx: 0,
                    resonances: vec![Resonance {
                        energy: 30000.0,
                        j: 1.5,
                        gn: 5.0,
                        gg: 1.0,
                        gfa: 0.0,
                        gfb: 0.0,
                    }],
                }],
                ap_table: None,
                r_external: vec![],
            }],
        };

        let e = 30000.0;
        let xs_naps0 = cross_sections_at_energy(&data_naps0, e);
        let xs_naps1 = cross_sections_at_energy(&data_naps1, e);

        // When AP equals the formula radius, NAPS=0 and NAPS=1 should
        // give identical results (both use the same radius everywhere).
        assert!(
            (xs_naps0.total - xs_naps1.total).abs() < 1e-10 * xs_naps1.total.abs().max(1.0),
            "NAPS=0 total={} vs NAPS=1 total={}: should match when AP=formula",
            xs_naps0.total,
            xs_naps1.total,
        );
        assert!(
            (xs_naps0.capture - xs_naps1.capture).abs() < 1e-10 * xs_naps1.capture.abs().max(1.0),
            "NAPS=0 capture={} vs NAPS=1 capture={}: should match when AP=formula",
            xs_naps0.capture,
            xs_naps1.capture,
        );

        // Verify finite and positive cross-sections (no NaN from formula)
        assert!(xs_naps0.total.is_finite() && xs_naps0.total > 0.0);
        assert!(xs_naps0.capture.is_finite() && xs_naps0.capture > 0.0);

        // Also verify the formula value is sane
        assert!(
            (formula_radius - 5.49).abs() < 0.1,
            "Expected ~5.49 fm for Fe-56-like, got {formula_radius}"
        );
    }

    // ─── Issue #465: batch vs per-point equivalence across formalisms ────────
    //
    // These tests lock the contract that `cross_sections_on_grid` must
    // produce element-wise bit-exact output matching
    // `cross_sections_at_energy` on the same grid.  Two paths diverging
    // silently is the class of bug that #465 exposed — the batch path was
    // routing MLBW ranges through the SLBW incoherent-sum evaluator,
    // producing up to 55 % relative error on real Hf isotopes — so this
    // harness must cover every formalism that has a distinct evaluator.
    //
    // The synthetic MLBW test below uses the smallest configuration that
    // exercises the coherent-vs-incoherent divergence (a single J-group
    // with ≥2 resonances).  The Hf-177 test provides a real-world anchor
    // against ENDF-B-VIII.1 data committed under `tests/data/endf/`.

    /// Synthetic MLBW fixture: the batch API and per-point API MUST agree
    /// bit-exactly.  This is the root cause of #465 — the batch
    /// dispatcher lumped MLBW into the SLBW evaluator (incoherent sum)
    /// while the per-point dispatcher correctly routed through the
    /// coherent-sum MLBW evaluator.  Both paths now share
    /// `slbw::mlbw_evaluate_with_cached_jgroups`.
    ///
    /// Uses the Hf-177 high-J fixture from `test_support`: two s-wave
    /// resonances at 2.386 eV and 5.89 eV in the same J = 4.0 group
    /// (`hf177_mlbw_two_resonances_high_j`).
    #[test]
    fn test_batch_matches_per_point_mlbw_synthetic() {
        let data = nereids_endf::resonance::test_support::hf177_mlbw_two_resonances_high_j();
        // Sample densely across the 2.386 eV and 5.89 eV resonances and
        // their interference region; also cover tail behaviour.
        let energies: Vec<f64> = (0..101).map(|i| 0.5 + (i as f64) * 0.1).collect();

        let per_point: Vec<CrossSections> = energies
            .iter()
            .map(|&e| cross_sections_at_energy(&data, e))
            .collect();
        let batch = cross_sections_on_grid(&data, &energies);

        assert_eq!(per_point.len(), batch.len());
        for (i, (pp, b)) in per_point.iter().zip(batch.iter()).enumerate() {
            // Bit-exact equality — same math, same inputs, same
            // accumulation order, so f64::to_bits() must match.
            assert_eq!(
                pp.total.to_bits(),
                b.total.to_bits(),
                "total mismatch at E[{i}]={}: per_point={:.17e}  batch={:.17e}  rel_diff={:.3e}",
                energies[i],
                pp.total,
                b.total,
                if pp.total != 0.0 {
                    (pp.total - b.total).abs() / pp.total.abs()
                } else {
                    (pp.total - b.total).abs()
                },
            );
            assert_eq!(
                pp.elastic.to_bits(),
                b.elastic.to_bits(),
                "elastic mismatch at E[{i}]={}: per_point={:.17e}  batch={:.17e}",
                energies[i],
                pp.elastic,
                b.elastic,
            );
            assert_eq!(pp.capture.to_bits(), b.capture.to_bits());
            assert_eq!(pp.fission.to_bits(), b.fission.to_bits());
        }
    }

    /// Real-world anchor: ENDF-B-VIII.1 Hf-177 is MLBW (LRF=2) with 180
    /// resonances.  Locks the #465 production symptom (up to 55 %
    /// relative divergence on the VENUS analysis grid).
    ///
    /// The ENDF file is committed under `tests/data/endf/Hf-177.endf`
    /// at the workspace root (~2.8 MB, public domain, see the README
    /// there).  The fixture lives outside `crates/nereids-physics/` so
    /// it is not included by `cargo package` when this crate is
    /// published in isolation.  When the fixture is absent the test
    /// skips with a note rather than panicking, so standalone crate
    /// checkouts still see a clean `cargo test` run.  In the full
    /// workspace (where the fixture is always present) the gate runs
    /// as normal.
    #[test]
    fn test_batch_matches_per_point_hf177_real_endf() {
        use nereids_endf::parser::parse_endf_file2;
        let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .parent()
            .unwrap()
            .parent()
            .unwrap()
            .join("tests/data/endf/Hf-177.endf");
        let text = match std::fs::read_to_string(&path) {
            Ok(t) => t,
            Err(e) => {
                eprintln!(
                    "skipping test_batch_matches_per_point_hf177_real_endf: \
                     fixture not available at {path:?}: {e}. \
                     Run from the full NEREIDS workspace to exercise this regression gate."
                );
                return;
            }
        };
        let data = parse_endf_file2(&text).unwrap();

        // 500 points spanning the resolved MLBW range (up to 250 eV).
        // Focuses coverage on the analysis-relevant VENUS grid segment.
        let n = 500;
        let energies: Vec<f64> = (0..n)
            .map(|i| 0.5 + (i as f64) * ((200.0 - 0.5) / (n - 1) as f64))
            .collect();

        let per_point: Vec<f64> = energies
            .iter()
            .map(|&e| cross_sections_at_energy(&data, e).total)
            .collect();
        let batch: Vec<f64> = cross_sections_on_grid(&data, &energies)
            .into_iter()
            .map(|cs| cs.total)
            .collect();

        // On MLBW data the legacy batch path differs by up to 55 %.  The
        // fix must bring this to bit-exact.  To keep the assertion
        // focused (and avoid cascade failures that hide other issues),
        // count mismatches and report the worst offender before asserting.
        let mut max_rel = 0.0f64;
        let mut worst_idx = 0usize;
        let mut mismatches = 0usize;
        for (i, (&a, &b)) in per_point.iter().zip(batch.iter()).enumerate() {
            if a.to_bits() != b.to_bits() {
                mismatches += 1;
            }
            let rel = if a != 0.0 {
                (a - b).abs() / a.abs()
            } else {
                (a - b).abs()
            };
            if rel > max_rel {
                max_rel = rel;
                worst_idx = i;
            }
        }
        assert_eq!(
            mismatches, 0,
            "{mismatches}/{n} points differ; worst at E[{worst_idx}]={} — \
             per_point={:.6e}  batch={:.6e}  rel_diff={max_rel:.3e}",
            energies[worst_idx], per_point[worst_idx], batch[worst_idx],
        );
    }

    // ─── Top-level pub-fn energy-validation guards ─────────────────────────
    //
    // Mirrors the SLBW / Reich-Moore test patterns: NaN, ±Inf, 0, and -1 must
    // each panic with the canonical "expected positive finite energy_ev"
    // message.  These tests gate the symmetric defense-in-depth contract on
    // both Reich-Moore top-level pub fns so a future refactor cannot
    // silently drop the assert.

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_at_energy_panics_on_nan() {
        let data = u238_single_resonance();
        let _ = cross_sections_at_energy(&data, f64::NAN);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_at_energy_panics_on_infinity() {
        let data = u238_single_resonance();
        let _ = cross_sections_at_energy(&data, f64::INFINITY);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_at_energy_panics_on_zero() {
        let data = u238_single_resonance();
        let _ = cross_sections_at_energy(&data, 0.0);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_at_energy_panics_on_negative() {
        let data = u238_single_resonance();
        let _ = cross_sections_at_energy(&data, -1.0);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_on_grid_panics_on_nan() {
        let data = u238_single_resonance();
        let _ = cross_sections_on_grid(&data, &[1.0, f64::NAN, 2.0]);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_on_grid_panics_on_infinity() {
        let data = u238_single_resonance();
        let _ = cross_sections_on_grid(&data, &[1.0, f64::INFINITY]);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_on_grid_panics_on_zero() {
        let data = u238_single_resonance();
        let _ = cross_sections_on_grid(&data, &[1.0, 0.0, 2.0]);
    }

    #[test]
    #[should_panic(expected = "expected positive finite energy_ev")]
    fn cross_sections_on_grid_panics_on_negative() {
        let data = u238_single_resonance();
        let _ = cross_sections_on_grid(&data, &[-1.0, 1.0, 2.0]);
    }

    /// The plan, the grid function and the per-point function are one
    /// evaluator: every channel agrees bit for bit, for every formalism.
    #[test]
    fn plan_grid_and_per_point_evaluators_are_bit_identical() {
        let energies: Vec<f64> = (1..=400).map(|i| 5.0 + 0.01 * f64::from(i)).collect();
        for formalism in [
            ResonanceFormalism::SLBW,
            ResonanceFormalism::MLBW,
            ResonanceFormalism::ReichMoore,
        ] {
            let data = u238_with_formalism(formalism);
            assert_bit_identical(&data, &energies, formalism);
        }
    }

    /// The same three evaluators across a bound SHARED by two evaluable
    /// ranges, which is the only arrangement that makes
    /// `upper_bound_is_half_open` answer `true`. The single-range fixture
    /// above can never reach that branch, so without this the consolidated
    /// `covers` predicate would be pinned only on its inclusive side and the
    /// three dispatch sites could drift apart exactly where the rule bites.
    #[test]
    fn the_three_evaluators_agree_on_a_shared_range_bound() {
        let mut data = u238_with_formalism(ResonanceFormalism::SLBW);
        data.ranges[0].energy_high = 100.0;
        let mut upper = u238_with_formalism(ResonanceFormalism::MLBW).ranges[0].clone();
        upper.energy_low = 100.0;
        upper.l_groups[0].resonances[0].energy = 150.0;
        data.ranges.push(upper);
        assert!(
            upper_bound_is_half_open(&data, 0),
            "the fixture must share a bound"
        );

        // Non-vacuity: the two ranges must disagree at the shared bound, or
        // agreeing on which one owns it would prove nothing. Each is asked
        // alone, with its own bound widened so that it does cover 100 eV.
        let alone = |index: usize| {
            let mut one = data.clone();
            one.ranges = vec![one.ranges[index].clone()];
            one.ranges[0].energy_low = 1e-5;
            one.ranges[0].energy_high = 1e4;
            cross_sections_at_energy(&one, 100.0)
        };
        let (lower_only, upper_only) = (alone(0), alone(1));
        assert_ne!(lower_only.total.to_bits(), upper_only.total.to_bits());

        let energies = [99.0, 99.999_999, 100.0, 100.000_001, 101.0];
        assert_bit_identical(&data, &energies, ResonanceFormalism::SLBW);

        // The bound itself belongs to the upper range alone.
        assert_eq!(
            cross_sections_at_energy(&data, 100.0).total.to_bits(),
            upper_only.total.to_bits()
        );
    }

    /// Every channel of the plan, the grid function and the per-point
    /// function, compared to the bit.
    fn assert_bit_identical(data: &ResonanceData, energies: &[f64], label: ResonanceFormalism) {
        let plan = CrossSectionPlan::new(data);
        let on_grid = cross_sections_on_grid(data, energies);
        for (&energy, grid) in energies.iter().zip(&on_grid) {
            let point = cross_sections_at_energy(data, energy);
            let one = plan.evaluate_one(energy);
            for (a, b, c) in [
                (grid.total, point.total, one.total),
                (grid.elastic, point.elastic, one.elastic),
                (grid.capture, point.capture, one.capture),
                (grid.fission, point.fission, one.fission),
            ] {
                assert_eq!(a.to_bits(), b.to_bits(), "{label:?} at {energy} eV");
                assert_eq!(a.to_bits(), c.to_bits(), "{label:?} at {energy} eV");
            }
        }
    }
}