nereids_endf/

parser.rs

//! ENDF-6 File 2 resonance parameter, File 3 cross section and File 32
//! covariance parser.
//!
//! Parses the fixed-width 80-character ENDF-6 format to extract resolved
//! resonance region (RRR) parameters.
//!
//! ## ENDF-6 Line Format
//! Each line is exactly 80 characters:
//! - Cols 1-11:  Field 1 (floating point or integer)
//! - Cols 12-22: Field 2
//! - Cols 23-33: Field 3
//! - Cols 34-44: Field 4
//! - Cols 45-55: Field 5
//! - Cols 56-66: Field 6
//! - Cols 67-70: MAT number
//! - Cols 71-72: MF (file number)
//! - Cols 73-75: MT (section number)
//! - Cols 76-80: Line sequence
//!
//! ## SAMMY Reference
//! - SAMMY manual Section 9 (ENDF-6 format)
//! - SAMMY source: `sammy/src/endf/` module

use crate::resonance::*;
use nereids_core::elements::isotope_from_za;

/// ENDF radius unit conversion factor.
///
/// ENDF-6 Formats Manual §2.1: "AP, APE, APT are in units of 10⁻¹² cm."
/// Physics convention: 1 fm = 10⁻¹³ cm, so 10⁻¹² cm = 10 fm.
/// All ENDF radii (AP, APL, APE, APT, AP(E) tables) are multiplied by this
/// factor at parse time so that downstream physics uses true femtometers.
///
/// SAMMY applies the identical ×10 conversion when reading ENDF:
///   `FillSammyRmatrixFromRMat.cpp` line 422: `newChannel->getApe() * 10.0`
const ENDF_RADIUS_TO_FM: f64 = 10.0;

/// Parse ENDF-6 File 2 resonance parameters, and File 3's cross sections
/// when the text holds them, from raw ENDF text.
///
/// Extracts all MF=2, MT=151 lines and parses the resolved resonance region;
/// File 3's sections MT=1, 2, 102 and 18 become [`ResonanceData::smooth`].
/// An evaluation File 1 marks LRP=2 is refused.
///
/// # Arguments
/// * `endf_text` — Full ENDF file contents as a string.
///
/// # Returns
/// `ResonanceData` containing all parsed resonance parameters.
pub fn parse_endf_file2(endf_text: &str) -> Result<ResonanceData, EndfParseError> {
    let data = parse_file2_records(endf_text)?;

    // Fail loudly when the evaluation carries NO evaluable content. A file
    // whose every range is a parse-and-skip placeholder (LRU=2, or an
    // LRU=0 scattering-radius-only stanza) would otherwise "load" with zero
    // resonances and add no resonance cross section anywhere, with no signal
    // to the caller. Files with at least one evaluable range still load:
    // parse-and-skip exists for real mixed tapes
    // (e.g. Ta-181, U-238, Pu-240) whose resolved ranges remain fully usable.
    if !data.has_evaluable_range() {
        let detail = if data.ranges.is_empty() {
            "the file carries no resonance-parameter ranges at all".to_string()
        } else {
            let skipped: Vec<String> = data.ranges.iter().map(|r| r.skip_description()).collect();
            format!(
                "every range in this file is a parse-and-skip placeholder: {}",
                skipped.join("; ")
            )
        };
        return Err(EndfParseError::UnsupportedFormat(format!(
            "No evaluable resonance ranges: NEREIDS evaluates resolved LRF=3 \
             (Reich-Moore) ranges, LRF=1/2 (SLBW/MLBW) ranges without a competitive \
             width (LRX=0), and LRF=7 (R-Matrix Limited) ranges in the Reich-Moore \
             form (KRM=3); {detail}. No resonance would contribute a cross section \
             anywhere on the energy grid."
        )));
    }

    Ok(data)
}

/// Parse File 2 (and File 3) records without the evaluable-range check.
fn parse_file2_records(endf_text: &str) -> Result<ResonanceData, EndfParseError> {
    let lines = section_151(endf_text, "2");

    if lines.is_empty() {
        return Err(EndfParseError::MissingSection(
            "No MF=2, MT=151 data found".to_string(),
        ));
    }

    let mut pos = 0;

    // HEAD record: ZA, AWR, 0, 0, NIS, 0
    let head = parse_cont(&lines, &mut pos)?;
    let za = head.c1 as u32;
    let awr = head.c2;
    let nis = checked_count(head.n1, "NIS")?; // number of isotopes (usually 1)

    // ENDF-6 §2.1 requires NIS >= 1 for a valid resonance evaluation. NIS=0
    // would leave the parser with no isotope subsection to read and fall
    // through to a confusing "unconsumed data lines" downstream failure;
    // reject up-front with a clear message.
    if nis == 0 {
        return Err(EndfParseError::UnsupportedFormat(
            "MF=2 NIS=0: no isotopes declared. ENDF-6 §2.1 requires NIS >= 1 \
             for a valid resonance evaluation."
                .into(),
        ));
    }
    // ENDF-6 §2.1: a material with NIS>1 contains multiple isotope subsections,
    // each carrying its own ZAI, ABN, LFW, and NER ranges (e.g. natural-element
    // evaluations such as nat-C with ZAI={6012,6013}). The reference reader
    // OpenScale (File2.cpp:71-87) stores these in a Vec<ResonanceIsotope>
    // tagged with per-isotope ABN, and downstream physics combines the per-
    // isotope cross sections with that ABN weighting. NEREIDS's `ResonanceData`
    // has no per-isotope discriminator and no abundance field; silently
    // flattening multi-isotope subsections into one flat range list would
    // discard ZAI/ABN and produce abundance-blind cross sections. Reject
    // the case explicitly until a proper multi-isotope container is wired
    // through. Single-isotope (NIS=1) ENDF evaluations remain fully supported.
    if nis > 1 {
        return Err(EndfParseError::UnsupportedFormat(format!(
            "MF=2 NIS={nis} > 1 multi-isotope materials are not supported. \
             Each NEREIDS ResonanceData represents a single isotope. \
             For multi-isotope ENDF evaluations, split the material into per-isotope \
             files or use sammy_to_resonance_data_multi from nereids_endf::sammy."
        )));
    }

    let isotope = isotope_from_za(za)?;
    let mut all_ranges = Vec::new();

    for _ in 0..nis {
        // Isotope CONT: ZAI, ABN, 0, LFW, NER, 0
        let iso_cont = parse_cont(&lines, &mut pos)?;
        let _zai = iso_cont.c1 as u32;
        let _abn = iso_cont.c2; // abundance
        let lfw = iso_cont.l2; // fission width flag (LFW=1 → energy-dependent fission widths in URR)
        let ner = checked_count(iso_cont.n1, "NER")?; // number of energy ranges

        for _ in 0..ner {
            // Range CONT: EL, EH, LRU, LRF, NRO, NAPS
            let range_cont = parse_cont(&lines, &mut pos)?;
            let energy_low = range_cont.c1;
            let energy_high = range_cont.c2;
            let lru = range_cont.l1; // 1=resolved, 2=unresolved
            let lrf = range_cont.l2; // resonance formalism

            if lru == 2 {
                // Unresolved resonance region (LRU=2): parsed and SKIPPED.
                //
                // NEREIDS does not compute URR average cross sections — the
                // width-fluctuation correction, a systematically wrong average).
                // We still consume the URR body to keep the line stream aligned
                // so resolved ranges in the same evaluation remain accessible,
                // and tag the range Unresolved (non-evaluable → Skip in the
                // physics layer). Structural guards inside the skip helpers
                // detect malformed records that would otherwise misalign the
                // cursor. Reference: ENDF-6 §2.2.2.
                let nro_urr = range_cont.n1; // energy-dependent scattering radius flag
                let naps_urr = range_cont.n2; // scattering radius calculation flag
                // A TAB1 AP(E) record precedes the URR body when NRO != 0.
                // Store it (converted to fm) on the placeholder: the placeholder
                // preserves everything the header carries even though the range
                // is never evaluated.
                let ap_table_urr = if nro_urr != 0 {
                    let mut tab = parse_tab1(&lines, &mut pos)?;
                    for pt in &mut tab.points {
                        pt.1 *= ENDF_RADIUS_TO_FM;
                    }
                    Some(tab)
                } else {
                    None
                };

                // LRF must be 1 or 2 for the URR (ENDF-6 §2.2.2). Any other
                // value is a malformed URR record; reject it as an unsupported
                // format — consistent with the sibling malformed-record guards
                // (INT range, N1 relations) — rather than heuristically
                // consuming a guessed body and silently dropping the span.
                if lrf != 1 && lrf != 2 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRU=2 (URR) with LRF={lrf}: ENDF-6 §2.2.2 restricts \
                         the unresolved region to LRF=1 or LRF=2"
                    )));
                }

                // LFW=1/LRF=1 uses the "Case B" shared-energy-grid layout
                // (ENDF-6 §2.2.2.1); LFW=1/LRF=2 is byte-identical to LFW=0/LRF=2.
                let (spi_urr, ap_urr_fm) = if lfw != 0 && lrf == 1 {
                    skip_urr_lfw1_lrf1(&lines, &mut pos)?
                } else {
                    skip_urr_range(&lines, &mut pos, lrf)?
                };

                all_ranges.push(ResonanceRange {
                    energy_low,
                    energy_high,
                    resolved: false,
                    formalism: ResonanceFormalism::Unresolved,
                    target_spin: spi_urr,
                    scattering_radius: ap_urr_fm,
                    naps: naps_urr,
                    ap_table: ap_table_urr,
                    l_groups: Vec::new(),
                    r_external: vec![],
                    r_matrix_limited: None,
                });
                continue;
            }

            if lru == 0 {
                // LRU=0: scattering-radius-only range (no resonance parameters).
                // ENDF-6 §2.1: after the range CONT (and optional TAB1 if
                // NRO!=0), a single CONT record follows: [SPI, AP, 0, 0,
                // NLS=0, 0]. We consume the NRO TAB1 if present, then the CONT,
                // and store a non-evaluable placeholder (rather than dropping
                // the range) so the LRU=0 span is named by skip_description and
                // the no-evaluable-content error — a file whose only range is
                // LRU=0 is still rejected (NEREIDS cannot evaluate a
                // resonance-parameter-free stanza, and loading it would yield
                // silent zero cross-sections), but with an accurate message.
                let nro_lru0 = range_cont.n1;
                let naps_lru0 = range_cont.n2;
                let ap_table_lru0 = if nro_lru0 != 0 {
                    // TAB1 AP(E) precedes the SPI/AP CONT; consume and store it.
                    let mut tab = parse_tab1(&lines, &mut pos)?;
                    for pt in &mut tab.points {
                        pt.1 *= ENDF_RADIUS_TO_FM;
                    }
                    Some(tab)
                } else {
                    None
                };
                // CONT: SPI, AP, 0, 0, NLS=0, 0
                // Validate NLS=0: a non-zero NLS in an LRU=0 range is
                // malformed and would cause the parser to look for L-groups that
                // don't exist, misaligning the cursor for subsequent ranges.
                let spi_cont = parse_cont(&lines, &mut pos)?;
                // ENDF-6 §2.1: the SPI/AP CONT is [SPI, AP, 0, 0, NLS=0, 0].
                // Validate that L1 and L2 are both zero — non-zero values
                // indicate a malformed or mis-identified record.
                if spi_cont.l1 != 0 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRU=0 range: L1={} in SPI/AP CONT record must be 0",
                        spi_cont.l1
                    )));
                }
                if spi_cont.l2 != 0 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRU=0 range: L2={} in SPI/AP CONT record must be 0",
                        spi_cont.l2
                    )));
                }
                if spi_cont.n1 != 0 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRU=0 range: NLS={} in SPI/AP CONT record must be 0 \
                         (scattering-radius-only ranges have no L-groups)",
                        spi_cont.n1
                    )));
                }
                if spi_cont.n2 != 0 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRU=0 range: N2={} in SPI/AP CONT record must be 0",
                        spi_cont.n2
                    )));
                }
                all_ranges.push(ResonanceRange {
                    energy_low,
                    energy_high,
                    resolved: false,
                    formalism: ResonanceFormalism::ScatteringRadiusOnly,
                    target_spin: spi_cont.c1,
                    scattering_radius: spi_cont.c2 * ENDF_RADIUS_TO_FM,
                    naps: naps_lru0,
                    ap_table: ap_table_lru0,
                    l_groups: Vec::new(),
                    r_external: vec![],
                    r_matrix_limited: None,
                });
                continue;
            }

            if lru != 1 {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "LRU={} not supported (expected 0=scattering-radius-only, 1=resolved, or 2=unresolved)",
                    lru
                )));
            }

            let nro = range_cont.n1; // energy-dependent scattering radius flag
            let naps = range_cont.n2; // scattering radius calculation flag

            // If NRO != 0, a TAB1 record immediately follows giving AP(E).
            // Parse and store it; scattering_radius_at(E) will interpolate it
            // at each energy point.  Reference: ENDF-6 §2.2.1; SAMMY mlb/mmlb1.f90.
            let ap_table = if nro != 0 {
                let mut tab = parse_tab1(&lines, &mut pos)?;
                // AP(E) y-values are in 10⁻¹² cm; convert to fm.
                for pt in &mut tab.points {
                    pt.1 *= ENDF_RADIUS_TO_FM;
                }
                Some(tab)
            } else {
                None
            };

            // ENDF-6 Formats Manual: LRF values for resolved resonance region
            // LRF=1: Single-Level Breit-Wigner (SLBW)
            // LRF=2: Multi-Level Breit-Wigner (MLBW)
            // LRF=3: Reich-Moore
            // LRF=4: Adler-Adler (deprecated, not supported)
            // LRF=7: R-Matrix Limited (general)
            let formalism = match lrf {
                1 => ResonanceFormalism::SLBW,
                2 => ResonanceFormalism::MLBW,
                3 => ResonanceFormalism::ReichMoore,
                7 => ResonanceFormalism::RMatrixLimited,
                _ => {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "LRF={} not yet supported",
                        lrf
                    )));
                }
            };

            let mut ctx = RangeParseContext {
                lines: &lines,
                pos: &mut pos,
                energy_low,
                energy_high,
                naps,
                ap_table,
            };
            let range = match formalism {
                ResonanceFormalism::MLBW | ResonanceFormalism::SLBW => {
                    parse_bw_range(&mut ctx, formalism)?
                }
                ResonanceFormalism::ReichMoore => parse_reich_moore_range(&mut ctx)?,
                ResonanceFormalism::RMatrixLimited => parse_rmatrix_limited_range(&mut ctx)?,
                ResonanceFormalism::Unresolved => {
                    // Unreachable: Unresolved is only assigned in the LRU=2 branch above.
                    unreachable!("Unresolved formalism should not appear in LRU=1 dispatch");
                }
                ResonanceFormalism::ScatteringRadiusOnly => {
                    // Unreachable: ScatteringRadiusOnly is only assigned in the
                    // LRU=0 branch above.
                    unreachable!(
                        "ScatteringRadiusOnly formalism should not appear in LRU=1 dispatch"
                    );
                }
            };
            all_ranges.push(range);
        }
    }

    if pos < lines.len() {
        return Err(EndfParseError::UnsupportedFormat(
            "Multiple materials detected in MF=2/MT=151: unconsumed data lines \
             remain after parsing the first material. Multi-MAT files are not \
             supported; split the file into single-material ENDF files."
                .to_string(),
        ));
    }

    let mat = lines[0].get(66..70).unwrap_or("").trim();
    if let Some(lrp) = section(endf_text, mat, "1", "451")
        .first()
        .map(|line| parse_endf_int(line, 2))
        .transpose()?
        && lrp == 2
    {
        return Err(EndfParseError::UnsupportedFormat(
            "LRP=2: File 3 holds the whole cross section and its resonance parameters are not \
             to be added to it; only LRP=1 evaluations are supported"
                .into(),
        ));
    }
    Ok(ResonanceData {
        isotope,
        za,
        awr,
        smooth: parse_file3(endf_text, mat)?,
        ranges: all_ranges,
    })
}

fn section_151<'a>(endf_text: &'a str, mf: &str) -> Vec<&'a str> {
    endf_text
        .lines()
        .filter(|line| line.len() >= 75 && line[70..72].trim() == mf && &line[72..75] == "151")
        .collect()
}

/// The lines of section `mt` of file `mf` of material `mat`.
fn section<'a>(endf_text: &'a str, mat: &str, mf: &str, mt: &str) -> Vec<&'a str> {
    endf_text
        .lines()
        .filter(|line| {
            line.len() >= 75
                && line[66..70].trim() == mat
                && line[70..72].trim() == mf
                && line[72..75].trim() == mt
        })
        .collect()
}

/// File 3's cross sections of material `mat` for the reactions
/// [`SmoothCrossSections`] holds, each the one TAB1 record that follows its
/// section's HEAD record.
fn parse_file3(endf_text: &str, mat: &str) -> Result<SmoothCrossSections, EndfParseError> {
    let table = |mt: &str| -> Result<Option<Tab1>, EndfParseError> {
        let lines = section(endf_text, mat, "3", mt);
        if lines.is_empty() {
            return Ok(None);
        }
        let mut pos = 1;
        let table = parse_tab1_with(&lines, &mut pos, true)?;
        if pos < lines.len() {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "MF=3 MT={mt} of MAT {mat} holds records after its table"
            )));
        }
        Ok(Some(table))
    };
    Ok(SmoothCrossSections {
        total: table("1")?,
        elastic: table("2")?,
        capture: table("102")?,
        fission: table("18")?,
    })
}

const FILE32_TOLERANCE: f64 = 1e-5;

/// The covariance, from File 32 (MT=151), of the energy, neutron width and
/// capture width of each resonance of `data` with an energy in `window_ev`
/// that File 32 lists with the same spin and an energy and widths that agree
/// to a relative 1e-5.  A resonance listed in several short-range
/// sections has their contributions summed; resonances outside the window,
/// or that File 32 does not list, are left out.  `None` when the text has no
/// File 32 or it lists no resonance in the window.
///
/// # Errors
/// [`EndfParseError::InvalidNumber`] for an empty window or one with a NaN
/// bound.  [`EndfParseError::UnsupportedFormat`] for a File 32 that is not
/// of the isotope `data` holds, or not of its resolved Reich–Moore ranges in
/// the general form (LCOMP=1) over each resonance's energy, neutron width and
/// capture width (MPAR=3), or that gives a scattering-radius covariance
/// (NRO≠0 or ISR≠0) or long-range sections, or an unresolved range other
/// than a last one of LRF=1, holds a number that is not finite or records
/// past its last range, or lists a resonance in the window
/// that matches no Reich–Moore resonance of `data`, more than one, or one its
/// block already lists.  The record errors of [`parse_endf_file2`].
///
/// Reference: ENDF-6 Formats Manual (ENDF-102, 2023), §32.2.2.2 for the
/// layout and §32.3 for the correspondence with File 2 and the summing of a
/// resonance listed in several blocks.
pub fn parse_resonance_covariance(
    endf_text: &str,
    data: &ResonanceData,
    window_ev: std::ops::RangeInclusive<f64>,
) -> Result<Option<ResonanceCovariance>, EndfParseError> {
    if window_ev.is_empty() {
        return Err(EndfParseError::InvalidNumber(format!(
            "the energy window {window_ev:?} must run from a lower to a higher energy"
        )));
    }
    let lines = section_151(endf_text, "32");
    if lines.is_empty() {
        return Ok(None);
    }
    let known: Vec<([usize; 3], &Resonance)> = data
        .ranges
        .iter()
        .enumerate()
        .filter(|(_, range)| range.formalism == ResonanceFormalism::ReichMoore)
        .flat_map(|(a, range)| {
            range
                .l_groups
                .iter()
                .enumerate()
                .flat_map(move |(b, group)| {
                    let resonances = group.resonances.iter().enumerate();
                    resonances.map(move |(c, r)| ([a, b, c], r))
                })
        })
        .collect();
    let agree = |a: [f64; 3], b: [f64; 3]| {
        a.iter()
            .zip(b)
            .all(|(a, b)| (a - b).abs() <= FILE32_TOLERANCE * a.abs())
    };
    let pos = &mut 0;
    let (head, isotope) = (parse_cont(&lines, pos)?, parse_cont(&lines, pos)?);
    if [head.c1, isotope.c1] != [f64::from(data.za); 2] || head.n1 != 1 {
        return Err(EndfParseError::UnsupportedFormat(format!(
            "File 32 is read for the one isotope `data` holds, ZA={}; got ZA={}, ZAI={}, NIS={}",
            data.za, head.c1, isotope.c1, head.n1
        )));
    }
    let ranges = checked_count(isotope.n1, "NER")?;
    let (mut resonances, mut entries) = (Vec::new(), Vec::new());
    for r in 0..ranges {
        let range = parse_cont(&lines, pos)?;
        if range.l1 == 2 && range.l2 == 1 && r + 1 == ranges {
            let groups = checked_count(parse_cont(&lines, pos)?.n1, "NLS")?;
            for _ in 0..=groups {
                let list = parse_cont(&lines, pos)?;
                parse_list_values(&lines, pos, checked_count(list.n1, "NPL")?)?;
            }
            break;
        }
        let (form, sections) = (parse_cont(&lines, pos)?, parse_cont(&lines, pos)?);
        if range.l1 != 1
            || range.l2 != 3
            || range.n1 != 0
            || form.l2 != 1
            || form.n2 != 0
            || sections.n2 != 0
        {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "File 32 is read for one isotope's resolved Reich-Moore ranges in the \
                 general form, without a scattering-radius covariance or long-range \
                 sections, followed by at most one unresolved range (LRF=1); got LRU={}, \
                 LRF={}, NRO={}, LCOMP={}, ISR={}, NLRS={}",
                range.l1, range.l2, range.n1, form.l2, form.n2, sections.n2
            )));
        }
        for _ in 0..checked_count(sections.n1, "NSRS")? {
            let list = parse_cont(&lines, pos)?;
            let count = checked_count(list.n2, "NRB")?;
            let n = 3 * count;
            let length = 6 * count + n * (n + 1) / 2;
            if list.l1 != 3 || usize::try_from(list.n1) != Ok(length) {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "File 32 is read over each resonance's energy, neutron width and capture \
                     width (MPAR=3, NPL=6*NRB+NVS); got MPAR={}, NPL={}, NRB={count}",
                    list.l1, list.n1
                )));
            }
            let values = parse_list_values(&lines, pos, length)?;
            if !values.iter().all(|v| v.is_finite()) {
                return Err(EndfParseError::UnsupportedFormat(
                    "File 32 holds a number that is not finite".into(),
                ));
            }
            let (table, packed) = values.split_at(6 * count);
            let mut slots = Vec::new();
            for (k, row) in table.chunks(6).enumerate() {
                let (energy, spin) = (row[0], row[1]);
                let mut same = known.iter().filter(|(_, r)| {
                    r.j == spin && agree([r.energy, r.gn, r.gg], [row[0], row[2], row[3]])
                });
                let found = match (same.next(), same.next()) {
                    (Some(&(at, r)), None) => Some((at, r.energy)),
                    _ => None,
                };
                if !window_ev.contains(&found.map_or(energy, |(_, e)| e)) {
                    continue;
                }
                let found = found.map(|(at, _)| at);
                let Some(at) = found.filter(|at| slots.iter().all(|&(_, s)| resonances[s] != *at))
                else {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "File 32 lists a resonance at {energy} eV of spin {spin} that matches no \
                         Reich-Moore resonance of File 2 of that spin to a relative \
                         {FILE32_TOLERANCE:e} in energy and widths, more than one, or one its \
                         block already lists"
                    )));
                };
                if !resonances.contains(&at) {
                    resonances.push(at);
                }
                slots.push((k, resonances.iter().position(|&r| r == at).expect("listed")));
            }
            for &(k, a) in &slots {
                for &(l, b) in &slots {
                    for (x, y) in (0..3).flat_map(|x| (0..3).map(move |y| (x, y))) {
                        let (u, v) = ((3 * k + x).min(3 * l + y), (3 * k + x).max(3 * l + y));
                        entries.push((3 * a + x, 3 * b + y, packed[u * n - u * (u + 1) / 2 + v]));
                    }
                }
            }
        }
    }
    if *pos != lines.len() {
        return Err(EndfParseError::UnsupportedFormat(
            "File 32 holds records past its last range".into(),
        ));
    }
    let n = 3 * resonances.len();
    if n == 0 {
        return Ok(None);
    }
    let mut covariance = vec![0.0; n * n];
    for (x, y, value) in entries {
        covariance[x * n + y] += value;
    }
    Ok(Some(ResonanceCovariance {
        resonances,
        covariance: covariance.into(),
    }))
}

/// Shared context for ENDF range parsers.
///
/// Groups the file-position state (`lines`, `pos`) with the fields from the
/// range CONT record that every range parser needs, eliminating long argument
/// lists.
struct RangeParseContext<'a> {
    lines: &'a [&'a str],
    pos: &'a mut usize,
    energy_low: f64,
    energy_high: f64,
    naps: i32,
    ap_table: Option<Tab1>,
}

/// Parse a Breit-Wigner (SLBW or MLBW) resolved resonance range.
///
/// ENDF-6 File 2, LRF=1 (SLBW) / LRF=2 (MLBW):
/// - CONT: SPI, AP, 0, 0, NLS, 0
/// - For each L-value:
///   - CONT: AWRI, 0.0, L, 0, 6*NRS, NRS
///   - LIST: NRS resonances, each 6 values: ER, AJ, GT, GN, GG, GF
///
/// Reference: ENDF-6 Formats Manual Section 2.2.1.1
fn parse_bw_range(
    ctx: &mut RangeParseContext<'_>,
    formalism: ResonanceFormalism,
) -> Result<ResonanceRange, EndfParseError> {
    // CONT: SPI, AP, 0, 0, NLS, 0
    // ENDF AP is in 10⁻¹² cm; convert to fm (×10).
    let cont = parse_cont(ctx.lines, ctx.pos)?;
    let target_spin = cont.c1;
    let scattering_radius = cont.c2 * ENDF_RADIUS_TO_FM;
    let nls = checked_count(cont.n1, "NLS")?; // number of L-values

    // Zero NLS means a resolved range with no L-groups: out of spec
    // (ENDF-6 §2.2.1.1 requires NLS >= 1; NLS=0 belongs to LRU=0
    // scattering-radius-only stanzas) and silently inert if accepted —
    // the range would count as evaluable yet contribute zero cross-section
    // everywhere. Mirrors the NJS=0 guards in the URR skip paths.
    if nls == 0 {
        return Err(EndfParseError::UnsupportedFormat(format!(
            "{formalism:?} range: NLS=0 (ENDF-6 §2.2.1.1 requires at least one \
             L-group in a resolved range; NLS=0 is reserved for LRU=0)"
        )));
    }

    let mut l_groups = Vec::with_capacity(nls);

    for _ in 0..nls {
        // CONT: AWRI, QX, L, LRX, 6*NRS, NRS
        let l_cont = parse_cont(ctx.lines, ctx.pos)?;
        let awr_l = l_cont.c1;
        let qx = l_cont.c2; // Q-value for competitive width (eV)
        // Validate L is non-negative: negative L1 wraps to a huge u32.
        if l_cont.l1 < 0 {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "BW range: negative L={}",
                l_cont.l1
            )));
        }
        let l_val = l_cont.l1 as u32;
        let lrx = l_cont.l2; // competitive width flag
        let n1 = checked_count(l_cont.n1, "N1")?; // should be 6*NRS
        let nrs = checked_count(l_cont.n2, "NRS")?; // number of resonances

        // Validate N1 == 6*NRS: a mismatch means the record is malformed
        // and reading N1 values would over-/under-consume lines.
        if n1 != 6 * nrs {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "BW range L={l_val}: N1={n1} != 6*NRS={} (NRS={nrs})",
                6 * nrs
            )));
        }

        let mut resonances = Vec::with_capacity(nrs);

        // Each resonance is 6 values on one line (or spanning lines).
        // In ENDF format, LIST records pack 6 values per line.
        let total_values = nrs * 6;
        let values = parse_list_values(ctx.lines, ctx.pos, total_values)?;

        for i in 0..nrs {
            let base = i * 6;
            resonances.push(Resonance {
                energy: values[base],  // ER
                j: values[base + 1],   // AJ
                gn: values[base + 3],  // GN (neutron width)
                gg: values[base + 4],  // GG (gamma width)
                gfa: values[base + 5], // GF (fission width)
                gfb: 0.0,              // Not used in BW
                                       // GT (values[base+2]) is not stored: with LRX = 0 it is
                                       // GN + GG + GF, and a range with LRX ≠ 0 is not evaluated.
            });
        }

        l_groups.push(LGroup {
            l: l_val,
            awr: awr_l,
            apl: 0.0, // Not in BW format
            qx,
            lrx,
            resonances,
        });
    }

    // A resolved range whose L-groups are all empty (every NRS=0) is accepted
    // by reference readers (SAMMY, OpenScale), so it is NOT a hard error here.
    // NEREIDS treats it as a placeholder, although SAMMY gives its channels
    // hard-sphere scattering: `is_evaluable()` classifies it non-evaluable,
    // a warn-and-skip placeholder like
    // LRF=7/LRU=2/LRU=0, so a mixed file keeps its populated ranges, and a
    // file where this is the only range fails the no-evaluable-content guard
    // (fail-closed beats silently reporting zero as physics).

    Ok(ResonanceRange {
        energy_low: ctx.energy_low,
        energy_high: ctx.energy_high,
        resolved: true,
        formalism,
        target_spin,
        scattering_radius,
        naps: ctx.naps,
        ap_table: ctx.ap_table.take(),
        l_groups,
        r_external: vec![],
        r_matrix_limited: None,
    })
}

/// Parse a Reich-Moore resolved resonance range.
///
/// ENDF-6 File 2, LRF=2:
/// ENDF-6 File 2, LRF=3 (Reich-Moore):
/// - CONT: SPI, AP, 0, 0, NLS, 0
/// - For each L-value:
///   - CONT: AWRI, APL, L, 0, 6*NRS, NRS
///   - LIST: NRS resonances, each 6 values: ER, AJ, GN, GG, GFA, GFB
///
/// Reference: ENDF-6 Formats Manual Section 2.2.1.3
/// Reference: SAMMY manual Section 2 (R-matrix theory)
fn parse_reich_moore_range(
    ctx: &mut RangeParseContext<'_>,
) -> Result<ResonanceRange, EndfParseError> {
    // CONT: SPI, AP, 0, 0, NLS, 0
    // ENDF AP is in 10⁻¹² cm; convert to fm (×10).
    let cont = parse_cont(ctx.lines, ctx.pos)?;
    let target_spin = cont.c1;
    let scattering_radius = cont.c2 * ENDF_RADIUS_TO_FM;
    let nls = checked_count(cont.n1, "NLS")?; // number of L-values

    // Zero NLS means a resolved range with no L-groups: out of spec
    // (ENDF-6 §2.2.1.1 requires NLS >= 1; NLS=0 belongs to LRU=0
    // scattering-radius-only stanzas) and silently inert if accepted —
    // the range would count as evaluable yet contribute zero cross-section
    // everywhere. Mirrors the NJS=0 guards in the URR skip paths.
    if nls == 0 {
        return Err(EndfParseError::UnsupportedFormat(
            "Reich-Moore range: NLS=0 (ENDF-6 §2.2.1.1 requires at least one \
             L-group in a resolved range; NLS=0 is reserved for LRU=0)"
                .to_string(),
        ));
    }

    let mut l_groups = Vec::with_capacity(nls);

    for _ in 0..nls {
        // CONT: AWRI, APL, L, 0, 6*NRS, NRS
        let l_cont = parse_cont(ctx.lines, ctx.pos)?;
        let awr_l = l_cont.c1;
        let apl = l_cont.c2 * ENDF_RADIUS_TO_FM; // L-dependent scattering radius
        // Validate L is non-negative: negative L1 wraps to a huge u32.
        if l_cont.l1 < 0 {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "Reich-Moore range: negative L={}",
                l_cont.l1
            )));
        }
        let l_val = l_cont.l1 as u32;
        let n1 = checked_count(l_cont.n1, "N1")?; // should be 6*NRS
        let nrs = checked_count(l_cont.n2, "NRS")?; // number of resonances

        // Validate N1 == 6*NRS: a mismatch means the record is malformed
        // and reading N1 values would over-/under-consume lines.
        if n1 != 6 * nrs {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "Reich-Moore range L={l_val}: N1={n1} != 6*NRS={} (NRS={nrs})",
                6 * nrs
            )));
        }

        let mut resonances = Vec::with_capacity(nrs);

        // Each resonance is 6 values: ER, AJ, GN, GG, GFA, GFB
        let total_values = nrs * 6;
        let values = parse_list_values(ctx.lines, ctx.pos, total_values)?;

        for i in 0..nrs {
            let base = i * 6;
            resonances.push(Resonance {
                energy: values[base],  // ER (eV)
                j: values[base + 1],   // AJ (total J)
                gn: values[base + 2],  // GN (neutron width, eV)
                gg: values[base + 3],  // GG (gamma width, eV)
                gfa: values[base + 4], // GFA (fission width 1, eV)
                gfb: values[base + 5], // GFB (fission width 2, eV)
            });
        }

        l_groups.push(LGroup {
            l: l_val,
            awr: awr_l,
            apl,
            qx: 0.0, // Not used in Reich-Moore
            lrx: 0,  // Not used in Reich-Moore
            resonances,
        });
    }

    // A resolved range whose L-groups are all empty (every NRS=0) is accepted
    // by reference readers (SAMMY, OpenScale), so it is NOT a hard error here.
    // NEREIDS treats it as a placeholder, although SAMMY gives its channels
    // hard-sphere scattering: `is_evaluable()` classifies it non-evaluable,
    // a warn-and-skip placeholder like
    // LRF=7/LRU=2/LRU=0, so a mixed file keeps its populated ranges, and a
    // file where this is the only range fails the no-evaluable-content guard
    // (fail-closed beats silently reporting zero as physics).

    Ok(ResonanceRange {
        energy_low: ctx.energy_low,
        energy_high: ctx.energy_high,
        resolved: true,
        formalism: ResonanceFormalism::ReichMoore,
        target_spin,
        scattering_radius,
        naps: ctx.naps,
        ap_table: ctx.ap_table.take(),
        l_groups,
        r_external: vec![],
        r_matrix_limited: None,
    })
}

/// Parse an R-Matrix Limited (LRF=7) resolved range in the Reich-Moore form.
///
/// Accepts KRM=3 with widths in eV (IFG=0), KRL=0, hard-sphere phase shifts
/// (KPS=0), constant AP (NRO=0), and background R-functions with LBK=0 or 2.
///
/// # Errors
/// Returns [`EndfParseError::UnsupportedFormat`] for a flag outside that set, a
/// fission pair, a shift factor on a channel that closes inside the range, and
/// any record whose counts or values are inconsistent.
///
/// # References
/// ENDF-6 Formats Manual, ENDF-102 (2023), §2.2.1.6 and Appendix D.1.7.
fn parse_rmatrix_limited_range(
    ctx: &mut RangeParseContext<'_>,
) -> Result<ResonanceRange, EndfParseError> {
    if ctx.ap_table.is_some() {
        return Err(rml_error(
            "NRO=1 (energy-dependent AP) is not supported; LRF=7 gives radii per channel",
        ));
    }

    // CONT: [SPI, AP, IFG, KRM, NJS, KRL]; AP in 10⁻¹² cm.
    let cont = parse_cont(ctx.lines, ctx.pos)?;
    let target_spin = cont.c1;
    let scattering_radius = cont.c2 * ENDF_RADIUS_TO_FM;
    if cont.l1 != 0 {
        return Err(rml_error(format!(
            "IFG={} (reduced-width amplitudes) is not supported; only IFG=0 (widths in eV)",
            cont.l1
        )));
    }
    if cont.l2 != 3 {
        return Err(rml_error(format!(
            "KRM={} is not supported; only KRM=3 (Reich-Moore) is evaluated",
            cont.l2
        )));
    }
    let njs = checked_count(cont.n1, "NJS")?;
    if cont.n2 != 0 {
        return Err(rml_error(format!(
            "KRL={} (relativistic kinematics) is not supported; only KRL=0",
            cont.n2
        )));
    }

    // LIST: [0, 0, NPP, 0, 12*NPP, 2*NPP] with 12 values per pair.
    let pair_head = parse_cont(ctx.lines, ctx.pos)?;
    let npp = checked_count(pair_head.l1, "NPP")?;
    let npl = checked_count(pair_head.n1, "NPL")?;
    if npl != 12 * npp {
        return Err(rml_error(format!(
            "particle-pair LIST: NPL={npl} != 12*NPP={}",
            12 * npp
        )));
    }
    let pair_values = parse_list_values(ctx.lines, ctx.pos, npl)?;
    let pairs = pair_values
        .chunks_exact(12)
        .enumerate()
        .map(|(i, values)| parse_particle_pair(i + 1, values))
        .collect::<Result<Vec<_>, _>>()?;

    let photon_pair = sole_index(&pairs, |p| p.mass_a == 0.0 && p.mt == 102).ok_or_else(|| {
        rml_error("KRM=3 needs exactly one photon pair (MA=0, MT=102) for the eliminated channel")
    })?;
    let incident_pair = sole_index(&pairs, |p| p.mt == 2)
        .ok_or_else(|| rml_error("exactly one particle pair must have MT=2 (neutron + target)"))?;
    let incident = &pairs[incident_pair];
    if incident.mass_a != 1.0 || incident.charge_a != 0.0 || incident.q_value != 0.0 {
        return Err(rml_error(format!(
            "the MT=2 pair must be a neutron with Q=0 (MA={}, ZA={}, Q={})",
            incident.mass_a, incident.charge_a, incident.q_value
        )));
    }
    for (i, pair) in pairs.iter().enumerate() {
        let threshold = pair.threshold(incident);
        if pair.shift && threshold > ctx.energy_low {
            return Err(rml_error(format!(
                "particle pair {}: SHF=1 on a channel with threshold {threshold:.6e} eV inside \
                 the range is not supported (closed-channel shift factor)",
                i + 1
            )));
        }
    }

    let mut spin_groups = Vec::with_capacity(njs);
    for group in 1..=njs {
        spin_groups.push(parse_rml_spin_group(
            ctx,
            group,
            &pairs,
            photon_pair,
            incident_pair,
        )?);
    }

    Ok(ResonanceRange {
        energy_low: ctx.energy_low,
        energy_high: ctx.energy_high,
        resolved: true,
        formalism: ResonanceFormalism::RMatrixLimited,
        target_spin,
        scattering_radius,
        naps: ctx.naps,
        ap_table: None,
        l_groups: Vec::new(),
        r_external: vec![],
        r_matrix_limited: Some(RmlRange {
            pairs,
            incident_pair,
            photon_pair,
            spin_groups,
        }),
    })
}

/// Error for an LRF=7 record NEREIDS does not accept.
fn rml_error(message: impl std::fmt::Display) -> EndfParseError {
    EndfParseError::UnsupportedFormat(format!("LRF=7 {message}"))
}

/// Index of the only item matching `predicate`, or `None` for zero or several.
fn sole_index<T>(items: &[T], predicate: impl Fn(&T) -> bool) -> Option<usize> {
    let mut found = None;
    for (i, item) in items.iter().enumerate() {
        if predicate(item) {
            if found.is_some() {
                return None;
            }
            found = Some(i);
        }
    }
    found
}

/// Read an integer stored in an ENDF float field; fractions and non-finite values fail.
fn rml_integer(value: f64, field: &str, location: &str) -> Result<i64, EndfParseError> {
    if !value.is_finite() || value.fract() != 0.0 || value.abs() > f64::from(MAX_ENDF_COUNT) {
        return Err(rml_error(format!(
            "{location}: {field}={value} is not an integer"
        )));
    }
    Ok(value as i64)
}

/// Parse one 12-value particle-pair entry, `number` counting from 1.
///
/// PNT=0 takes the ENDF-6 default: no penetrability for MT=102 and MT=19,
/// penetrability for every other MT.
fn parse_particle_pair(number: usize, values: &[f64]) -> Result<ParticlePair, EndfParseError> {
    let location = format!("particle pair {number}");
    let &[ma, mb, za, zb, ia, ib, q, pnt, shf, mt, _pa, _pb] = values else {
        return Err(rml_error(format!("{location}: expected 12 values")));
    };
    for (field, value) in [
        ("MA", ma),
        ("MB", mb),
        ("ZA", za),
        ("ZB", zb),
        ("IA", ia),
        ("IB", ib),
        ("Q", q),
    ] {
        if !value.is_finite() {
            return Err(rml_error(format!(
                "{location}: {field}={value} is not finite"
            )));
        }
    }
    let pnt = rml_integer(pnt, "PNT", &location)?;
    let shf = rml_integer(shf, "SHF", &location)?;
    let mt = rml_integer(mt, "MT", &location)?;
    if mt <= 0 {
        return Err(rml_error(format!(
            "{location}: MT={mt} is not a reaction type"
        )));
    }
    let mt = mt as u32;
    let penetrability = match pnt {
        1 => true,
        -1 => false,
        0 => !(mt == 102 || mt == 19),
        _ => {
            return Err(rml_error(format!(
                "{location}: PNT={pnt} is not defined (ENDF-6 allows -1, 0 and 1)"
            )));
        }
    };
    let shift = match shf {
        0 => false,
        1 => true,
        _ => {
            return Err(rml_error(format!(
                "{location}: SHF={shf} is not defined (ENDF-6 allows 0 and 1)"
            )));
        }
    };
    if ma == 0.0 && penetrability {
        return Err(rml_error(format!(
            "{location}: a photon pair (MA=0) cannot calculate the penetrability (PNT={pnt}, MT={mt})"
        )));
    }
    if shift && !penetrability {
        return Err(rml_error(format!(
            "{location}: SHF=1 needs the penetrability calculated (PNT={pnt}, MT={mt})"
        )));
    }
    if penetrability {
        let reduced_mass = ma * mb / (ma + mb);
        if !(ma > 0.0 && mb > 0.0 && reduced_mass.is_finite() && reduced_mass > 0.0) {
            return Err(rml_error(format!(
                "{location}: a penetrating pair needs positive masses (MA={ma}, MB={mb})"
            )));
        }
    }
    if matches!(mt, 18..=21 | 38) {
        return Err(rml_error(format!(
            "{location}: fission pair (MT={mt}) is not supported"
        )));
    }
    Ok(ParticlePair {
        mass_a: ma,
        mass_b: mb,
        charge_a: za,
        charge_b: zb,
        spin_a: ia,
        spin_b: ib,
        q_value: q,
        mt,
        penetrability,
        shift,
    })
}

/// Parse one spin group: channel LIST, resonance LIST and KBK background records.
fn parse_rml_spin_group(
    ctx: &mut RangeParseContext<'_>,
    group: usize,
    pairs: &[ParticlePair],
    photon_pair: usize,
    incident_pair: usize,
) -> Result<RmlSpinGroup, EndfParseError> {
    let location = format!("spin group {group}");

    // LIST: [AJ, PJ, KBK, KPS, 6*NCH, NCH] with 6 values per channel.
    let head = parse_cont(ctx.lines, ctx.pos)?;
    let (aj, pj) = (head.c1, head.c2);
    if head.l2 != 0 {
        return Err(rml_error(format!(
            "{location}: KPS={} (tabulated phase shifts) is not supported",
            head.l2
        )));
    }
    let kbk = checked_count(head.l1, "KBK")?;
    let nch = checked_count(head.n2, "NCH")?;
    let npl = checked_count(head.n1, "NPL")?;
    if npl != 6 * nch {
        return Err(rml_error(format!(
            "{location}: channel LIST NPL={npl} != 6*NCH={}",
            6 * nch
        )));
    }
    if !aj.is_finite() || (2.0 * aj).fract() != 0.0 {
        return Err(rml_error(format!(
            "{location}: AJ={aj} is not a multiple of 1/2"
        )));
    }
    let parity = if aj < 0.0 || (aj == 0.0 && pj < 0.0) {
        -1
    } else {
        1
    };

    let channel_values = parse_list_values(ctx.lines, ctx.pos, npl)?;
    let mut channels = Vec::with_capacity(nch);
    for (c, v) in channel_values.chunks_exact(6).enumerate() {
        let channel_location = format!("{location} channel {}", c + 1);
        let ppi = rml_integer(v[0], "PPI", &channel_location)?;
        if ppi < 1 || ppi as usize > pairs.len() {
            return Err(rml_error(format!(
                "{channel_location}: PPI={ppi} is out of range 1..={}",
                pairs.len()
            )));
        }
        let l = rml_integer(v[1], "L", &channel_location)?;
        if l < 0 {
            return Err(rml_error(format!("{channel_location}: L={l} is negative")));
        }
        let (sch, bnd, ape, apt) = (v[2], v[3], v[4], v[5]);
        if ![sch, bnd, ape, apt].iter().all(|x| x.is_finite()) || ape < 0.0 || apt < 0.0 {
            return Err(rml_error(format!(
                "{channel_location}: SCH={sch}, BND={bnd}, APE={ape}, APT={apt} must be finite \
                 with non-negative radii"
            )));
        }
        let pair = ppi as usize - 1;
        if pairs[pair].penetrability && apt == 0.0 {
            return Err(rml_error(format!(
                "{channel_location}: APT=0 for a pair that calculates the penetrability"
            )));
        }
        channels.push(RmlChannel {
            pair,
            l: l as u32,
            channel_spin: sch,
            boundary: bnd,
            radius_phase: ape * ENDF_RADIUS_TO_FM,
            radius_true: apt * ENDF_RADIUS_TO_FM,
        });
    }
    if channels.iter().filter(|c| c.pair == photon_pair).count() != 1 {
        return Err(rml_error(format!(
            "{location}: needs exactly one channel of the photon pair (the eliminated channel)"
        )));
    }
    if !channels.iter().any(|c| c.pair == incident_pair) {
        return Err(rml_error(format!(
            "{location}: has no neutron entrance channel"
        )));
    }

    // LIST: [0, 0, 0, NRS, 6*NX, NX]; each resonance holds ER and NCH widths on
    // NCH/6+1 lines (ENDF-102 §2.2.1.6), and NRS=0 still carries one line.
    let res_head = parse_cont(ctx.lines, ctx.pos)?;
    let nrs = checked_count(res_head.l2, "NRS")?;
    let nx = checked_count(res_head.n2, "NX")?;
    let res_npl = checked_count(res_head.n1, "NPL")?;
    let lines_per_resonance = nch / 6 + 1;
    let expected_nx = if nrs == 0 {
        1
    } else {
        nrs * lines_per_resonance
    };
    if nx != expected_nx || res_npl != 6 * nx {
        return Err(rml_error(format!(
            "{location}: resonance LIST NX={nx}, NPL={res_npl}; NRS={nrs} and NCH={nch} need \
             NX={expected_nx} and NPL={}",
            6 * expected_nx
        )));
    }
    let res_values = parse_list_values(ctx.lines, ctx.pos, res_npl)?;
    let mut resonances = Vec::with_capacity(nrs);
    if nrs > 0 {
        for row in res_values.chunks_exact(6 * lines_per_resonance) {
            let energy = row[0];
            let widths = row[1..=nch].to_vec();
            if !energy.is_finite() || !widths.iter().all(|w| w.is_finite()) {
                return Err(rml_error(format!(
                    "{location}: resonance at ER={energy} has a non-finite value"
                )));
            }
            // A width at the channel threshold has P = 0 and no finite amplitude.
            for (channel, &width) in channels.iter().zip(&widths) {
                let pair = &pairs[channel.pair];
                if pair.penetrability
                    && width != 0.0
                    && energy == pair.threshold(&pairs[incident_pair])
                {
                    return Err(rml_error(format!(
                        "{location}: resonance at ER={energy} eV sits on the threshold of a \
                         penetrating channel with width {width} eV"
                    )));
                }
            }
            resonances.push(RmlResonance { energy, widths });
        }
    }

    // KBK records: [0, 0, LCH, LBK, 0, 0] then the LBK body.
    let mut background = Vec::new();
    let mut seen = vec![false; nch];
    for _ in 0..kbk {
        let record = parse_cont(ctx.lines, ctx.pos)?;
        let lch = checked_count(record.l1, "LCH")?;
        if lch == 0 || lch > nch {
            return Err(rml_error(format!(
                "{location}: background LCH={lch} is out of range 1..={nch}"
            )));
        }
        let channel = lch - 1;
        if seen[channel] {
            return Err(rml_error(format!(
                "{location}: background given twice for channel {lch}"
            )));
        }
        seen[channel] = true;
        if channels[channel].pair == photon_pair {
            return Err(rml_error(format!(
                "{location}: background on the eliminated photon channel {lch}"
            )));
        }
        match record.l2 {
            0 => {}
            2 => {
                // LIST: [ED, EU, 0, 0, 5, 0 / R0, R1, R2, S0, S1].
                let list = parse_cont(ctx.lines, ctx.pos)?;
                if list.n1 != 5 {
                    return Err(rml_error(format!(
                        "{location}: LBK=2 LIST has NPL={}, expected 5",
                        list.n1
                    )));
                }
                let v = parse_list_values(ctx.lines, ctx.pos, 5)?;
                let (e_down, e_up) = (list.c1, list.c2);
                let valid = e_down < ctx.energy_low
                    && e_up > ctx.energy_high
                    && v.iter().all(|x| x.is_finite());
                if !valid {
                    return Err(rml_error(format!(
                        "{location}: LBK=2 needs finite terms and ED < EL, EU > EH \
                         (ED={e_down}, EU={e_up})"
                    )));
                }
                background.push(RmlBackground {
                    channel,
                    e_down,
                    e_up,
                    r_con: v[0],
                    r_lin: v[1],
                    r_quad: v[2],
                    s_con: v[3],
                    s_lin: v[4],
                });
            }
            1 | 3 => {
                return Err(rml_error(format!(
                    "{location}: LBK={} background is not supported (only LBK=0 and 2)",
                    record.l2
                )));
            }
            lbk => {
                return Err(rml_error(format!("{location}: LBK={lbk} is not defined")));
            }
        }
    }

    Ok(RmlSpinGroup {
        j: aj.abs(),
        parity,
        channels,
        resonances,
        background,
    })
}

/// Maximum sane ENDF count value.
///
/// ENDF files in practice never contain more than ~100k resonances per section.
/// Accepting `i32::MAX` would cause enormous allocations (gigabytes) on
/// malformed files.  This cap is generous enough for any real evaluation while
/// protecting against allocation bombs.
const MAX_ENDF_COUNT: i32 = 1_000_000;

/// Validate that an ENDF integer count is non-negative and return as `usize`.
///
/// Malformed records can contain negative counts which, if cast directly to
/// `usize`, wrap to huge values and cause OOM panics in `Vec::with_capacity`
/// or `parse_list_values`.
fn checked_count(value: i32, label: &str) -> Result<usize, EndfParseError> {
    if value < 0 {
        return Err(EndfParseError::UnsupportedFormat(format!(
            "Negative ENDF count: {label}={value}"
        )));
    }
    if value > MAX_ENDF_COUNT {
        return Err(EndfParseError::UnsupportedFormat(format!(
            "ENDF count too large: {label}={value} (maximum {MAX_ENDF_COUNT})"
        )));
    }
    Ok(value as usize)
}

// ---------------------------------------------------------------------------
// Low-level ENDF line parsing helpers
// ---------------------------------------------------------------------------

/// Parsed CONT (control) record with 2 floats, 4 integers.
struct ContRecord {
    c1: f64,
    c2: f64,
    l1: i32,
    l2: i32,
    n1: i32,
    n2: i32,
}

/// Parse a CONT record from the current line.
fn parse_cont(lines: &[&str], pos: &mut usize) -> Result<ContRecord, EndfParseError> {
    if *pos >= lines.len() {
        return Err(EndfParseError::UnexpectedEof(
            "Expected CONT record but reached end of data".to_string(),
        ));
    }
    let line = lines[*pos];
    *pos += 1;

    Ok(ContRecord {
        c1: parse_endf_float(line, 0)?,
        c2: parse_endf_float(line, 1)?,
        l1: parse_endf_int(line, 2)?,
        l2: parse_endf_int(line, 3)?,
        n1: parse_endf_int(line, 4)?,
        n2: parse_endf_int(line, 5)?,
    })
}

/// Parse a LIST of floating-point values spanning multiple lines.
///
/// ENDF packs 6 values per line. We read ceil(n/6) lines.
fn parse_list_values(
    lines: &[&str],
    pos: &mut usize,
    n_values: usize,
) -> Result<Vec<f64>, EndfParseError> {
    let mut values = Vec::with_capacity(n_values.min(6 * lines.len().saturating_sub(*pos)));
    let n_lines = n_values.div_ceil(6);

    for _ in 0..n_lines {
        if *pos >= lines.len() {
            return Err(EndfParseError::UnexpectedEof(
                "Expected LIST data but reached end".to_string(),
            ));
        }
        let line = lines[*pos];
        *pos += 1;

        let remaining = n_values - values.len();
        let fields_on_line = remaining.min(6);

        for field in 0..fields_on_line {
            values.push(parse_endf_float(line, field)?);
        }
    }

    Ok(values)
}

/// Parse a floating-point value from an 11-character ENDF field.
///
/// ENDF uses Fortran-style floats that may omit 'E', e.g.:
/// - " 1.234567+2" means 1.234567e+2
/// - "-3.456789-1" means -3.456789e-1
/// - " 0.000000+0" means 0.0
fn parse_endf_float(line: &str, field_index: usize) -> Result<f64, EndfParseError> {
    let start = field_index * 11;
    let end = start + 11;

    if line.len() < end {
        // Short line — treat as zero.
        return Ok(0.0);
    }

    let field = &line[start..end];
    let trimmed = field.trim();

    if trimmed.is_empty() {
        return Ok(0.0);
    }

    // Try standard Rust float parsing first.
    if let Ok(v) = trimmed.parse::<f64>() {
        return Ok(v);
    }

    // Handle Fortran-style: "1.234567+2" or "-3.456789-1"
    // Look for +/- that is NOT the first character and NOT preceded by 'e'/'E'/'d'/'D'.
    let bytes = trimmed.as_bytes();
    for i in 1..bytes.len() {
        if (bytes[i] == b'+' || bytes[i] == b'-')
            && bytes[i - 1] != b'e'
            && bytes[i - 1] != b'E'
            && bytes[i - 1] != b'd'
            && bytes[i - 1] != b'D'
            && bytes[i - 1] != b'+'
            && bytes[i - 1] != b'-'
        {
            let mantissa = &trimmed[..i];
            let exp_slice = &trimmed[i..];
            // Strip spaces from the exponent only when present (some ENDF files
            // write "+ 4" not "+4").  Avoid allocation on the common path.
            let with_e = if exp_slice.contains(' ') {
                let exponent: String = exp_slice.chars().filter(|c| !c.is_whitespace()).collect();
                format!("{}E{}", mantissa, exponent)
            } else {
                format!("{}E{}", mantissa, exp_slice)
            };
            if let Ok(v) = with_e.parse::<f64>() {
                return Ok(v);
            }
        }
    }

    Err(EndfParseError::InvalidNumber(format!(
        "Cannot parse ENDF float: '{}'",
        field
    )))
}

/// Parse an integer from an 11-character ENDF field.
fn parse_endf_int(line: &str, field_index: usize) -> Result<i32, EndfParseError> {
    let start = field_index * 11;
    let end = start + 11;

    if line.len() < end {
        return Ok(0);
    }

    let field = &line[start..end];
    let trimmed = field.trim();

    if trimmed.is_empty() {
        return Ok(0);
    }

    // ENDF integers may have a decimal point (e.g., "1.000000+0" for 1).
    // Try integer parse first, then float-with-integral-value.
    if let Ok(v) = trimmed.parse::<i32>() {
        return Ok(v);
    }

    // Try parsing as float, but reject non-integral values (e.g., "1.9e+0")
    // rather than silently truncating.  Use the same ε=1e-6 tolerance as
    // `parse_tab1`'s NBT/INT validation so that integer fields stored as
    // ENDF floats ("1.000000+0") still round-trip, but a malformed
    // "1.900000+0" is surfaced as an InvalidNumber rather than parsed as 1.
    if let Ok(v) = parse_endf_float(line, field_index) {
        if (v - v.round()).abs() > 1e-6 {
            return Err(EndfParseError::InvalidNumber(format!(
                "Non-integral value in ENDF int field: '{}' (={})",
                field, v
            )));
        }
        return Ok(v.round() as i32);
    }

    Err(EndfParseError::InvalidNumber(format!(
        "Cannot parse ENDF int: '{}'",
        field
    )))
}

/// Consume (skip) a URR LFW=1/LRF=1 "Case B" range body (ENDF-6 §2.2.2.1),
/// advancing the line cursor past the shared fission-width energy LIST and each
/// per-(L,J) LIST record. NEREIDS does not evaluate URR cross sections, so the
/// parameter values are discarded; only the record structure is read, which is
/// what keeps the line stream aligned for any following range or SEND record.
/// Returns the header's `(SPI, AP)` — AP converted to fm — so the caller can
/// preserve them on the placeholder range.
///
/// The per-J control line (carrying `N1 = NE+6`) MUST be consumed before its
/// body — omitting it misaligns the stream by one record per J-group. The
/// `N1 = NE+6` guard is preserved to catch malformed records.
///
/// Reference: ENDF-6 §2.2.2.1 Case B; OpenScale `File2.cpp`
/// (`lfw==1 && lrf==1` branch) and AMPX `File2Unres.f90` — both OpenScale/AMPX
/// readers vendored under SAMMY's external tree, not SAMMY proper.
fn skip_urr_lfw1_lrf1(lines: &[&str], pos: &mut usize) -> Result<(f64, f64), EndfParseError> {
    // CONT: SPI, AP, LSSF, 0, NE, NLS
    let header = parse_cont(lines, pos)?;
    let spi = header.c1;
    let ap_fm = header.c2 * ENDF_RADIUS_TO_FM;
    let ne = checked_count(header.n1, "NE")?;
    let nls = checked_count(header.n2, "NLS")?;

    // LIST: NE shared fission-width-grid energies.
    parse_list_values(lines, pos, ne)?;

    for _ in 0..nls {
        // CONT: AWRI, 0, L, 0, NJS, 0
        let l_cont = parse_cont(lines, pos)?;
        let njs = checked_count(l_cont.n1, "NJS")?;

        // Zero NJS means no J-groups for this L-value, which is malformed
        // (ENDF §2.2.2.1 requires at least one J-group per L-group) and would
        // under-consume the body, surfacing as a confusing downstream error.
        // Mirrors the NJS=0 guards in `skip_urr_range`.
        if njs == 0 {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "URR LFW=1/LRF=1 L={}: NJS=0 (at least one J-group required)",
                l_cont.l1
            )));
        }

        for _ in 0..njs {
            // Per-J LIST control: 0.0, 0.0, L, MUF, NE+6, 0.  Consuming this
            // control record keeps the line stream aligned — the body follows.
            let j_cont = parse_cont(lines, pos)?;
            let n1 = checked_count(j_cont.n1, "N1")?;
            // SCALE validates this exact relation (File2.cpp: N1-6 == NE).
            if n1 != ne + 6 {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "URR LFW=1/LRF=1: per-J N1={n1} ≠ NE+6={} (NE={ne})",
                    ne + 6
                )));
            }
            // LIST body: [D, AJ, AMUN, GNO, GG, 0] + NE fission widths (discarded).
            parse_list_values(lines, pos, n1)?;
        }
    }

    Ok((spi, ap_fm))
}

/// Consume (skip) a URR (LRU=2) range body of the given LRF, advancing the line
/// cursor past every L- and (L,J)-record. NEREIDS does not evaluate URR cross
/// sections, so the parameter values are discarded; only the record structure
/// is read, which keeps the line stream aligned for any following range or SEND
/// record. Returns the header's `(SPI, AP)` — AP converted to fm — so the
/// caller can preserve them on the placeholder range.
///
/// Handles LFW=0/LRF=1 (energy-independent widths) and LRF=2 (tabulated
/// widths); LFW=1/LRF=2 is byte-identical to LFW=0/LRF=2 and routes here too.
/// LFW=1/LRF=1 (shared-grid "Case B") is handled by `skip_urr_lfw1_lrf1`.
///
/// The structural guards (N1=6*NJS for LRF=1; INT∈1..=5 and N1=6*(NE+1) for
/// LRF=2; NJS>0) are preserved: they catch malformed records that would
/// otherwise over-/under-consume lines and misalign the cursor.
///
/// Reference: ENDF-6 Formats Manual §2.2.2
fn skip_urr_range(lines: &[&str], pos: &mut usize, lrf: i32) -> Result<(f64, f64), EndfParseError> {
    // CONT: SPI, AP, 0, 0, NLS, 0
    let spi_cont = parse_cont(lines, pos)?;
    let spi = spi_cont.c1;
    let ap_fm = spi_cont.c2 * ENDF_RADIUS_TO_FM;
    let nls = checked_count(spi_cont.n1, "NLS")?;

    if lrf == 1 {
        // LRF=1: energy-independent widths, one LIST block per L covering all J.
        for _ in 0..nls {
            // CONT: AWRI, 0, L, 0, N1=6*NJS, N2=NJS
            let l_cont = parse_cont(lines, pos)?;
            if l_cont.l1 < 0 {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "URR LRF=1: negative L={}",
                    l_cont.l1
                )));
            }
            let l = l_cont.l1 as u32;
            let n1 = checked_count(l_cont.n1, "N1")?; // 6*NJS
            let njs = checked_count(l_cont.n2, "NJS")?;
            if njs == 0 || n1 != 6 * njs {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "URR LRF=1 L={l}: N1={n1} ≠ 6×NJS={} (NJS={njs})",
                    6 * njs
                )));
            }
            // LIST body: [D, AJ, AMUN, GNO, GG, GF] × NJS (discarded).
            parse_list_values(lines, pos, n1)?;
        }
    } else {
        // LRF=2: energy-dependent width tables, one LIST per (L, J).
        for _ in 0..nls {
            // CONT: AWRI, 0, L, 0, NJS, 0
            let l_cont = parse_cont(lines, pos)?;
            if l_cont.l1 < 0 {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "URR LRF=2: negative L={}",
                    l_cont.l1
                )));
            }
            let l = l_cont.l1 as u32;
            let njs = checked_count(l_cont.n1, "NJS")?; // N1 = NJS for LRF=2

            // Zero NJS means no J-groups for this L-value, which is malformed
            // (ENDF §2.2.2.2 requires at least one J-group per L-group).
            if njs == 0 {
                return Err(EndfParseError::UnsupportedFormat(format!(
                    "URR LRF=2 L={l}: NJS=0 (at least one J-group required)"
                )));
            }

            for _ in 0..njs {
                // CONT: AJ, 0, INT, 0, N1=6*(NE+1), N2=NE
                let j_cont = parse_cont(lines, pos)?;
                let int_code = j_cont.l1; // interpolation law (L1 field)
                // ENDF-6 §0.5 defines INT codes 1..=5; anything outside that
                // range — including negative values and INT=0 or INT≥6 — is a
                // malformed record.
                if !(1..=5).contains(&int_code) {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "URR LRF=2: INT={int_code} out of spec (expected 1..=5)"
                    )));
                }
                let n1 = checked_count(j_cont.n1, "N1")?; // 6*(NE+1)
                let ne = checked_count(j_cont.n2, "NE")?; // NE (number of energy points)

                // Validate N1 = 6*(NE+1) before consuming the LIST body so a
                // malformed record cannot over-/under-consume lines.
                let expected_n1 = 6 * (ne + 1);
                if n1 != expected_n1 {
                    return Err(EndfParseError::UnsupportedFormat(format!(
                        "URR LRF=2: N1={n1} ≠ 6*(NE+1)={expected_n1} (NE={ne})"
                    )));
                }
                // LIST body: DOF row + NE width rows (discarded).
                parse_list_values(lines, pos, n1)?;
            }
        }
    }

    Ok((spi, ap_fm))
}

/// Parse a TAB1 record into a `Tab1` interpolation table.
///
/// ENDF TAB1 layout (Reference: ENDF-6 Formats Manual §0.5):
/// ```text
/// CONT: [C1, C2, L1, L2, NR, NP]
/// NR×2 integer values: (NBT_i, INT_i) pairs  — 6 per line
/// NP×2 float values:   (x_i,   y_i)   pairs  — 6 per line
/// ```
///
/// INT codes: 1=histogram, 2=lin-lin, 3=log-x/lin-y, 4=lin-x/log-y, 5=log-log.
fn parse_tab1(lines: &[&str], pos: &mut usize) -> Result<Tab1, EndfParseError> {
    parse_tab1_with(lines, pos, false)
}

/// A TAB1 record, with `jumps` read as File 3 is.
///
/// With `jumps`, points at one x are a jump from the first value to the last,
/// which holds from x on, and a point repeated exactly is read once. A region
/// left with no interval is dropped.
///
/// # References
/// ENDF-6 Formats Manual (ENDF-102, 2023), §0.5.2.1 and §3.3.1.
fn parse_tab1_with(lines: &[&str], pos: &mut usize, jumps: bool) -> Result<Tab1, EndfParseError> {
    let cont = parse_cont(lines, pos)?;
    let nr = checked_count(cont.n1, "NR")?; // number of interpolation regions
    let np = checked_count(cont.n2, "NP")?; // number of data points

    // NR=0 is valid ENDF: it means a single implicit interpolation region
    // covering all NP points with no explicit boundary record.  The
    // evaluate() call will fall through to the `unwrap_or(2)` default in
    // interp_code_for_interval(), which correctly returns INT=2 (lin-lin).
    // When NR=0, the loop below is a no-op and the interp_raw vec stays empty.

    // Read NR×2 integers: (NBT, INT) pairs packed as ENDF floats.
    // Validate that values are integers, INT codes are in 1..=5, boundaries
    // are strictly increasing, and the last boundary equals NP.
    let interp_raw = parse_list_values(lines, pos, nr * 2)?;
    let mut boundaries = Vec::with_capacity(nr);
    let mut interp_codes = Vec::with_capacity(nr);
    for i in 0..nr {
        let nbt_raw = interp_raw[i * 2];
        let int_raw = interp_raw[i * 2 + 1];

        // ENDF stores integers as floats (e.g. "2.000000+0").  They must be
        // exact whole numbers.  Use a small epsilon (1e-6) rather than the
        // half-unit tolerance 0.5, which would silently accept 1.4 or 2.49.
        // NBT is a 1-based index (ENDF §0.5), so 0 is invalid.
        if (nbt_raw - nbt_raw.round()).abs() > 1e-6 || nbt_raw < 1.0 {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 NBT[{}] is not a positive integer: {}",
                i, nbt_raw
            )));
        }
        if (int_raw - int_raw.round()).abs() > 1e-6 {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 INT[{}] is not an integer: {}",
                i, int_raw
            )));
        }
        let int_code = int_raw.round() as u32;
        if !(1..=5).contains(&int_code) {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 INT[{}]={} is out of range 1..=5",
                i, int_code
            )));
        }
        let nbt = nbt_raw.round() as usize;

        // Boundaries must be strictly increasing (ENDF §0.5).
        if let Some(&prev) = boundaries.last()
            && nbt <= prev
        {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 NBT[{}]={} is not greater than NBT[{}]={}",
                i,
                nbt,
                i - 1,
                prev
            )));
        }
        boundaries.push(nbt);
        interp_codes.push(int_code);
    }

    // The final boundary must equal NP (ENDF §0.5: last NBT is 1-based index of last point).
    if nr > 0 {
        let last_nbt = *boundaries.last().unwrap();
        if last_nbt != np {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 last NBT={} does not equal NP={}",
                last_nbt, np
            )));
        }
    }

    if np == 0 {
        return Err(EndfParseError::UnsupportedFormat(
            "TAB1 NP=0: table must have at least one point".to_string(),
        ));
    }

    // Read NP×2 floats: (E, AP) pairs.
    let data_raw = parse_list_values(lines, pos, np * 2)?;
    let mut points = Vec::with_capacity(np);
    for i in 0..np {
        let x = data_raw[i * 2];
        let y = data_raw[i * 2 + 1];
        if !(x.is_finite() && y.is_finite()) {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 point {i} ({x}, {y}) is not finite"
            )));
        }
        // A jump has two values; one between them holds over no interval.
        if jumps
            && let &[.., (x_first, _), (x_last, y_last)] = points.as_slice()
            && x_first == x
            && x_last == x
            && y_last != y
        {
            points.pop();
            drop_point(&mut boundaries, &mut interp_codes, points.len());
        }
        if jumps && points.last() == Some(&(x, y)) {
            drop_point(&mut boundaries, &mut interp_codes, points.len());
            continue;
        }
        // x-values must be strictly increasing, or repeat once at a jump;
        // Tab1::evaluate() relies on this.
        let jump = jumps && points.last().is_some_and(|&(x_prev, _)| x == x_prev);
        if let Some(&(x_prev, _)) = points.last()
            && x <= x_prev
            && !jump
        {
            return Err(EndfParseError::UnsupportedFormat(format!(
                "TAB1 x[{}]={} is not greater than x[{}]={} (x must be strictly increasing)",
                i,
                x,
                i - 1,
                x_prev
            )));
        }
        points.push((x, y));
    }

    Ok(Tab1 {
        boundaries,
        interp_codes,
        points,
    })
}

/// Remove the 0-based point `index` from a TAB1's interpolation regions, and
/// drop a region it leaves with no interval.
fn drop_point(boundaries: &mut Vec<usize>, interp_codes: &mut Vec<u32>, index: usize) {
    for nbt in boundaries.iter_mut().filter(|nbt| **nbt > index) {
        *nbt -= 1;
    }
    // A region left with no interval: NBT = 1 or NBT equal to the one before.
    let empty = if boundaries.first() == Some(&1) {
        Some(0)
    } else {
        (1..boundaries.len()).find(|&r| boundaries[r] == boundaries[r - 1])
    };
    if let Some(empty) = empty {
        boundaries.remove(empty);
        interp_codes.remove(empty);
    }
}

/// Errors from ENDF parsing.
#[derive(Debug, thiserror::Error)]
pub enum EndfParseError {
    #[error("Missing section: {0}")]
    MissingSection(String),

    #[error("Unsupported format: {0}")]
    UnsupportedFormat(String),

    #[error("Invalid number: {0}")]
    InvalidNumber(String),

    #[error("Unexpected end of file: {0}")]
    UnexpectedEof(String),

    #[error("Invalid isotope: {0}")]
    InvalidIsotope(#[from] nereids_core::error::NereidsError),
}

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

    #[test]
    fn test_parse_endf_float_standard() {
        // ENDF fields are exactly 11 chars wide, no separators.
        // " 1.23456+2" in 11 chars = " 1.23456+02" (Fortran E11.4 style)
        //  01234567890  (field 0: cols 0-10, field 1: cols 11-21, etc.)
        let line = " 1.23456+02 2.34567-01 0.00000+00                                            ";
        assert!((parse_endf_float(line, 0).unwrap() - 123.456).abs() < 0.01);
        assert!((parse_endf_float(line, 1).unwrap() - 0.234567).abs() < 1e-6);
        assert!((parse_endf_float(line, 2).unwrap() - 0.0).abs() < 1e-10);
    }

    #[test]
    fn test_parse_endf_float_with_e() {
        // 11-char fields: "1.23456E+02" "2.34567E-01"
        let line = "1.23456E+022.34567E-01                                                       ";
        assert!((parse_endf_float(line, 0).unwrap() - 123.456).abs() < 0.01);
        assert!((parse_endf_float(line, 1).unwrap() - 0.234567).abs() < 1e-6);
    }

    #[test]
    fn test_parse_endf_float_negative() {
        let line = "-1.23456+02-2.34567-01                                                       ";
        assert!((parse_endf_float(line, 0).unwrap() - (-123.456)).abs() < 0.01);
        assert!((parse_endf_float(line, 1).unwrap() - (-0.234567)).abs() < 1e-6);
    }

    /// Fortran exponents with a space between the sign and digit — e.g. "9.22330+ 4"
    /// — appear in some older ENDF evaluations (observed in SAMMY tr149/t149a.endf
    /// for U-233).  The parser strips the space before parsing the exponent.
    #[test]
    fn test_parse_endf_float_spaced_exponent() {
        // " 9.22330+ 4" occupies 11 chars: space before mantissa, space before digit
        let line =
            " 9.22330+ 4 1.23400- 2                                                         ";
        assert!((parse_endf_float(line, 0).unwrap() - 92_233.0).abs() < 1.0);
        assert!((parse_endf_float(line, 1).unwrap() - 0.01234).abs() < 1e-6);
    }

    #[test]
    fn test_parse_endf_int() {
        let line = "          0          1          2          3          4          5            ";
        assert_eq!(parse_endf_int(line, 0).unwrap(), 0);
        assert_eq!(parse_endf_int(line, 1).unwrap(), 1);
        assert_eq!(parse_endf_int(line, 2).unwrap(), 2);
    }

    /// Parse the vendored U-238 ENDF file (Reich-Moore, LRF=3).
    ///
    /// This test validates against the public-domain U-238 ENDF/B-VIII.0
    /// evaluation shipped at `examples/data/u238_ex027.endf` (the same
    /// file SAMMY distributes as `samexm_new/ex027_new/ex027.endf`). The
    /// first positive-energy resonance of U-238 is at 6.674 eV.
    ///
    /// Vendored under public-domain ENDF/B redistribution, so this gate
    /// runs unconditionally on CI — no `Skipping…` fall-through.
    #[test]
    fn test_parse_u238_sammy_endf() {
        // Crate-local copy so the test works when nereids-endf 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.
        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 = parse_endf_file2(&endf_text).unwrap();

        // Basic structure checks.
        assert_eq!(data.za, 92238, "Should be U-238");
        assert!((data.awr - 236.006).abs() < 0.01, "AWR should be ~236");
        assert!(!data.ranges.is_empty(), "Should have at least one range");

        let range = &data.ranges[0];
        assert!(range.resolved, "First range should be resolved");
        assert_eq!(
            range.formalism,
            ResonanceFormalism::ReichMoore,
            "U-238 ENDF uses Reich-Moore (LRF=3)"
        );
        assert!(
            (range.target_spin - 0.0).abs() < 1e-10,
            "U-238 target spin I=0"
        );
        assert!(
            (range.scattering_radius - 9.4285).abs() < 0.01,
            "Scattering radius ~9.4285 fm (ENDF 0.94285 × 10)"
        );
        assert_eq!(range.l_groups.len(), 2, "Should have L=0 and L=1 groups");

        // Check first L-group (L=0).
        let l0 = &range.l_groups[0];
        assert_eq!(l0.l, 0, "First group should be L=0");
        assert!(
            l0.resonances.len() > 500,
            "L=0 should have hundreds of resonances"
        );

        // Find the famous 6.674 eV resonance of U-238.
        let first_positive = l0
            .resonances
            .iter()
            .find(|r| r.energy > 0.0)
            .expect("Should have positive-energy resonances");
        assert!(
            (first_positive.energy - 6.674).abs() < 0.01,
            "First positive resonance should be at 6.674 eV, got {}",
            first_positive.energy
        );
        assert!(
            (first_positive.j - 0.5).abs() < 1e-10,
            "6.674 eV resonance has J=0.5"
        );

        // The 6.674 eV resonance neutron width: ~1.493e-3 eV
        assert!(
            (first_positive.gn - 1.493e-3).abs() < 1e-5,
            "Neutron width should be ~1.493e-3 eV, got {}",
            first_positive.gn
        );
        // Gamma width: ~2.3e-2 eV
        assert!(
            (first_positive.gg - 2.3e-2).abs() < 1e-3,
            "Gamma width should be ~2.3e-2 eV, got {}",
            first_positive.gg
        );

        // The aggregate accessor must cover at least the resolved range's
        // resonances (it sums over every range, including any unresolved
        // range the file carries).
        let total = data.total_resonance_count();
        let resolved_sum: usize = range.l_groups.iter().map(|g| g.resonances.len()).sum();
        assert!(
            total >= resolved_sum && resolved_sum > 500,
            "total_resonance_count ({total}) must cover the resolved range's \
             {resolved_sum} resonances"
        );
    }

    /// Pin the Ta-181 ENDF/B-VIII.0 resonance count at 76 (genuinely-sparse RRR).
    ///
    /// Ta-181 (MAT 7328) in ENDF/B-VIII.0 has NER=2 ranges:
    ///   range 0: LRU=1 resolved, LRF=2 (MLBW), one L-group, 76 discrete
    ///            resonances, resolved region only to 330 eV;
    ///   range 1: LRU=2 unresolved (URR), 0 discrete resonances.
    /// So `total_resonance_count()` == 76 is **faithful**, not a dropped range —
    /// the parser reads every NER range and errors on unconsumed MF2/MT151 data.
    /// ENDF/B-VIII.1 later extended the resolved region (RRR to 2554 eV, 565
    /// resonances); a low VIII.0 count reflects that evaluation's resolved-region
    /// extent, not a parser bug. This test is a regression guard for that fact.
    ///
    /// Vendored under public-domain ENDF/B redistribution (73-Ta-181 LLNL EVAL,
    /// same evaluation NNDC/IAEA distribute) at the workspace-root
    /// `tests/data/endf/Ta-181.endf` (Hf-177 precedent). Inside the full NEREIDS
    /// workspace — where CI runs — the fixture is always present and the gate
    /// asserts fully; a standalone crate build (no workspace fixtures) skips.
    #[test]
    fn test_parse_ta181_endf8_0_resonance_count() {
        let endf_path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .parent()
            .unwrap()
            .parent()
            .unwrap()
            .join("tests/data/endf/Ta-181.endf");
        // The fixture is committed to the repo, so a missing file is a
        // packaging/path regression that MUST fail CI — not a skip that would
        // silently disable this resonance-count guard.
        let endf_text = std::fs::read_to_string(&endf_path).unwrap_or_else(|e| {
            panic!(
                "vendored Ta-181 regression fixture must be present at \
                 {endf_path:?} (committed test data): {e}"
            )
        });
        let data = parse_endf_file2(&endf_text).unwrap();

        assert_eq!(data.za, 73181, "Should be Ta-181");
        assert!(
            (data.awr - 179.3936).abs() < 0.01,
            "AWR should be ~179.3936, got {}",
            data.awr
        );

        // Faithful count: VIII.0's RRR is genuinely sparse (76), not a drop.
        assert_eq!(
            data.total_resonance_count(),
            76,
            "Ta-181 VIII.0 has 76 discrete resonances (sparse RRR to 330 eV)"
        );
        assert_eq!(
            data.ranges.len(),
            2,
            "Ta-181 VIII.0 has NER=2 ranges (resolved MLBW + unresolved URR)"
        );

        // Range 0: resolved MLBW (LRF=2), all 76 resonances live here.
        let resolved = &data.ranges[0];
        assert!(resolved.resolved, "Range 0 must be resolved (LRU=1)");
        assert_eq!(
            resolved.formalism,
            ResonanceFormalism::MLBW,
            "Ta-181 VIII.0 resolved range uses MLBW (LRF=2)"
        );
        assert_eq!(
            resolved.resonance_count(),
            76,
            "All 76 resonances belong to the resolved MLBW range"
        );

        // Range 1: unresolved (LRU=2, URR), zero discrete resonances.
        let unresolved = &data.ranges[1];
        assert!(!unresolved.resolved, "Range 1 must be unresolved (LRU=2)");
        assert_eq!(
            unresolved.formalism,
            ResonanceFormalism::Unresolved,
            "Ta-181 VIII.0 second range is the URR (LRU=2)"
        );
        assert_eq!(
            unresolved.resonance_count(),
            0,
            "URR range carries no discrete resonances"
        );

        // Mixed evaluation: the resolved MLBW range keeps the file loadable,
        // and the skipped URR range must be flagged as unevaluated so callers
        // can warn that its span contributes zero cross-section.
        assert!(
            data.has_unevaluated_ranges(),
            "Ta-181's URR range must be flagged as unevaluated"
        );
        assert_eq!(
            data.unevaluated_ranges().len(),
            1,
            "exactly the URR range is unevaluated"
        );
    }

    /// Assert that a fixture whose every range is a parse-and-skip placeholder
    /// is rejected with the no-evaluable-ranges error.
    ///
    /// This still pins cursor advancement: a misaligned skip surfaces as a
    /// *different* error (multi-material guard, unexpected EOF, or a
    /// structural-guard message), not this one — so reaching this specific
    /// message proves the skip consumed exactly the range body.
    #[track_caller]
    fn assert_rejected_no_evaluable_ranges(endf: &str) {
        let err = parse_endf_file2(endf).unwrap_err();
        match err {
            EndfParseError::UnsupportedFormat(ref msg) => {
                assert!(
                    msg.contains("No evaluable resonance ranges"),
                    "expected no-evaluable-ranges rejection, got: {msg}"
                );
            }
            other => panic!("expected UnsupportedFormat, got: {other:?}"),
        }
    }

    /// Parse a minimal hand-crafted ENDF snippet with NRO=1 (energy-dependent
    /// scattering radius).
    ///
    /// The fixture encodes:
    /// - LRF=3 (Reich-Moore), NRO=1
    /// - AP TAB1: 2 points — ENDF values 8.0 and 10.0 (10⁻¹² cm),
    ///   which become 80.0 fm and 100.0 fm after ×10 conversion
    /// - One L-group (L=0) with one resonance at 6.674 eV
    ///
    /// Verifies:
    /// - ap_table is Some after parsing
    /// - ap_table.evaluate(1.0) ≈ 80.0 fm (8.0 × ENDF_RADIUS_TO_FM)
    /// - ap_table.evaluate(500.5) ≈ 90.0 fm (midpoint, lin-lin)
    /// - ap_table.evaluate(1000.0) ≈ 100.0 fm
    /// - scattering_radius_at() delegates to the table
    #[test]
    fn test_parse_nro1_tab1() {
        // Each ENDF line is exactly 80 chars: 66 data chars + 14 MAT/MF/MT/SEQ.
        // Cols 67-70: MAT=9237, Cols 71-72: MF=2, Cols 73-75: MT=151, Cols 76-80: seq
        //
        // Line layout (11 chars per field × 6 fields = 66 chars, then 14 control chars):
        //   HEAD:  ZA=92238  AWR=236.006  0  0  NIS=1  0
        //   CONT:  ZAI=92238 ABN=1.0      0  LFW=0 NER=1  0
        //   CONT:  EL=1e-5   EH=1e4    LRU=1  LRF=3  NRO=1  NAPS=0
        //   TAB1 CONT: 0  0  0  0  NR=1  NP=2
        //   TAB1 interp: NBT=2, INT=2  (plus 4 padding zeros)
        //   TAB1 data:   (1.0, 8.0), (1000.0, 10.0)
        //   RM CONT:  SPI=0.0  AP=9.0  0  0  NLS=1  0
        //   L CONT:  AWRI=236.006  0  L=0  0  6*NRS=6  NRS=1
        //   Resonance: ER=6.674  AJ=0.5  GN=1.493e-3  GG=23e-3  GFA=0  GFB=0
        //   SEND: all zeros
        // Each ENDF line: 66 data chars + 4-char MAT(9237) + 2-char MF(" 2")
        //   + 3-char MT("151") + 5-char SEQ = 80 chars total.
        let endf = include_str!("../../../tests/data/synthetic/lrf3_nro1_tab1.endf");

        let data = parse_endf_file2(endf).expect("NRO=1 fixture must parse cleanly");
        assert_eq!(data.ranges.len(), 1, "one energy range");

        let range = &data.ranges[0];
        assert_eq!(
            range.formalism,
            ResonanceFormalism::ReichMoore,
            "must be LRF=3"
        );

        let table = range
            .ap_table
            .as_ref()
            .expect("NRO=1 range must have ap_table");
        assert_eq!(table.points.len(), 2, "TAB1 must have 2 points");

        // Exact boundary values (ENDF 8.0 × 10 = 80.0 fm).
        assert!(
            (table.evaluate(1.0) - 80.0).abs() < 1e-10,
            "AP(1 eV) = 80.0 fm"
        );
        assert!(
            (table.evaluate(1000.0) - 100.0).abs() < 1e-10,
            "AP(1000 eV) = 100.0 fm"
        );
        // Lin-lin midpoint: AP(500.5 eV) ≈ 90.0 fm.
        let mid = table.evaluate(500.5);
        assert!((mid - 90.0).abs() < 0.1, "AP midpoint ≈ 90.0 fm, got {mid}");

        // scattering_radius_at delegates to the table.
        assert!(
            (range.scattering_radius_at(1.0) - 80.0).abs() < 1e-10,
            "scattering_radius_at(1 eV) = 80.0"
        );
        assert!(
            (range.scattering_radius_at(1000.0) - 100.0).abs() < 1e-10,
            "scattering_radius_at(1000 eV) = 100.0"
        );

        // Resonance is still parsed correctly.
        assert_eq!(range.l_groups.len(), 1, "one L-group");
        let res = &range.l_groups[0].resonances[0];
        assert!((res.energy - 6.674).abs() < 1e-6);
    }

    /// LFW=1 with LRF=2 (tabulated widths, U-233-style record).
    ///
    /// SAMMY test tr149 (`t149a.endf`, MAT=9222, ZA=92233) has two ranges:
    ///   - Range 0: LRU=1 (resolved, Reich-Moore / LRF=3)
    ///   - Range 1: LRU=2, LRF=2, **LFW=1** (energy-dependent fission widths)
    ///
    /// ENDF-6 §2.2.2.2: for LFW=1/LRF=2 the per-(L,J) LIST layout is
    /// **identical to LFW=0/LRF=2** (the fission widths are already
    /// per-energy-point in the LIST tail), so the parser dispatches to the
    /// shared `skip_urr_range` path: the URR body is parsed-and-skipped
    /// (structural guards applied, values discarded) and the range is
    /// stored as a non-evaluable placeholder.
    ///
    /// sibling checkout; on CI (and on any clean clone) the file was
    /// absent and the test silently reported `ok` after a `Skipping…`
    /// print. Vendoring the full tr149 ENDF would be heavy; instead we
    /// synthesise a minimal but record-shape-faithful LFW=1/LRF=2
    /// fixture so the assertion runs unconditionally.
    #[test]
    fn test_parse_u233_lfw1_lrf2_urr_parsed() {
        // Minimal MF=2/MT=151 with two ranges, mirroring tr149 layout:
        //   Range 0: LRU=1, LRF=3 (Reich-Moore) — one trivial resonance.
        //   Range 1: LRU=2, LRF=2, LFW=1       — NLS=1, NJS=1, NE=2.
        // LFW=1 is flagged on the isotope CONT (L2 field).
        const ENDF: &str = include_str!("../../../tests/data/synthetic/u233_lfw1_lrf2_urr.endf");

        let data = parse_endf_file2(ENDF)
            .expect("U-233 LFW=1/LRF=2 fixture must parse (record layout = LFW=0/LRF=2)");

        // Both ranges must round-trip: resolved + URR.
        assert_eq!(data.ranges.len(), 2, "must have 2 ranges (resolved + URR)");

        let resolved_count = data.ranges.iter().filter(|r| r.resolved).count();
        assert_eq!(resolved_count, 1, "exactly one resolved range");

        // The URR range is parsed-and-skipped: it must appear as a
        // non-resolved Unresolved range. LFW=1/LRF=2 shares the LFW=0/LRF=2
        // record layout, so a mis-sized skip would trip the multi-material
        // guard or drop the resolved range — the counts above pin alignment.
        let urr_count = data
            .ranges
            .iter()
            .filter(|r| r.formalism == ResonanceFormalism::Unresolved)
            .count();
        assert_eq!(
            urr_count, 1,
            "LFW=1/LRF=2 URR range must be present (skipped)"
        );
        let urr_range = data
            .ranges
            .iter()
            .find(|r| r.formalism == ResonanceFormalism::Unresolved)
            .expect("URR range must exist");
        assert!(!urr_range.resolved, "URR range must not be resolved");

        // Mixed evaluation: the resolved range keeps the file loadable, and
        // the skipped URR range must be flagged as unevaluated.
        assert!(data.has_unevaluated_ranges());
        assert_eq!(data.unevaluated_ranges().len(), 1);
    }

    /// Hand-crafted LRF=1 URR roundtrip test.
    ///
    /// Verifies that a minimal synthetic ENDF snippet with LRU=2, LRF=1 is
    /// parsed correctly: one L-group (L=0), two J-groups with known D, AJ,
    /// AMUN, GNO, GG, GF values.
    #[test]
    fn test_parse_lrf1_urr_roundtrip() {
        // Minimal ENDF MF=2/MT=151 with one resolved range followed by one
        // LRU=2 LRF=1 unresolved range.
        //
        // Resolved range: a simple RM LRF=3 with one resonance (gives the
        // parser something valid to consume before the URR section).
        //
        // URR range: LRU=2, LRF=1, NLS=1 (L=0), NJS=2 J-groups.
        //   J=2.0: D=0.5 eV, AMUN=1, GNO=3e-4 eV, GG=3.5e-2 eV, GF=0
        //   J=3.0: D=0.4 eV, AMUN=1, GNO=2e-4 eV, GG=3.0e-2 eV, GF=1e-3 eV
        //
        // Each ENDF line: 66 data chars + MAT(4) MF(2) MT(3) SEQ(5) = 80 chars.
        const ENDF: &str = include_str!("../../../tests/data/synthetic/lrf1_urr_roundtrip.endf");

        let data = parse_endf_file2(ENDF).expect("LRF=1 URR fixture must parse cleanly");

        // Should have 2 ranges: one resolved + one URR.
        assert_eq!(data.ranges.len(), 2, "must have 2 ranges");

        let urr_range = &data.ranges[1];
        assert!(!urr_range.resolved, "URR range must not be resolved");
        assert_eq!(
            urr_range.formalism,
            ResonanceFormalism::Unresolved,
            "formalism must be Unresolved"
        );

        // The LRF=1 URR body is parsed-and-skipped; the range is present as
        // Unresolved (asserted above). ranges.len()==2 confirms the skip
        // consumed exactly the URR body and left the line stream aligned.

        // Mixed evaluation: the resolved range keeps the file loadable, and
        // the skipped URR range must be flagged as unevaluated.
        assert!(data.has_unevaluated_ranges());
        assert_eq!(data.unevaluated_ranges().len(), 1);
    }

    /// An LRF=2 URR range with INT=3 (log-x/lin-y) passes the per-J INT guard
    /// (§0.5: INT∈1..=5 — a stricter guard would spuriously reject valid
    /// evaluations), then the file is rejected because the URR range is its
    /// only range (no evaluable content).
    #[test]
    fn test_parse_lrf2_urr_int3_rejected_no_evaluable_ranges() {
        // Minimal ENDF MF=2/MT=151 with one LRU=2/LRF=2 range:
        // NLS=1 (L=0), NJS=1, NE=2 energy points, INT=3 (log-x/lin-y).
        // LIST layout: row 0 = [0, 0, 0, AMUN, 0, AMUF]; rows 1..=NE = (E,
        // D, GX, GN, GG, GF). Total = 6*(NE+1) = 18 floats = 3 lines.
        // Reaching the no-evaluable-ranges error (not an INT rejection or a
        // misalignment error) proves the INT=3 guard accepted the record and
        // the skip consumed exactly the URR body.
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lrf2_urr_int3_roundtrip.endf");
        assert_rejected_no_evaluable_ranges(ENDF);
    }

    /// LRF=2 URR with INT=0 is rejected as a hard error.
    ///
    /// ENDF-6 §0.5 defines INT codes 1..=5 only. INT=0 is malformed,
    /// not merely unsupported, so the parser surfaces it as
    /// `UnsupportedFormat("INT=0 out of spec (expected 1..=5)")` rather
    /// than panicking or silently defaulting to lin-lin.
    #[test]
    fn test_parse_lrf2_urr_int0_rejected() {
        // Same skeleton as INT=3 test, but with INT=0 in the J CONT.
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lrf2_urr_int0_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        // Match on the error variant structurally so the test only
        // depends on the parser surfacing `UnsupportedFormat`, not on
        // the exact `Display` wording of `EndfParseError`.  Then check
        // the payload string for the diagnostic substrings.
        match err {
            EndfParseError::UnsupportedFormat(ref msg) => {
                assert!(
                    msg.contains("INT=0") && msg.contains("out of spec"),
                    "expected INT-out-of-spec rejection in UnsupportedFormat payload, \
                     got: {msg}"
                );
            }
            other => panic!("expected EndfParseError::UnsupportedFormat for INT=0, got: {other:?}"),
        }
    }

    // -----------------------------------------------------------------------
    // robustness tests for malformed input validation
    // -----------------------------------------------------------------------

    /// `checked_count` rejects negative values.
    #[test]
    fn test_checked_count_negative() {
        let err = checked_count(-1, "NLS").unwrap_err();
        assert!(
            err.to_string().contains("Negative"),
            "expected negative error, got: {err}"
        );
    }

    /// `checked_count` rejects values above `MAX_ENDF_COUNT` to prevent
    /// allocation bombs from malformed files.
    #[test]
    fn test_checked_count_upper_bound() {
        // Just above the limit.
        let err = checked_count(MAX_ENDF_COUNT + 1, "NRS").unwrap_err();
        assert!(
            err.to_string().contains("too large"),
            "expected upper-bound error, got: {err}"
        );

        // At the limit: should succeed.
        assert_eq!(checked_count(MAX_ENDF_COUNT, "NRS").unwrap(), 1_000_000);

        // i32::MAX: should be rejected.
        let err = checked_count(i32::MAX, "NPL").unwrap_err();
        assert!(
            err.to_string().contains("too large"),
            "expected upper-bound error for i32::MAX, got: {err}"
        );
    }

    /// `parse_endf_int` rejects non-integral float values rather than
    /// silently truncating.
    ///
    /// Without the strict check, an INT-field value stored as
    /// `"1.900000+0"` would be cast as `1.9_f64 as i32 == 1`, masking
    /// a malformed evaluation.  After the strict check, the parser
    /// surfaces `InvalidNumber` immediately.
    ///
    /// The two field layouts:
    ///   • field 0: integral float "1.000000+0" → returns Ok(1)
    ///   • field 1: non-integral  "1.900000+0" → returns Err(InvalidNumber)
    /// must both round-trip the standard ENDF integer-as-float encoding.
    #[test]
    fn test_parse_endf_int_rejects_non_integral_float() {
        // 11-char fields, padded with leading space to match the ENDF column
        // width.  Field 0 is integral; field 1 is non-integral.
        const LINE: &str = " 1.000000+0 1.900000+0";

        // Integral float must parse cleanly.
        let ok = parse_endf_int(LINE, 0).expect("integral float must parse");
        assert_eq!(ok, 1);

        // Non-integral float must be rejected, not truncated to 1.
        let err = parse_endf_int(LINE, 1).expect_err("non-integral must be rejected");
        match err {
            EndfParseError::InvalidNumber(ref msg) => {
                assert!(
                    msg.contains("Non-integral"),
                    "expected Non-integral diagnostic, got: {msg}"
                );
            }
            other => panic!("expected InvalidNumber, got: {other:?}"),
        }
    }

    /// Negative L-value in a Breit-Wigner range is rejected.
    ///
    /// Constructs a minimal SLBW fixture with L=-1 in the L-group CONT record.
    /// Without the validation, `l_cont.l1 as u32` would wrap to `u32::MAX`.
    #[test]
    fn test_bw_negative_l_rejected() {
        // Minimal SLBW fixture: HEAD + isotope CONT + range CONT + SPI/AP CONT +
        // L-group CONT with L=-1 (field 3 = -1).
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lrf1_bw_negative_l_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("negative L"),
            "expected negative L error, got: {err}"
        );
    }

    /// Negative L-value in a Reich-Moore range is rejected.
    #[test]
    fn test_rm_negative_l_rejected() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lrf3_rm_negative_l_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("negative L"),
            "expected negative L error, got: {err}"
        );
    }

    /// LRU=0 range with non-zero NLS is rejected.
    ///
    /// ENDF-6 §2.2 says the SPI/AP CONT after an LRU=0 range must have
    /// NLS=0 (no L-groups for scattering-radius-only ranges).
    #[test]
    fn test_lru0_nonzero_nls_rejected() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lru0_nonzero_nls_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("NLS=3"),
            "expected LRU=0 NLS validation error, got: {err}"
        );
    }

    /// A resolved (LRU=1) MLBW range with NLS=0 is rejected.
    ///
    /// ENDF-6 §2.2.1.1 requires NLS >= 1 in a resolved range (NLS=0 is
    /// reserved for LRU=0 scattering-radius-only stanzas). Without the guard
    /// this file loads as "evaluable" with zero resonances — zero
    /// cross-section everywhere and no warning.
    #[test]
    fn test_mlbw_nls0_rejected() {
        const ENDF: &str = include_str!("../../../tests/data/synthetic/mlbw_nls0_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("NLS=0"),
            "expected resolved NLS=0 rejection, got: {err}"
        );
    }

    /// A file whose ONLY range is resolved with NLS=1 but NRS=0 in every
    /// L-group is rejected by the no-evaluable-content guard, naming the
    /// inert span. (The range itself is accepted as a warn-and-skip
    /// placeholder — SAMMY/OpenScale accept NRS=0 L-groups — so the
    /// rejection comes from the file-level guard, not a per-range error.)
    #[test]
    fn test_mlbw_nrs0_everywhere_rejected() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/mlbw_nrs0_everywhere_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("No evaluable resonance ranges"),
            "expected no-evaluable-content rejection, got: {msg}"
        );
        assert!(
            msg.contains("resolved range with no resonances"),
            "rejection must name the inert resolved span, got: {msg}"
        );
    }

    /// A mixed file [populated resolved range][all-empty resolved range]
    /// loads: the populated range stays evaluable, the inert range becomes a
    /// warn-and-skip placeholder — symmetric with LRF=7/LRU=2/LRU=0 handling,
    /// so the good range is not lost.
    #[test]
    fn test_mlbw_nrs0_plus_resolved_mixed_loads_and_flags() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/mlbw_nrs0_plus_resolved_mixed.endf");

        let data = parse_endf_file2(ENDF).expect("mixed populated + inert file must load");
        assert_eq!(data.ranges.len(), 2);
        assert!(data.has_evaluable_range());
        assert!(data.has_unevaluated_ranges());
        assert_eq!(data.total_resonance_count(), 1);
        let inert = data
            .ranges
            .iter()
            .find(|r| !r.is_evaluable())
            .expect("the inert range must be present");
        assert!(
            inert
                .skip_description()
                .contains("LRF=2 (MLBW) resolved range with no resonances"),
            "skip_description must name the inert shape, got: {}",
            inert.skip_description()
        );
    }

    /// An MLBW range with a competitive width (LRX=1) loads but is not
    /// evaluated; the range before it still is.
    #[test]
    fn test_mlbw_lrx1_range_is_parsed_but_not_evaluated() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/mlbw_lrx1_plus_resolved_mixed.endf");

        let data = parse_endf_file2(ENDF).expect("the LRX=0 range keeps the file loadable");
        assert_eq!(data.ranges.len(), 2);
        let competitive = &data.ranges[1].l_groups[0];
        assert_eq!((competitive.lrx, competitive.qx), (1, -4.49e4));
        assert!(data.ranges[0].is_evaluable());
        assert_eq!(data.unevaluated_ranges().len(), 1);
        assert!(
            data.ranges[1]
                .skip_description()
                .contains("LRF=2 (MLBW) resolved range with a competitive width (LRX=1)"),
            "got: {}",
            data.ranges[1].skip_description()
        );
    }

    /// A resolved (LRU=1) Reich-Moore range with NLS=0 is rejected.
    #[test]
    fn test_rm_nls0_rejected() {
        const ENDF: &str = include_str!("../../../tests/data/synthetic/lrf3_rm_nls0_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("NLS=0"),
            "expected resolved NLS=0 rejection, got: {err}"
        );
    }

    /// LFW=1/LRF=1 (energy-dependent fission widths) URR is skipped with the
    /// line stream left aligned.
    ///
    /// ENDF-6 §2.2.2.1 Case B: a shared NE-point energy grid is followed,
    /// for each (L, J), by a full LIST record — a control line
    /// `[0.0, 0.0, L, MUF, NE+6, 0]` and then a body
    /// `[D, AJ, AMUN, GNO, GG, 0] + GF(1..NE)`.  The per-J control line MUST
    /// be consumed before the body; otherwise the line stream misaligns by
    /// one record per J-group and the file fails to parse (the stray lines
    /// trip the multi-material guard).
    ///
    /// This fixture is standards-compliant (it includes the per-J control
    /// line); the assertion that it parses end-to-end is the regression guard
    /// for the Case-B skip's cursor advancement.
    ///
    /// The fixture has NE=2, NLS=1, NJS=1, MUF=1.
    #[test]
    fn test_lfw1_lrf1_urr_rejected_no_evaluable_ranges() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lfw1_urr_gracefully_skipped.endf");

        // The Case-B URR body is parsed-and-skipped; consuming the per-J
        // control line (N1=NE+6) keeps the stream aligned. The fixture's only
        // range is the URR, so the aligned parse then hits the
        // no-evaluable-ranges rejection — a misaligned skip would instead trip
        // the multi-material guard or EOF (see the fixture note above), so
        // reaching this specific error pins the cursor advancement.
        assert_rejected_no_evaluable_ranges(ENDF);
    }

    /// A per-J LIST control whose `N1 != NE+6` is a malformed Case-B record.
    /// The parser must reject it — this covers the per-J N1 validation guard
    /// (the SCALE `list.getN1()-6 == ener.getNtot()` relation). The fixture is
    /// byte-identical to the valid one except the per-J control N1 (8 -> 7).
    #[test]
    fn test_lfw1_lrf1_urr_rejects_bad_perj_n1() {
        const ENDF: &str = include_str!("../../../tests/data/synthetic/lfw1_urr_bad_perj_n1.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("N1=") && err.to_string().contains("NE+6"),
            "expected per-J N1 != NE+6 rejection, got: {err}"
        );
    }

    /// An LRU=2 (URR) range with LRF ∉ {1, 2} is a malformed record and is
    /// rejected with an error naming the offending LRF.
    ///
    /// ENDF-6 §2.2.2 restricts the unresolved region to LRF=1 or LRF=2. The
    /// span silently; it now hard-errors like the sibling malformed-record
    /// guards.
    #[test]
    fn test_urr_lrf3_rejected() {
        const ENDF: &str = include_str!("../../../tests/data/synthetic/urr_lrf3_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("LRF=3") && msg.contains("URR"),
            "expected URR LRF=3 rejection naming the LRF, got: {msg}"
        );
    }

    /// LRU=0 range with non-zero L1 in SPI/AP CONT is rejected.
    ///
    /// ENDF-6 §2.2: the SPI/AP CONT record for LRU=0 must be
    /// [SPI, AP, 0, 0, NLS=0, 0].  Non-zero L1 or L2 indicates a
    /// malformed or mis-identified record.
    #[test]
    fn test_lru0_nonzero_l1_rejected() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lru0_nonzero_l1_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("L1=5"),
            "expected LRU=0 L1 validation error, got: {err}"
        );
    }

    /// A well-formed LRU=0-only evaluation (scattering radius, no resonance
    /// parameters) is rejected — NEREIDS cannot evaluate it and loading it
    /// would yield silent zero cross-sections — but the message NAMES the
    /// LRU=0 span rather than misreporting an empty file.
    ///
    /// ENDF-6 §2.1: LRU=0 gives a scattering radius but no resonances. The
    /// fixture is byte-identical to `lru0_nonzero_nls_rejected` except the
    /// SPI/AP CONT carries NLS=0 (well-formed), so parsing reaches the
    /// no-evaluable-content guard rather than the NLS validation guard.
    #[test]
    fn test_lru0_only_rejected_names_lru0_span() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lru0_only_rejected_no_evaluable.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("No evaluable resonance ranges"),
            "expected no-evaluable-ranges rejection, got: {msg}"
        );
        assert!(
            msg.contains("LRU=0 (scattering-radius-only, no resonance parameters)"),
            "rejection must name the LRU=0 span, got: {msg}"
        );
    }

    /// A mixed evaluation carrying an LRU=0 range plus a resolved range still
    /// loads: the resolved range keeps the file evaluable, and the LRU=0 span
    /// is surfaced as a non-evaluable parse-and-skip placeholder (same handling
    /// as LRF=7 / LRU=2).
    #[test]
    fn test_lru0_plus_resolved_mixed_loads_and_flags_lru0() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lru0_plus_resolved_mixed.endf");

        let data = parse_endf_file2(ENDF).expect("mixed LRU=0 + resolved file must load");
        assert_eq!(data.ranges.len(), 2, "LRU=0 range + resolved range");
        assert!(
            data.has_evaluable_range(),
            "resolved range keeps it evaluable"
        );
        assert!(data.has_unevaluated_ranges(), "LRU=0 range is unevaluated");

        let lru0 = data
            .ranges
            .iter()
            .find(|r| r.formalism == ResonanceFormalism::ScatteringRadiusOnly)
            .expect("the LRU=0 range must be present");
        assert!(!lru0.resolved, "LRU=0 range must not be resolved");
        assert!(
            lru0.skip_description()
                .contains("LRU=0 (scattering-radius-only, no resonance parameters)"),
            "LRU=0 skip_description must name the span, got: {}",
            lru0.skip_description()
        );

        // The fixture gives the placeholder DIFFERENT SPI/AP (1.5 / 5.0 fm)
        // than the resolved range (2.5 / 9.6931 fm), so consumers that must
        // prefer the evaluable range are distinguishable from ones reading
        // the placeholder.
        assert_eq!(lru0.target_spin, 1.5);
        assert_eq!(lru0.scattering_radius, 5.0);
        let resolved = data
            .ranges
            .iter()
            .find(|r| r.is_evaluable())
            .expect("the resolved range must be evaluable");
        assert_eq!(resolved.target_spin, 2.5);
        assert!((resolved.scattering_radius - 9.6931).abs() < 1e-9);
    }

    /// N1 != 6*NRS in a BW range CONT is rejected.
    #[test]
    fn test_bw_n1_nrs_mismatch_rejected() {
        const ENDF: &str =
            include_str!("../../../tests/data/synthetic/lrf1_bw_n1_nrs_mismatch_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("N1=7"),
            "expected N1/NRS mismatch error, got: {err}"
        );
    }

    /// Multi-MAT detection: unconsumed MF=2/MT=151 lines after the first
    /// material are rejected.
    #[test]
    fn test_multi_mat_detection() {
        // A valid single-range SLBW file with an extra trailing data line
        // that still carries MF=2/MT=151 tags.
        const ENDF: &str = include_str!("../../../tests/data/synthetic/multi_mat_detection.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        assert!(
            err.to_string().contains("Multiple materials"),
            "expected multi-MAT error, got: {err}"
        );
    }

    /// MF=2 NIS>1 multi-isotope materials are rejected.
    ///
    /// ENDF-6 §2.1 allows a single material to carry several isotopes, each
    /// with its own ZAI/ABN/NER subsection. NEREIDS's `ResonanceData` cannot
    /// represent that hierarchy without losing per-isotope abundance weights;
    /// rather than silently flatten the ranges into one isotope, the parser
    /// returns `UnsupportedFormat` so a downstream consumer cannot be tricked
    /// into computing an abundance-blind cross section.
    ///
    /// The synthetic fixture is a minimal two-isotope SLBW material: Cu-63
    /// (ZAI=29063, ABN=0.6917) and Cu-65 (ZAI=29065, ABN=0.3083), each with
    /// a single L=0 resonance. The parser should reject the file as soon as
    /// it reads the HEAD record (NIS=2 > 1) and never advance into either
    /// isotope subsection.
    #[test]
    fn test_parse_endf_rejects_nis_gt_1() {
        // Minimal NIS=2 fixture. Both isotope subsections use LRF=1 SLBW so
        // that, if the NIS guard were ever removed, the parser would have a
        // valid stream to walk and the assertion below would still pin the
        // guard. The HEAD's ZA is set to the natural-element identifier
        // ZA=29000 (nat-Cu); for NIS=1 callers that would be rejected by
        // `isotope_from_za` with an "A=0" error, but the NIS=2 check runs
        // first and returns the expected UnsupportedFormat.
        const ENDF: &str = include_str!("../../../tests/data/synthetic/nis_gt_1_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        match &err {
            EndfParseError::UnsupportedFormat(msg) => {
                assert!(
                    msg.contains("NIS=2"),
                    "expected NIS=2 in error message, got: {msg}"
                );
                assert!(
                    msg.contains("multi-isotope"),
                    "expected 'multi-isotope' in error message, got: {msg}"
                );
            }
            other => panic!("expected UnsupportedFormat for NIS>1, got {other:?}"),
        }
    }

    /// MF=2 NIS=0 (no isotopes declared) is rejected up-front.
    ///
    /// ENDF-6 §2.1 requires NIS >= 1 for a valid resonance evaluation.
    /// Without the explicit reject, NIS=0 would fall through the per-isotope
    /// loop (zero iterations), leave the resonance section empty, and trip a
    /// confusing downstream "unconsumed data lines" / empty-range failure
    /// far from the actual root cause. The reject mirrors the NIS>1 guard
    /// pattern so both invalid extremes return a clear UnsupportedFormat.
    #[test]
    fn test_parse_endf_rejects_nis_zero() {
        // Minimal NIS=0 fixture: just the HEAD line with NIS=0. The HEAD's
        // ZA=74184 (W-184) is a valid identifier, so any error must come
        // from the NIS=0 guard, not from `isotope_from_za`.
        // HEAD: ZA=74184, AWR=182, NIS=0
        const ENDF: &str = include_str!("../../../tests/data/synthetic/nis_zero_rejected.endf");

        let err = parse_endf_file2(ENDF).unwrap_err();
        match &err {
            EndfParseError::UnsupportedFormat(msg) => {
                assert!(
                    msg.contains("NIS=0"),
                    "expected NIS=0 in error message, got: {msg}"
                );
                assert!(
                    msg.contains("NIS >= 1"),
                    "expected 'NIS >= 1' guidance in error message, got: {msg}"
                );
            }
            other => panic!("expected UnsupportedFormat for NIS=0, got {other:?}"),
        }
    }

    fn file32_fixture() -> String {
        let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("../../tests/data/endf/Ta-181-file32.endf");
        std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{path:?}: {e}"))
    }

    #[test]
    fn test_file32_covariance_of_the_resonances_in_a_window() {
        let text = file32_fixture();
        let data = parse_endf_file2(&text).unwrap();
        let read = |window| parse_resonance_covariance(&text, &data, window).unwrap();
        let all = read(0.0..=100.0).expect("listed");
        assert_eq!(all.resonances, [[0, 0, 0], [0, 0, 1], [0, 0, 2], [0, 0, 3]]);
        assert_eq!(all.covariance[0], 1.229472e-6);
        let window = read(14.0..=30.0).expect("listed");
        assert_eq!(window.resonances, [[0, 0, 1], [0, 0, 2], [0, 0, 3]]);
        let at = |x: usize, y: usize| window.covariance[9 * x + y];
        for (x, y, expected) in [
            (0, 0, 1.558071e-6),
            (1, 2, -1.216048e-9),
            (0, 3, 6.691704e-9),
            (2, 7, 4.216599e-9),
            (5, 6, 8.726475e-9),
            (8, 8, 1.540406e-6),
        ] {
            assert_eq!((at(x, y), at(y, x)), (expected, expected), "{x}, {y}");
        }
        assert!(read(100.0..=200.0).is_none());
    }

    #[test]
    fn test_file32_forms_not_read_are_refused() {
        let text = file32_fixture();
        let data = parse_endf_file2(&text).unwrap();
        let read = |text: &str| parse_resonance_covariance(text, &data, 0.0..=100.0);
        let edit = |text: &str, line: usize, field: usize, value: &str| {
            let change = |(k, l): (usize, &str)| match k + 1 == line {
                true => format!("{}{value:>11}{}\n", &l[..11 * field], &l[11 * field + 11..]),
                false => format!("{l}\n"),
            };
            text.lines().enumerate().map(change).collect::<String>()
        };
        let edited = |line: usize, field: usize, value: &str| edit(&text, line, field, value);
        let edits = |changes: &[(usize, usize, &str)]| {
            changes
                .iter()
                .fold(text.clone(), |text, &(l, f, v)| edit(&text, l, f, v))
        };
        let unresolved =
            " 2.554000+3 1.000000+5          2          1          0          0732832151
 3.500000+0 7.860000-1          0          0          1          0732832151
 1.793936+2 0.000000+0          0          0          6          1732832151
 1.000000+0 3.000000+0 1.000000-3 6.000000-2 0.000000+0 0.000000+0732832151
 0.000000+0 0.000000+0          1          0          1          1732832151
 1.000000-2                                                       732832151
";
        let head = |text: &str, n| text.split_inclusive('\n').take(n).collect::<String>();
        let only = head(&text, 11) + unresolved;
        for refused in [
            edit(&only, 11, 4, "2"),
            edit(&only, 12, 3, "2"),
            edits(&[
                (17, 0, "1.034682+1"),
                (17, 2, "4.816293-3"),
                (17, 3, "6.164143-2"),
            ]),
            edited(10, 0, "7.318200+4"),
            edited(11, 0, "7.318200+4"),
            edited(10, 4, "2"),
            edits(&[(10, 4, "2"), (12, 2, "2"), (12, 3, "1")]),
            edited(12, 2, "2"),
            edits(&[(11, 4, "2"), (12, 2, "2"), (12, 3, "1")]),
            edited(11, 4, "0"),
            edited(12, 2, "0"),
            edited(12, 3, "2"),
            edited(12, 4, "1"),
            edited(13, 3, "2"),
            edited(13, 5, "1"),
            edited(14, 5, "1"),
            edited(15, 2, "2"),
            edited(15, 4, "101"),
            edited(16, 0, "nan"),
            edited(17, 0, "2.031200+1"),
            edited(17, 1, "4.000000+0"),
            edited(17, 2, "1.066000-3"),
            edited(11, 4, "2") + unresolved + &unresolved[..76],
        ] {
            let refused = read(&refused);
            assert!(
                matches!(refused, Err(EndfParseError::UnsupportedFormat(_))),
                "{refused:?}"
            );
        }
        let all = read(&text).unwrap().expect("listed");
        for accepted in [
            edited(17, 0, "2.031171+1"),
            edited(17, 2, "1.065878-3"),
            edited(11, 4, "2") + unresolved,
        ] {
            assert_eq!(read(&accepted).unwrap().as_ref(), Some(&all));
        }
        assert!(read(&(edited(11, 4, "2") + &unresolved[..380])).is_err());
        let edge = parse_resonance_covariance(&edited(17, 0, "2.031171+1"), &data, 0.0..=20.3117);
        assert_eq!(
            edge.unwrap().expect("listed").resonances,
            all.resonances[..2]
        );
        let mut grouped: Vec<String> = edits(&[(4, 4, "2"), (5, 4, "12"), (5, 5, "2")])
            .lines()
            .map(String::from)
            .collect();
        let group =
            " 1.793936+2 7.860000-1          1          0         12          27328 2151    0";
        grouped.insert(7, group.to_string());
        let grouped = grouped.join("\n");
        let split = parse_endf_file2(&grouped).unwrap();
        let regrouped = parse_resonance_covariance(&grouped, &split, 0.0..=100.0).unwrap();
        let regrouped = regrouped.expect("listed");
        let indices = [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1]];
        assert_eq!(regrouped.resonances, indices);
        assert_eq!(regrouped.covariance, all.covariance);
        let doublet = edits(&[
            (8, 0, "2.031170+1"),
            (8, 2, "1.065877-3"),
            (8, 3, "6.370451-2"),
        ]);
        let mut multilevel = data.clone();
        multilevel.ranges[0].formalism = ResonanceFormalism::MLBW;
        for (file, file2, window) in [
            (&text, multilevel, 0.0..=100.0),
            (&doublet, parse_endf_file2(&doublet).unwrap(), 15.0..=21.0),
        ] {
            let refused = parse_resonance_covariance(file, &file2, window);
            assert!(matches!(refused, Err(EndfParseError::UnsupportedFormat(_))));
        }
        for window in [f64::NAN..=1.0, 30.0..=14.0] {
            let refused = parse_resonance_covariance(&text, &data, window);
            assert!(matches!(refused, Err(EndfParseError::InvalidNumber(_))));
        }
        for empty in [head(&text, 9), only] {
            assert!(read(&empty).unwrap().is_none());
        }
        let field = |v: f64| format!("{:>11}", format!("{v:e}").replace('e', ""));
        let members = [[0, 1], [1, 2]];
        let block = |rows: [usize; 2]| {
            let at = |p: usize| 3 * rows[p / 3] + p % 3;
            let packed = (0..6).flat_map(|p| (p..6).map(move |q| (p, q)));
            let values: Vec<String> = packed
                .map(|(p, q)| field(all.covariance[12 * at(p) + at(q)]))
                .collect();
            let head = "        0.0        0.0          3          0         33          2".into();
            let table = rows.map(|r| text.lines().nth(15 + r).unwrap()[..66].to_string());
            let body = values.chunks(6).map(|c| c.concat()).collect::<Vec<_>>();
            std::iter::once(head)
                .chain(table)
                .chain(body)
                .map(|l| format!("{l:<66}732832151    0\n"))
                .collect::<String>()
        };
        let sections = head(&edited(14, 4, "2"), 14);
        let overlapping = sections + &members.map(block).concat();
        let summed = read(&overlapping).unwrap().expect("listed");
        assert_eq!(summed.resonances, all.resonances[..3]);
        for (a, b) in (0..9).flat_map(|a| (0..9).map(move |b| (a, b))) {
            let shared = members
                .iter()
                .filter(|m| m.contains(&(a / 3)) && m.contains(&(b / 3)));
            let expected = shared.count() as f64 * all.covariance[12 * a + b];
            assert_eq!(summed.covariance[9 * a + b], expected, "{a}, {b}");
        }
    }
}

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

    const PHOTON: [f64; 12] = [
        0.0, 183.37, 0.0, 74.0, 1.0, 0.0, 0.0, 0.0, 0.0, 102.0, 0.0, 0.0,
    ];
    const NEUTRON: [f64; 12] = [
        1.0, 182.37, 0.0, 74.0, 0.5, 0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 1.0,
    ];

    /// One KBK record: LCH, LBK and the LBK=2 list [ED, EU, R0, R1, R2, S0, S1].
    struct Background {
        lch: i64,
        lbk: i64,
        list: Option<[f64; 7]>,
    }

    struct Group {
        aj: f64,
        pj: f64,
        kps: i64,
        channels: Vec<[f64; 6]>,
        resonances: Vec<Vec<f64>>,
        nx: Option<usize>,
        background: Vec<Background>,
    }

    /// One-range LRF=7 MF=2 section written in the ENDF-102 §2.2.1.6 layout.
    struct Fixture {
        krm: i64,
        ifg: i64,
        krl: i64,
        nro: bool,
        pairs: Vec<[f64; 12]>,
        groups: Vec<Group>,
    }

    impl Fixture {
        /// W-184-like evaluation: photon pair 1, neutron pair 2, one s-wave group.
        fn basic() -> Self {
            Self {
                krm: 3,
                ifg: 0,
                krl: 0,
                nro: false,
                pairs: vec![PHOTON, NEUTRON],
                groups: vec![Group {
                    aj: 0.5,
                    pj: 0.0,
                    kps: 0,
                    channels: vec![
                        [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                        [2.0, 0.0, 0.5, 0.0, 0.764, 0.75],
                    ],
                    resonances: vec![vec![-12.0, 0.05, 0.3], vec![101.9, 0.06, -0.02]],
                    nx: None,
                    background: vec![],
                }],
            }
        }

        fn text(&self) -> String {
            fn float(x: f64) -> String {
                format!("{:>11}", format!("{x:.5E}"))
            }
            fn int(n: i64) -> String {
                format!("{n:>11}")
            }
            fn list(record: &mut dyn FnMut(Vec<String>), values: &[f64]) {
                for chunk in values.chunks(6) {
                    let mut fields: Vec<String> = chunk.iter().map(|&v| float(v)).collect();
                    fields.resize(6, " ".repeat(11));
                    record(fields);
                }
            }
            let mut lines: Vec<String> = Vec::new();
            let mut record = |fields: Vec<String>| {
                let seq = lines.len() + 1;
                lines.push(format!("{:66}7437 2151{seq:>5}", fields.concat()));
            };
            let zero = || float(0.0);
            record(vec![
                float(74184.0),
                float(182.37),
                int(0),
                int(0),
                int(1),
                int(0),
            ]);
            record(vec![
                float(74184.0),
                float(1.0),
                int(0),
                int(0),
                int(1),
                int(0),
            ]);
            record(vec![
                float(1e-5),
                float(1e4),
                int(1),
                int(7),
                int(i64::from(self.nro)),
                int(0),
            ]);
            if self.nro {
                record(vec![zero(), zero(), int(0), int(0), int(1), int(2)]);
                record(vec![
                    int(2),
                    int(2),
                    " ".repeat(11),
                    " ".repeat(11),
                    " ".repeat(11),
                    " ".repeat(11),
                ]);
                list(&mut record, &[1e-5, 0.7, 1e4, 0.7]);
            }
            record(vec![
                zero(),
                zero(),
                int(self.ifg),
                int(self.krm),
                int(self.groups.len() as i64),
                int(self.krl),
            ]);
            let npp = self.pairs.len() as i64;
            record(vec![
                zero(),
                zero(),
                int(npp),
                int(0),
                int(12 * npp),
                int(2 * npp),
            ]);
            list(&mut record, &self.pairs.concat());
            for g in &self.groups {
                let nch = g.channels.len();
                record(vec![
                    float(g.aj),
                    float(g.pj),
                    int(g.background.len() as i64),
                    int(g.kps),
                    int(6 * nch as i64),
                    int(nch as i64),
                ]);
                list(&mut record, &g.channels.concat());
                let per = 6 * (nch / 6 + 1);
                let mut values: Vec<f64> = Vec::new();
                for r in &g.resonances {
                    let mut row = r.clone();
                    row.resize(per, 0.0);
                    values.extend(row);
                }
                if g.resonances.is_empty() {
                    values.resize(6, 0.0);
                }
                let nx = g.nx.unwrap_or(values.len() / 6);
                values.resize(6 * nx, 0.0);
                record(vec![
                    zero(),
                    zero(),
                    int(0),
                    int(g.resonances.len() as i64),
                    int(6 * nx as i64),
                    int(nx as i64),
                ]);
                list(&mut record, &values);
                for b in &g.background {
                    record(vec![zero(), zero(), int(b.lch), int(b.lbk), int(0), int(0)]);
                    if let Some(v) = b.list {
                        record(vec![
                            float(v[0]),
                            float(v[1]),
                            int(0),
                            int(0),
                            int(5),
                            int(0),
                        ]);
                        list(&mut record, &v[2..]);
                    }
                }
            }
            lines.join("\n")
        }

        fn parse(&self) -> Result<ResonanceData, EndfParseError> {
            parse_file2_records(&self.text())
        }
    }

    fn rml(data: &ResonanceData) -> &RmlRange {
        data.ranges[0]
            .r_matrix_limited
            .as_ref()
            .expect("parsed LRF=7 payload")
    }

    fn assert_rejects(fixture: &Fixture, needle: &str) {
        match fixture.parse() {
            Err(EndfParseError::UnsupportedFormat(msg)) => {
                assert!(msg.contains(needle), "expected '{needle}' in: {msg}");
            }
            other => panic!("expected UnsupportedFormat with '{needle}', got {other:?}"),
        }
    }

    #[test]
    fn parses_pairs_channels_and_widths() {
        let data = Fixture::basic().parse().unwrap();
        let range = &data.ranges[0];
        assert_eq!(range.formalism, ResonanceFormalism::RMatrixLimited);
        let rml = rml(&data);
        assert_eq!((rml.photon_pair, rml.incident_pair), (0, 1));
        assert!(
            !rml.pairs[0].penetrability,
            "PNT=0 on MT=102 means no penetrability"
        );
        assert!(rml.pairs[1].penetrability);
        let group = &rml.spin_groups[0];
        assert_eq!((group.j, group.parity), (0.5, 1));
        assert_eq!(group.channels[1].pair, 1);
        assert!((group.channels[1].radius_phase - 7.64).abs() < 1e-12);
        assert!((group.channels[1].radius_true - 7.5).abs() < 1e-12);
        assert_eq!(group.resonances.len(), 2);
        assert_eq!(group.resonances[1].energy, 101.9);
        assert_eq!(group.resonances[1].widths, vec![0.06, -0.02]);
        assert_eq!(rml.resonance_count(), 2);
    }

    /// ENDF-6 defines PNT=-1 as "do not calculate the penetrability".
    #[test]
    fn accepts_pnt_minus_one_on_photon_pair() {
        let mut fixture = Fixture::basic();
        fixture.pairs[0][7] = -1.0;
        let data = fixture.parse().unwrap();
        assert!(!rml(&data).pairs[0].penetrability);
    }

    #[test]
    fn pnt_zero_takes_the_mt_default() {
        let mut fixture = Fixture::basic();
        fixture.pairs[1][7] = 0.0;
        let data = fixture.parse().unwrap();
        assert!(rml(&data).pairs[1].penetrability, "MT=2 defaults to PNT=1");
    }

    #[test]
    fn rejects_undefined_or_fractional_pnt() {
        let mut fixture = Fixture::basic();
        fixture.pairs[1][7] = 2.0;
        assert_rejects(&fixture, "PNT=2 is not defined");
        fixture.pairs[1][7] = 1.5;
        assert_rejects(&fixture, "PNT=1.5 is not an integer");
    }

    #[test]
    fn rejects_penetrability_on_photon_pair() {
        let mut fixture = Fixture::basic();
        fixture.pairs[0][7] = 1.0;
        assert_rejects(
            &fixture,
            "photon pair (MA=0) cannot calculate the penetrability",
        );
    }

    #[test]
    fn rejects_penetrating_pair_without_mass() {
        let mut fixture = Fixture::basic();
        fixture.pairs[1][1] = 0.0;
        assert_rejects(&fixture, "needs positive masses");
    }

    #[test]
    fn rejects_unsupported_header_flags() {
        let mut fixture = Fixture::basic();
        fixture.krm = 2;
        assert_rejects(&fixture, "KRM=2 is not supported");
        let mut fixture = Fixture::basic();
        fixture.ifg = 1;
        assert_rejects(&fixture, "IFG=1");
        let mut fixture = Fixture::basic();
        fixture.krl = 1;
        assert_rejects(&fixture, "KRL=1");
        let mut fixture = Fixture::basic();
        fixture.nro = true;
        assert_rejects(&fixture, "NRO=1");
        let mut fixture = Fixture::basic();
        fixture.groups[0].kps = 1;
        assert_rejects(&fixture, "KPS=1");
    }

    #[test]
    fn accepts_charged_and_rejects_fission_pairs() {
        let mut fixture = Fixture::basic();
        let mut proton = NEUTRON;
        proton[0] = 0.9986235;
        proton[2] = 1.0;
        proton[3] = 73.0;
        proton[6] = 6.15e5;
        proton[9] = 600.0;
        fixture.pairs.push(proton);
        let data = fixture.parse().unwrap();
        let pair = &rml(&data).pairs[2];
        assert_eq!((pair.charge_a, pair.charge_b, pair.mt), (1.0, 73.0, 600));
        let mut fixture = Fixture::basic();
        let mut fission = PHOTON;
        fission[0] = 1.0;
        fission[9] = 18.0;
        fission[7] = -1.0;
        fixture.pairs.push(fission);
        assert_rejects(&fixture, "fission pair (MT=18)");
    }

    #[test]
    fn shift_factor_needs_an_open_channel() {
        let mut fixture = Fixture::basic();
        fixture.pairs[1][8] = 1.0;
        let data = fixture.parse().unwrap();
        assert!(
            rml(&data).pairs[1].shift,
            "SHF=1 on the incident pair is open everywhere"
        );

        let mut fixture = Fixture::basic();
        let mut inelastic = NEUTRON;
        inelastic[6] = -1.0e3;
        inelastic[8] = 1.0;
        inelastic[9] = 51.0;
        fixture.pairs.push(inelastic);
        assert_rejects(&fixture, "SHF=1 on a channel with threshold");
    }

    /// The eliminated channel is the photon pair's channel at any position.
    #[test]
    fn finds_photon_channel_by_pair_not_position() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].channels.swap(0, 1);
        for r in &mut fixture.groups[0].resonances {
            r.swap(1, 2);
        }
        let data = fixture.parse().unwrap();
        let group = &rml(&data).spin_groups[0];
        assert_eq!(group.channels[1].pair, 0);
        assert_eq!(group.resonances[1].widths, vec![-0.02, 0.06]);
    }

    #[test]
    fn rejects_missing_or_duplicate_photon_and_entrance_channels() {
        let mut fixture = Fixture::basic();
        fixture.pairs[0][7] = -1.0;
        fixture.pairs[0][9] = 103.0;
        assert_rejects(&fixture, "exactly one photon pair");
        let mut fixture = Fixture::basic();
        fixture.groups[0].channels.truncate(1);
        for r in &mut fixture.groups[0].resonances {
            r.truncate(2);
        }
        assert_rejects(&fixture, "has no neutron entrance channel");
        let mut fixture = Fixture::basic();
        fixture.groups[0].channels[0] = fixture.groups[0].channels[1];
        assert_rejects(&fixture, "exactly one channel of the photon pair");
        let mut fixture = Fixture::basic();
        fixture.groups[0].channels[1][0] = 1.0;
        assert_rejects(&fixture, "needs exactly one channel of the photon pair");
        let mut fixture = Fixture::basic();
        fixture.groups[0].channels[1][0] = 3.0;
        assert_rejects(&fixture, "PPI=3 is out of range");
    }

    #[test]
    fn rejects_resonance_on_a_penetrating_threshold() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].resonances[0][0] = 0.0;
        assert_rejects(&fixture, "sits on the threshold");
    }

    #[test]
    fn negative_aj_and_pj_give_negative_parity() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].aj = -1.5;
        let data = fixture.parse().unwrap();
        assert_eq!(
            (
                rml(&data).spin_groups[0].j,
                rml(&data).spin_groups[0].parity
            ),
            (1.5, -1)
        );
        let mut fixture = Fixture::basic();
        fixture.groups[0].aj = 0.0;
        fixture.groups[0].pj = -1.0;
        let data = fixture.parse().unwrap();
        assert_eq!(rml(&data).spin_groups[0].parity, -1);
    }

    #[test]
    fn empty_group_keeps_one_filler_line() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].resonances.clear();
        let data = fixture.parse().unwrap();
        assert!(rml(&data).spin_groups[0].resonances.is_empty());
        fixture.groups[0].nx = Some(2);
        assert_rejects(&fixture, "NX=2");
    }

    /// Six or more channels put each resonance on NCH/6+1 lines.
    #[test]
    fn wide_group_reads_two_lines_per_resonance() {
        let mut fixture = Fixture::basic();
        let mut channels = vec![[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]];
        for l in 0..5 {
            channels.push([2.0, f64::from(l), 0.5, 0.0, 0.7, 0.7]);
        }
        let group = &mut fixture.groups[0];
        group.channels = channels;
        group.resonances = vec![
            vec![5.0, 0.04, 1.0, 2.0, 3.0, 4.0, 5.0],
            vec![9.0, 0.05, 6.0, 7.0, 8.0, 9.0, 10.0],
        ];
        let data = fixture.parse().unwrap();
        let group = &rml(&data).spin_groups[0];
        assert_eq!(group.resonances[1].energy, 9.0);
        assert_eq!(
            group.resonances[1].widths,
            vec![0.05, 6.0, 7.0, 8.0, 9.0, 10.0]
        );
        fixture.groups[0].nx = Some(2);
        assert_rejects(&fixture, "need NX=4");
    }

    #[test]
    fn reads_logarithmic_background() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].background = vec![Background {
            lch: 2,
            lbk: 2,
            list: Some([0.0, 2.0e4, -0.043, 1.0e-6, 2.0e-9, 0.01, 3.0e-7]),
        }];
        let data = fixture.parse().unwrap();
        let background = &rml(&data).spin_groups[0].background;
        assert_eq!(background.len(), 1);
        let b = &background[0];
        assert_eq!(b.channel, 1);
        assert_eq!((b.e_down, b.e_up), (0.0, 2.0e4));
        assert_eq!(
            [b.r_con, b.r_lin, b.r_quad, b.s_con, b.s_lin],
            [-0.043, 1.0e-6, 2.0e-9, 0.01, 3.0e-7]
        );
    }

    #[test]
    fn lbk_zero_adds_no_term() {
        let mut fixture = Fixture::basic();
        fixture.groups[0].background = vec![Background {
            lch: 2,
            lbk: 0,
            list: None,
        }];
        let data = fixture.parse().unwrap();
        assert!(rml(&data).spin_groups[0].background.is_empty());
    }

    #[test]
    fn rejects_unsupported_or_misplaced_background() {
        let background = |lch, lbk, list| {
            let mut fixture = Fixture::basic();
            fixture.groups[0].background = vec![Background { lch, lbk, list }];
            fixture
        };
        let list = Some([0.0, 2.0e4, 0.1, 0.0, 0.0, 0.0, 0.0]);
        assert_rejects(&background(2, 1, None), "LBK=1 background is not supported");
        assert_rejects(&background(2, 3, None), "LBK=3 background is not supported");
        assert_rejects(&background(1, 2, list), "eliminated photon channel");
        assert_rejects(&background(3, 2, list), "LCH=3 is out of range");
        assert_rejects(
            &background(2, 2, Some([1.0, 2.0e4, 0.1, 0.0, 0.0, 0.0, 0.0])),
            "ED < EL, EU > EH",
        );
        let mut fixture = background(2, 0, None);
        fixture.groups[0].background.push(Background {
            lch: 2,
            lbk: 2,
            list,
        });
        assert_rejects(&fixture, "given twice for channel 2");
    }

    fn mf2_fixture(name: &str) -> String {
        let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
            .join("../../tests/data/endf/mf2")
            .join(format!("{name}.endf"));
        std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{path:?}: {e}"))
    }

    /// Each ENDF/B-VIII.1 LRF=7 fixture parses to the end with the counts its
    /// records state (resonances per isotope, from the files' NRS fields).
    #[test]
    fn parses_endf_b_viii_1_lrf7_evaluations() {
        let cases = [
            ("W-182", 3, 312),
            ("W-183", 5, 374),
            ("W-184", 3, 218),
            ("W-186", 3, 175),
            ("Fe-54", 5, 737),
            ("Fe-57", 5, 123),
            ("Si-28", 7, 55),
            ("Si-29", 10, 31),
            ("Si-30", 5, 30),
            ("Rh-103", 5, 555),
            ("Cu-63", 6, 1172),
            ("Cu-65", 6, 969),
            ("V-51", 6, 136),
            ("Sr-88", 7, 443),
            ("Cl-35", 8, 272),
            ("Ca-40", 5, 246),
        ];
        for (name, groups, resonances) in cases {
            let data =
                parse_file2_records(&mf2_fixture(name)).unwrap_or_else(|e| panic!("{name}: {e}"));
            let rml = data
                .ranges
                .iter()
                .find_map(|r| r.r_matrix_limited.as_ref())
                .unwrap_or_else(|| panic!("{name}: no LRF=7 range"));
            assert_eq!(rml.spin_groups.len(), groups, "{name}");
            assert_eq!(rml.resonance_count(), resonances, "{name}");
            assert_eq!(rml.pairs[rml.photon_pair].mt, 102, "{name}");
            assert!(!rml.pairs[rml.photon_pair].penetrability, "{name}");
        }
    }

    /// Sr-88 gives a logarithmic background (KBK=1, LBK=2) in every group.
    #[test]
    fn parses_sr88_backgrounds() {
        let data = parse_file2_records(&mf2_fixture("Sr-88")).unwrap();
        let rml = rml(&data);
        assert_eq!(rml.pairs[0].mt, 102);
        assert!(
            !rml.pairs[0].penetrability,
            "Sr-88 writes PNT=-1 on the photon pair"
        );
        for group in &rml.spin_groups {
            assert_eq!(group.background.len(), 1);
            assert_eq!(group.channels[group.background[0].channel].pair, 1);
        }
        let first = &rml.spin_groups[0].background[0];
        assert_eq!((first.e_down, first.e_up), (0.0, 9.55e5));
        assert_eq!((first.r_con, first.s_con, first.s_lin), (-0.043, 0.01, 0.0));
    }

    /// Cu-63 calculates the shift factor on its neutron pair with B = -1 for l = 1.
    #[test]
    fn parses_cu63_shift_and_boundary() {
        let data = parse_file2_records(&mf2_fixture("Cu-63")).unwrap();
        let rml = rml(&data);
        assert!(rml.pairs[rml.incident_pair].shift);
        let p_wave = rml
            .spin_groups
            .iter()
            .flat_map(|g| &g.channels)
            .find(|c| c.pair == rml.incident_pair && c.l == 1)
            .unwrap();
        assert_eq!(p_wave.boundary, -1.0);
    }

    /// Fe-57's inelastic pair (MT=51, Q = -14.413 keV) is a closed channel
    /// over the VENUS band and carries no shift factor.
    #[test]
    fn parses_fe57_inelastic_pair() {
        let data = parse_file2_records(&mf2_fixture("Fe-57")).unwrap();
        let rml = rml(&data);
        let inelastic = rml.pairs.iter().find(|p| p.mt == 51).unwrap();
        assert_eq!(inelastic.q_value, -14413.0);
        assert!(inelastic.penetrability && !inelastic.shift);
    }
}