nucleide-nuclei 0.7.0

//! Static nuclear reference data: atomic masses (AME2020 plus ENDF-derived
//! isomer masses), natural abundances, radioactive half-lives, decay branches,
//! screening-level cross sections, neutron scattering lengths, mean decay
//! energies, and dose factors.
//!
//! # Provenance
//!
//! - Atomic masses: the [AME2020 atomic mass evaluation][ame] (Huang et al.,
//!   *Chinese Physics C* **45**, 030002/030003, 2021; data courtesy the
//!   IAEA-supported AMDC), condensed into the compact [`crate::data`]
//!   tables below. Isomer rows add the ENDF/B-VIII.0 File-1 MT451 `ELIS`
//!   excitation energy: `m = m_ground + E*/931.49410242 u` (ENDF excitation
//!   energies only; no NUBASE import; isomer tapes with unset `ELIS` carry
//!   the ground-state mass).
//! - Natural abundances: standard isotopic compositions from the
//!   ENDF/B-VIII.0 evaluation, expressed as fractions in 0..1.
//! - Half-lives: ENDF/B-VIII.0 decay evaluations (distributed via the IAEA
//!   and BNL/NNDC), expressed in seconds.
//! - Decay branches (`decay_branches.tsv`): per-branch daughters from the
//!   ENDF/B-VIII.0 decay sublibrary (MF8/MT457 NDK records: RTYP decay-mode
//!   code, RFS daughter state flag, BR branching fraction). RTYP digits
//!   apply in emission order (ENDF-102 §8.4; OpenMC `decay.py` digit table):
//!   beta- gives Z+1, EC/beta+ gives Z-1, alpha gives Z-2/A-4, IT is
//!   unchanged, delayed neutrons/protons subtract the emitted nucleons.
//!   Spontaneous-fission and fission-family branches are dropped (depletion
//!   matrices skip `sf` gains); zero-half-life and stable-flagged tapes
//!   yield no rows, so effectively-stable entries stay absent, matching the
//!   half-life table's stable-absent convention.
//! - Simple cross sections (`simple_xs.tsv`): **total** microscopic cross
//!   sections in barns. Thermal values combine NIST NCNR 2200 m/s bound
//!   scattering/absorption converted to free-atom totals via
//!   `xs·(A/(A+1))² + xs_a` (ENDF File 3 alone is background-only inside
//!   resonance ranges); 14-MeV values are ENDF/B-VII.1 MF3/MT1 interpolated
//!   per the tape INT laws. Screening-level; use evaluated libraries for
//!   transport, never for safety calculations.
//! - Scattering lengths (`scattering_lengths.tsv`): bound coherent /
//!   incoherent lengths in femtometres from the NIST NCNR tabulation of
//!   Sears, *Neutron News* **3**(3) (1992) and its update; complex
//!   (absorbing) lengths enter by real part, incoherent lengths otherwise
//!   derive from the tabulated spin-incoherent cross section via
//!   `b_inc = 10·sqrt(σ_i/4π)` fm. Isotopes the table does not tabulate
//!   (and non-monoisotopic element rows) are absent.
//! - Decay energies (`decay_energy.tsv`): mean *prompt* recoverable energy
//!   per decay in MeV from ENDF/B-VII.1 decay tapes (MF8/MT457 summary
//!   components, paired uncertainties skipped; delayed-neutron kinetic
//!   energy included, neutrinos never appear). Daughter gammas belong to
//!   the daughter row — chain codes must sum members (e.g. Cs137 prompt
//!   plus Ba137_m1 662 keV). Not for spectroscopy or safety use.
//! - Dose factors (`dose_factors.tsv`): external-air/soil, ingestion, and
//!   inhalation factors from the BSD-3 PyNE `dbgen/dosefactors*.csv` tables
//!   (HNF-SD-WM-TI-707 Rev.1 / HNF-5636 App. O; GENII/EPA/DOE are 3 parallel
//!   evaluations). `+D` folds into the parent. Air is EPA-only: GENII/DOE
//!   air rows are `-1` sentinels (PyNE convention); accessors return the
//!   stored `-1` and dose analytics treat negative factors as missing.
//!   Not for safety decisions (upstream PyNE disclaimer).
//!
//! The screening tables are generated by
//! `scripts/gen-nuclear-data.py` (stdlib-only; run `--help` for the
//! upstream download URLs) — never hand-edit a value, regenerate.
//!
//! All tables are vendored as tab-separated text under `src/data/` and
//! embedded with `include_str!`; no runtime dependencies beyond `std`.
//! Parsing happens lazily on first lookup into fixed static maps.
//!
//! Table contents:
//!
//! - `data/ame2020.tsv`: one row per ground-state nuclide (`nucid`,
//!   `mass_u`, `uncertainty_u`), covering all 3 557 nuclides with `Z >= 1`
//!   (the free-neutron row of the source file is dropped), plus one row per
//!   ENDF/B-VIII.0 isomer tape (full state-bearing nucid, 738 rows):
//!   `m = m_ground + ELIS/931.49410242 u` (4 295 rows total).
//! - `data/natural_abundance.tsv`: `GNDS name` → `fraction`, including the
//!   lone naturally occurring isomer Ta180_m1.
//! - `data/half_life.tsv`: `GNDS name` → `half_life_seconds`, for every
//!   radionuclide in the evaluation (stable nuclides are simply absent).
//! - `data/decay_branches.tsv`: `parent_GNDS` → (`progeny_GNDS`, `bf`,
//!   `mode`) per kept branch (5 068 rows over 3 541 parents; SF/fission
//!   branches dropped, so strong SF emitters sum to `1 - BR(SF)`).
//! - `data/simple_xs.tsv`: `GNDS name` → (`thermal_barn`, `fast14mev_barn`)
//!   total cross sections (241 rows; resonance nuclides without a NIST row
//!   and isomers are absent by construction).
//! - `data/scattering_lengths.tsv`: `GNDS name` → (`b_coherent_fm`,
//!   `b_incoherent_fm`) bound scattering lengths (267 rows, NIST-tabulated
//!   isotopes plus monoisotopic element attributions).
//! - `data/decay_energy.tsv`: `GNDS name` → `mev_per_decay` mean prompt
//!   recoverable decay energy (3 557 rows covering every ENDF/B-VII.1 decay
//!   tape with nonzero heat data, isomers with prompt gammas such as
//!   Ba137_m1 included). Basis note: this table stays on ENDF/B-VII.1
//!   while `half_life.tsv` and `decay_branches.tsv` use ENDF/B-VIII.0.
//! - `data/dose_factors.tsv`: (`GNDS name`, `pathway`, `source`) → (`factor`,
//!   `f1`, `lung_model`) dose factors (1 116 rows: 93 folded nuclides × 4
//!   pathways × 3 sources; `f1` set only on ingest rows, `lung_model` only
//!   on inhale rows).
//!
//! Masses are keyed by the canonical [`NuclideId`] nucid layout, including
//! the metastable state: isomers with their own ENDF/B-VIII.0 tape resolve
//! to `m_ground + E*/931.49410242 u`, while isomers without a tape fall back
//! to the ground-state mass (the missing excitation energy is at most MeV
//! against a GeV-scale mass). Q-value helpers stay ground-state-only by
//! design (see [`q_value_neutron_capture`]): they combine tabulated atomic
//! masses per the textbook formulas, and isomer-resolved Q-values would need
//! excitation bookkeeping beyond this screening-level scope.
//!
//! [ame]: https://doi.org/10.1088/1674-1137/abddaf
//!
//! Ownership note: this module + its data files only — do NOT edit `lib.rs`.

use std::collections::BTreeMap;
use std::sync::OnceLock;

use crate::NuclideId;

const AME2020_TSV: &str = include_str!("data/ame2020.tsv");
const NATURAL_ABUNDANCE_TSV: &str = include_str!("data/natural_abundance.tsv");
const HALF_LIFE_TSV: &str = include_str!("data/half_life.tsv");
const SIMPLE_XS_TSV: &str = include_str!("data/simple_xs.tsv");
const SCATTERING_LENGTHS_TSV: &str = include_str!("data/scattering_lengths.tsv");
const DECAY_ENERGY_TSV: &str = include_str!("data/decay_energy.tsv");
const DECAY_BRANCHES_TSV: &str = include_str!("data/decay_branches.tsv");
const DOSE_FACTORS_TSV: &str = include_str!("data/dose_factors.tsv");

static MASSES: OnceLock<BTreeMap<u32, f64>> = OnceLock::new();
static ABUNDANCES: OnceLock<BTreeMap<u32, f64>> = OnceLock::new();
static HALF_LIVES: OnceLock<BTreeMap<u32, f64>> = OnceLock::new();
static SIMPLE_XS: OnceLock<BTreeMap<u32, (f64, f64)>> = OnceLock::new();
static SCATTERING_LENGTHS: OnceLock<BTreeMap<u32, (f64, f64)>> = OnceLock::new();
static DECAY_ENERGIES: OnceLock<BTreeMap<u32, f64>> = OnceLock::new();
static DECAY_BRANCHES: OnceLock<BTreeMap<u32, Vec<DecayBranch>>> = OnceLock::new();
static DOSE_FACTORS: OnceLock<BTreeMap<(u32, DosePathway, DoseSource), DoseEntry>> =
    OnceLock::new();

/// MeV per unified atomic mass unit `c²` (2022 CODATA consistent with
/// AME2020 usage).
pub const MEV_PER_U: f64 = 931.494_102_42;

/// Free-neutron atomic mass in u (AME2020: 1.00866491595 u).
pub const NEUTRON_MASS_U: f64 = 1.008_664_915_95;

/// Helium-4 atomic mass in u (AME2020), for alpha-decay Q-values.
pub const HELIUM4_MASS_U: f64 = 4.002_603_254_13;

/// The free neutron's mass, in u.
///
/// See [`NEUTRON_MASS_U`].
pub const fn neutron_mass_u() -> f64 {
    NEUTRON_MASS_U
}

/// Parse `nucid \t mass_u [\t uncertainty_u]` rows, skipping comment lines.
fn parse_masses(tsv: &str) -> BTreeMap<u32, f64> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let nucid = cols.next()?.parse().ok()?;
            let mass = cols.next()?.parse().ok()?;
            Some((nucid, mass))
        })
        .collect()
}

/// Parse `GNDS name \t fraction` rows into nucid-keyed fractions.
///
/// Malformed rows are skipped silently, matching [`parse_masses`]; the
/// embedded files are compile-time constants and a corrupt row should not
/// panic library consumers.
fn parse_abundances(tsv: &str) -> BTreeMap<u32, f64> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let fraction: f64 = cols.next()?.parse().ok()?;
            let nucid = NuclideId::from_name(name).ok()?;
            Some((nucid.nucid(), fraction))
        })
        .collect()
}

/// Parse `GNDS name \t half_life_seconds` rows into nucid-keyed seconds.
///
/// Malformed rows are skipped silently, matching [`parse_masses`]; the
/// embedded files are compile-time constants and a corrupt row should not
/// panic library consumers.
fn parse_half_lives(tsv: &str) -> BTreeMap<u32, f64> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let seconds: f64 = cols.next()?.parse().ok()?;
            let nucid = NuclideId::from_name(name).ok()?;
            Some((nucid.nucid(), seconds))
        })
        .collect()
}

fn masses() -> &'static BTreeMap<u32, f64> {
    MASSES.get_or_init(|| parse_masses(AME2020_TSV))
}

fn abundances() -> &'static BTreeMap<u32, f64> {
    ABUNDANCES.get_or_init(|| parse_abundances(NATURAL_ABUNDANCE_TSV))
}

fn half_lives() -> &'static BTreeMap<u32, f64> {
    HALF_LIVES.get_or_init(|| parse_half_lives(HALF_LIFE_TSV))
}

/// Parse `GNDS name \t thermal_barn \t fast14mev_barn` rows into
/// nucid-keyed `(thermal, fast)` pairs.
///
/// Malformed rows are skipped silently, matching [`parse_masses`].
fn parse_simple_xs(tsv: &str) -> BTreeMap<u32, (f64, f64)> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let thermal: f64 = cols.next()?.parse().ok()?;
            let fast: f64 = cols.next()?.parse().ok()?;
            let nucid = NuclideId::from_name(name).ok()?;
            Some((nucid.nucid(), (thermal, fast)))
        })
        .collect()
}

/// Parse `GNDS name \t b_coherent_fm \t b_incoherent_fm` rows into
/// nucid-keyed `(coherent, incoherent)` pairs.
///
/// Malformed rows are skipped silently, matching [`parse_masses`].
fn parse_scattering_lengths(tsv: &str) -> BTreeMap<u32, (f64, f64)> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let coherent: f64 = cols.next()?.parse().ok()?;
            let incoherent: f64 = cols.next()?.parse().ok()?;
            let nucid = NuclideId::from_name(name).ok()?;
            Some((nucid.nucid(), (coherent, incoherent)))
        })
        .collect()
}

/// Parse `GNDS name \t mev_per_decay` rows into nucid-keyed MeV.
///
/// Malformed rows are skipped silently, matching [`parse_masses`].
fn parse_decay_energies(tsv: &str) -> BTreeMap<u32, f64> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let mev: f64 = cols.next()?.parse().ok()?;
            let nucid = NuclideId::from_name(name).ok()?;
            Some((nucid.nucid(), mev))
        })
        .collect()
}

fn simple_xs_map() -> &'static BTreeMap<u32, (f64, f64)> {
    SIMPLE_XS.get_or_init(|| parse_simple_xs(SIMPLE_XS_TSV))
}

fn scattering_length_map() -> &'static BTreeMap<u32, (f64, f64)> {
    SCATTERING_LENGTHS.get_or_init(|| parse_scattering_lengths(SCATTERING_LENGTHS_TSV))
}

fn decay_energy_map() -> &'static BTreeMap<u32, f64> {
    DECAY_ENERGIES.get_or_init(|| parse_decay_energies(DECAY_ENERGY_TSV))
}

/// The full AME2020 ground-state mass table, keyed by nucid (in u).
pub fn mass_table() -> &'static BTreeMap<u32, f64> {
    masses()
}

/// The full natural-abundance table, keyed by nucid (fraction in 0..1).
pub fn abundance_table() -> &'static BTreeMap<u32, f64> {
    abundances()
}

/// The full half-life table in seconds, keyed by nucid.
pub fn half_life_table() -> &'static BTreeMap<u32, f64> {
    half_lives()
}

/// The full simple cross-section table, keyed by nucid:
/// `(thermal_barn, fast14mev_barn)` total cross sections.
pub fn simple_xs_table() -> &'static BTreeMap<u32, (f64, f64)> {
    simple_xs_map()
}

/// The full bound scattering-length table, keyed by nucid:
/// `(b_coherent_fm, b_incoherent_fm)`.
pub fn scattering_length_table() -> &'static BTreeMap<u32, (f64, f64)> {
    scattering_length_map()
}

/// The full mean-recoverable-decay-energy table in MeV, keyed by nucid.
pub fn decay_energy_table() -> &'static BTreeMap<u32, f64> {
    decay_energy_map()
}

/// Total microscopic cross sections `(thermal_barn, fast14mev_barn)` for the
/// given `nucid`.
///
/// Screening-level values (see [`simple_xs_table`]); returns `None` for
/// nuclides outside the table — notably resonance absorbers without a NIST
/// row (Cs137, Co60, I135, Xe135) and isomers.
pub fn simple_xs(nucid: u32) -> Option<(f64, f64)> {
    simple_xs_map().get(&nucid).copied()
}

/// Approximate total microscopic cross sections of a named nuclide (see
/// [`NuclideId::from_name`]): `(thermal_barn, fast14mev_barn)`.
pub fn simple_xs_by_name(name: &str) -> Option<(f64, f64)> {
    simple_xs(NuclideId::from_name(name).ok()?.nucid())
}

/// Bound neutron scattering lengths `(b_coherent_fm, b_incoherent_fm)` for
/// the given `nucid`.
///
/// Returns `None` for nuclides outside the NIST-tabulated set (see
/// [`scattering_length_table`]).
pub fn scattering_length(nucid: u32) -> Option<(f64, f64)> {
    scattering_length_map().get(&nucid).copied()
}

/// Bound neutron scattering lengths of a named nuclide (see
/// [`NuclideId::from_name`]): `(b_coherent_fm, b_incoherent_fm)`.
pub fn scattering_length_by_name(name: &str) -> Option<(f64, f64)> {
    scattering_length(NuclideId::from_name(name).ok()?.nucid())
}

/// Mean prompt recoverable decay energy per decay of the given `nucid`, in
/// MeV (see [`decay_energy_table`]).
///
/// Returns `None` for nuclides outside the table, including stable nuclides
/// such as O16 and Fe56.
pub fn decay_energy_mev(nucid: u32) -> Option<f64> {
    decay_energy_map().get(&nucid).copied()
}

/// Mean recoverable decay energy per decay of a named nuclide (see
/// [`NuclideId::from_name`]), in MeV.
pub fn decay_energy_mev_by_name(name: &str) -> Option<f64> {
    decay_energy_mev(NuclideId::from_name(name).ok()?.nucid())
}

// ---------------------------------------------------------------------------
// Decay branches
// ---------------------------------------------------------------------------

/// One decay mode token, matching the depletion-chain vocabulary (`beta-`,
/// `ec/beta+`, `alpha`, `IT`, `sf`) plus the direct nucleon-emission modes
/// (`n`, `p`) for particle-unbound light nuclei.
///
/// EC-only and beta+/EC evaluations both fold into [`Self::EcBetaPlus`].
/// `sf` never appears in the vendored table (fission branches are dropped
/// at generation) but parses for robustness.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub enum DecayBranchMode {
    /// Negative beta decay (Z+1 daughter, incl. delayed-neutron branches).
    BetaMinus,
    /// Electron capture and/or positron emission (Z-1 daughter).
    EcBetaPlus,
    /// Alpha decay (Z-2/A-4 daughter).
    Alpha,
    /// Isomeric transition (same Z/A daughter).
    It,
    /// Spontaneous fission (dropped at generation; parses only).
    Sf,
    /// Direct neutron emission (A-1 daughter).
    Neutron,
    /// Direct proton emission (Z-1/A-1 daughter).
    Proton,
}

impl DecayBranchMode {
    /// Canonical table token (`beta-`/`ec/beta+`/`alpha`/`IT`/`sf`/`n`/`p`).
    pub fn as_str(self) -> &'static str {
        match self {
            Self::BetaMinus => "beta-",
            Self::EcBetaPlus => "ec/beta+",
            Self::Alpha => "alpha",
            Self::It => "IT",
            Self::Sf => "sf",
            Self::Neutron => "n",
            Self::Proton => "p",
        }
    }

    /// Parse a mode token (case-insensitive).
    pub fn parse(s: &str) -> Option<Self> {
        match s.trim().to_ascii_lowercase().as_str() {
            "beta-" | "beta" | "b-" => Some(Self::BetaMinus),
            "ec/beta+" | "ec" | "beta+" => Some(Self::EcBetaPlus),
            "alpha" | "a" => Some(Self::Alpha),
            "it" => Some(Self::It),
            "sf" => Some(Self::Sf),
            "n" => Some(Self::Neutron),
            "p" => Some(Self::Proton),
            _ => None,
        }
    }
}

/// One evaluated decay branch: daughter nuclide, branching fraction, mode.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DecayBranch {
    /// Daughter nucid (full state-bearing id; RFS sets the state).
    pub progeny: u32,
    /// Branching fraction in 0..1 (evaluated value, verbatim).
    pub branching_fraction: f64,
    /// Decay mode token (initial event for multi-particle branches).
    pub mode: DecayBranchMode,
}

/// Parse `parent_GNDS \t progeny_GNDS \t bf \t mode` rows into
/// parent-nucid-keyed branch lists.
///
/// Malformed rows are skipped silently, matching [`parse_masses`].
fn parse_decay_branches(tsv: &str) -> BTreeMap<u32, Vec<DecayBranch>> {
    let mut map: BTreeMap<u32, Vec<DecayBranch>> = BTreeMap::new();
    for line in tsv
        .lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
    {
        let mut cols = line.split('\t');
        let (Some(parent), Some(progeny), Some(bf), Some(mode)) =
            (cols.next(), cols.next(), cols.next(), cols.next())
        else {
            continue;
        };
        let (Ok(p), Ok(d), Ok(b), Some(m)) = (
            NuclideId::from_name(parent).map(|id| id.nucid()),
            NuclideId::from_name(progeny).map(|id| id.nucid()),
            bf.parse::<f64>(),
            DecayBranchMode::parse(mode),
        ) else {
            continue;
        };
        map.entry(p).or_default().push(DecayBranch {
            progeny: d,
            branching_fraction: b,
            mode: m,
        });
    }
    for branches in map.values_mut() {
        branches.sort_by_key(|b| (b.progeny, b.mode));
    }
    map
}

fn decay_branch_map() -> &'static BTreeMap<u32, Vec<DecayBranch>> {
    DECAY_BRANCHES.get_or_init(|| parse_decay_branches(DECAY_BRANCHES_TSV))
}

/// The full decay-branch table, keyed by parent nucid.
///
/// 5 068 rows over 3 541 parents from ENDF/B-VIII.0; SF/fission branches
/// are dropped, so strong SF emitters sum to `1 - BR(SF)`.
pub fn decay_branch_table() -> &'static BTreeMap<u32, Vec<DecayBranch>> {
    decay_branch_map()
}

/// Evaluated decay branches of the given parent `nucid`.
///
/// Returns `None` for nuclides with no branch rows (stable nuclides and
/// zero-half-life evaluation dummies such as Te123).
pub fn decay_branches(nucid: u32) -> Option<Vec<DecayBranch>> {
    decay_branch_map().get(&nucid).cloned()
}

/// Evaluated decay branches of a named nuclide (see [`NuclideId::from_name`]).
pub fn decay_branches_by_name(name: &str) -> Option<Vec<DecayBranch>> {
    decay_branches(NuclideId::from_name(name).ok()?.nucid())
}

/// Branching fraction from `parent` to `progeny` (both nucids), if tabulated.
pub fn branching_fraction(parent: u32, progeny: u32) -> Option<f64> {
    decay_branch_map()
        .get(&parent)?
        .iter()
        .find(|b| b.progeny == progeny)
        .map(|b| b.branching_fraction)
}

/// Branching fraction from `parent` to `progeny` (GNDS names), if tabulated.
pub fn branching_fraction_by_name(parent: &str, progeny: &str) -> Option<f64> {
    branching_fraction(
        NuclideId::from_name(parent).ok()?.nucid(),
        NuclideId::from_name(progeny).ok()?.nucid(),
    )
}

// ---------------------------------------------------------------------------
// Dose factors
// ---------------------------------------------------------------------------

/// Dose pathway: external air/soil, ingestion, or inhalation.
///
/// `Air` is `ext_air` in PyNE (`mrem/h per Ci/m^3`); `Soil` is `ext_soil`
/// (`mrem/h per Ci/m^2`, 15 cm slab); `Ingest`/`Inhale` are `mrem/pCi`.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub enum DosePathway {
    /// External dose from air submersion.
    Air,
    /// External dose from a 15 cm soil slab.
    Soil,
    /// Committed dose from ingestion.
    Ingest,
    /// Committed dose from inhalation.
    Inhale,
}

impl DosePathway {
    /// Canonical lowercase name (`air`/`soil`/`ingest`/`inhale`).
    pub fn as_str(self) -> &'static str {
        match self {
            Self::Air => "air",
            Self::Soil => "soil",
            Self::Ingest => "ingest",
            Self::Inhale => "inhale",
        }
    }

    /// Parse a pathway name (case-insensitive; `ext_air`/`ext_soil` aliases
    /// accepted for PyNE compatibility).
    pub fn parse(s: &str) -> Option<Self> {
        match s.trim().to_ascii_lowercase().as_str() {
            "air" | "ext_air" | "ext-air" => Some(Self::Air),
            "soil" | "ext_soil" | "ext-soil" => Some(Self::Soil),
            "ingest" | "ingestion" => Some(Self::Ingest),
            "inhale" | "inhalation" => Some(Self::Inhale),
            _ => None,
        }
    }
}

/// Dose-factor evaluation source: the three parallel HNF-5636 evaluations.
///
/// PyNE source ids are 0 = EPA (default), 1 = DOE, 2 = GENII.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub enum DoseSource {
    /// EPA evaluation (PyNE source id 0, the default).
    Epa,
    /// DOE evaluation (PyNE source id 1).
    Doe,
    /// GENII evaluation (PyNE source id 2).
    Genii,
}

impl DoseSource {
    /// Canonical uppercase name (`EPA`/`DOE`/`GENII`).
    pub fn as_str(self) -> &'static str {
        match self {
            Self::Epa => "EPA",
            Self::Doe => "DOE",
            Self::Genii => "GENII",
        }
    }

    /// Parse a source name (case-insensitive).
    pub fn parse(s: &str) -> Option<Self> {
        match s.trim().to_ascii_uppercase().as_str() {
            "EPA" => Some(Self::Epa),
            "DOE" => Some(Self::Doe),
            "GENII" => Some(Self::Genii),
            _ => None,
        }
    }

    /// PyNE source id (0 = EPA default, 1 = DOE, 2 = GENII).
    pub fn to_int(self) -> u8 {
        match self {
            Self::Epa => 0,
            Self::Doe => 1,
            Self::Genii => 2,
        }
    }

    /// Inverse of [`DoseSource::to_int`]; anything else is `None`.
    pub fn from_int(v: u8) -> Option<Self> {
        match v {
            0 => Some(Self::Epa),
            1 => Some(Self::Doe),
            2 => Some(Self::Genii),
            _ => None,
        }
    }
}

/// One dose-factor row: the table factor plus aux columns.
///
/// `factor` is the raw table value (`-1` for GENII/DOE air, matching PyNE's
/// missing-air sentinel); `f1` (fraction to body fluids) is set only on
/// ingest rows; `lung_model` (`D`/`W`/`Y`/`V`/`O`) only on inhale rows.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DoseEntry {
    /// Raw table factor (see [`DosePathway`] for units).
    pub factor: f64,
    /// `f1` aux column (ingest rows only).
    pub f1: Option<f64>,
    /// Lung-model aux column (inhale rows only).
    pub lung_model: Option<char>,
}

/// Parse `GNDS name \t pathway \t source \t factor [\t f1 [\t lung]]` rows.
///
/// Malformed rows are skipped silently, matching [`parse_masses`].
fn parse_dose_factors(tsv: &str) -> BTreeMap<(u32, DosePathway, DoseSource), DoseEntry> {
    tsv.lines()
        .filter(|line| !line.is_empty() && !line.starts_with('#'))
        .filter_map(|line| {
            let mut cols = line.split('\t');
            let name = cols.next()?;
            let pathway = DosePathway::parse(cols.next()?)?;
            let source = DoseSource::parse(cols.next()?)?;
            let factor: f64 = cols.next()?.parse().ok()?;
            // Trailing aux columns may be absent after `strip`-style writes;
            // missing means None.
            let f1_txt = cols.next().unwrap_or("");
            let lung_txt = cols.next().unwrap_or("");
            let f1 = if f1_txt.trim().is_empty() {
                None
            } else {
                Some(f1_txt.parse().ok()?)
            };
            let lung_model = {
                let t = lung_txt.trim();
                if t.is_empty() {
                    None
                } else {
                    t.chars().next()
                }
            };
            let nucid = NuclideId::from_name(name).ok()?;
            Some((
                (nucid.nucid(), pathway, source),
                DoseEntry {
                    factor,
                    f1,
                    lung_model,
                },
            ))
        })
        .collect()
}

fn dose_factor_map() -> &'static BTreeMap<(u32, DosePathway, DoseSource), DoseEntry> {
    DOSE_FACTORS.get_or_init(|| parse_dose_factors(DOSE_FACTORS_TSV))
}

/// The full dose-factor table, keyed by `(nucid, pathway, source)`.
///
/// 1 116 rows (93 folded nuclides × 4 pathways × 3 sources).
pub fn dose_table() -> &'static BTreeMap<(u32, DosePathway, DoseSource), DoseEntry> {
    dose_factor_map()
}

/// Full dose-factor entry for the given `nucid`, `pathway`, and `source`.
///
/// Returns the stored `-1` sentinel for GENII/DOE air (PyNE convention for
/// missing air data); `None` only when the nuclide has no row at all.
pub fn dose_entry(nucid: u32, pathway: DosePathway, source: DoseSource) -> Option<DoseEntry> {
    dose_factor_map().get(&(nucid, pathway, source)).copied()
}

/// Raw dose factor for the given `nucid`, `pathway`, and `source`.
///
/// See [`dose_entry`]: GENII/DOE air resolve to `Some(-1.0)` (missing-air
/// sentinel); unknown nuclides resolve to `None` (PyNE's `-1` unknown
/// sentinel maps to `None` here).
pub fn dose_factor(nucid: u32, pathway: DosePathway, source: DoseSource) -> Option<f64> {
    dose_entry(nucid, pathway, source).map(|e| e.factor)
}

/// Raw dose factor for a named nuclide (see [`NuclideId::from_name`]).
pub fn dose_factor_by_name(name: &str, pathway: DosePathway, source: DoseSource) -> Option<f64> {
    dose_factor(NuclideId::from_name(name).ok()?.nucid(), pathway, source)
}

/// `f1` (fraction to body fluids) for an ingest row, if present.
pub fn dose_f1(nucid: u32, source: DoseSource) -> Option<f64> {
    dose_entry(nucid, DosePathway::Ingest, source)?.f1
}

/// `f1` for a named nuclide's ingest row.
pub fn dose_f1_by_name(name: &str, source: DoseSource) -> Option<f64> {
    dose_f1(NuclideId::from_name(name).ok()?.nucid(), source)
}

/// Lung model (`D`/`W`/`Y`/`V`/`O`) for an inhale row, if present.
pub fn dose_lung_model(nucid: u32, source: DoseSource) -> Option<char> {
    dose_entry(nucid, DosePathway::Inhale, source)?.lung_model
}

/// Lung model for a named nuclide's inhale row.
pub fn dose_lung_model_by_name(name: &str, source: DoseSource) -> Option<char> {
    dose_lung_model(NuclideId::from_name(name).ok()?.nucid(), source)
}

/// Atomic mass of the nuclide with the given `nucid`, in u.
///
/// Ground states resolve to the AME2020 value. Metastable states resolve to
/// their ENDF-derived isomer mass (`m_ground + E*/931.49410242 u`) when the
/// isomer has its own tape, and fall back to the ground-state mass otherwise
/// (the missing excitation energy is at most MeV against a GeV-scale mass).
/// Returns `None` for unknown nuclides and non-nuclide ids (`Z = 0`).
pub fn atomic_mass(nucid: u32) -> Option<f64> {
    if let Some(mass) = masses().get(&nucid) {
        return Some(*mass);
    }
    let id = NuclideId::from_nucid(nucid);
    if id.state() != 0 {
        let ground = (id.z() * 1000 + id.a()) * 10_000;
        return masses().get(&ground).copied();
    }
    None
}

/// Atomic mass of a named nuclide (see [`NuclideId::from_name`]), in u.
pub fn atomic_mass_by_name(name: &str) -> Option<f64> {
    atomic_mass(NuclideId::from_name(name).ok()?.nucid())
}

/// Natural abundance of the nuclide with the given `nucid`, as a fraction
/// in 0..1.
///
/// Only naturally occurring nuclides have entries, including the isomer
/// Ta180_m1; everything else returns `None`.
pub fn natural_abundance(nucid: u32) -> Option<f64> {
    abundances().get(&nucid).copied()
}

/// Natural abundance of a named nuclide, as a fraction in 0..1.
pub fn natural_abundance_by_name(name: &str) -> Option<f64> {
    natural_abundance(NuclideId::from_name(name).ok()?.nucid())
}

/// Radioactive half-life of the nuclide with the given `nucid`, in seconds.
///
/// Values come from the ENDF/B-VIII.0-derived evaluation; stable
/// nuclides and unknown ids have no entry and return `None`.
pub fn half_life(nucid: u32) -> Option<f64> {
    half_lives().get(&nucid).copied()
}

/// Half-life in seconds of a named nuclide (see [`NuclideId::from_name`]).
pub fn half_life_by_name(name: &str) -> Option<f64> {
    half_life(NuclideId::from_name(name).ok()?.nucid())
}

/// Decay constant λ = ln(2) / t½ of the given `nucid`, in inverse seconds.
///
/// Computed from [`half_life`]; `None` wherever the half-life is unknown.
pub fn decay_constant(nucid: u32) -> Option<f64> {
    half_life(nucid).map(|t_half| std::f64::consts::LN_2 / t_half)
}

/// Decay constant λ = ln(2) / t½ (inverse seconds) of a named nuclide.
pub fn decay_constant_by_name(name: &str) -> Option<f64> {
    decay_constant(NuclideId::from_name(name).ok()?.nucid())
}

/// Q-value of neutron radiative capture X(n,γ)X', in MeV.
///
/// Ground-state-only by design: the formula combines tabulated atomic
/// masses per the textbook definition, and isomer-resolved Q-values would
/// need explicit excitation-energy bookkeeping beyond this
/// screening-level scope (metastable-state and non-nuclide inputs return
/// `None` even though [`atomic_mass`] now resolves isomers).
///
/// From atomic masses (AME2020, in u):
///
/// ```text
/// Q = [m(X + A) + m_n − m(X' = Z, A+1)] · c²   (MeV/u × 931.49410242)
/// ```
///
/// using the AME2020 free-neutron mass ([`NEUTRON_MASS_U`]). Electron
/// binding differences are neglected (standard for capture Q-values).
/// Returns `None` when either ground-state mass is missing, or for
/// metastable-state or non-nuclide inputs.
pub fn q_value_neutron_capture(nucid: u32) -> Option<f64> {
    let id = NuclideId::from_nucid(nucid);
    if id.state() != 0 || id.z() == 0 {
        return None;
    }
    let product = (id.z() * 1000 + id.a() + 1) * 10_000;
    let q_u = atomic_mass(nucid)? + NEUTRON_MASS_U - atomic_mass(product)?;
    Some(q_u * MEV_PER_U)
}

/// Q-value of neutron capture on a named nuclide, in MeV.
pub fn q_value_neutron_capture_by_name(name: &str) -> Option<f64> {
    q_value_neutron_capture(NuclideId::from_name(name).ok()?.nucid())
}

/// Q-value of alpha decay X → Y + He4, in MeV.
///
/// Ground-state-only by design, like [`q_value_neutron_capture`]:
/// metastable and non-nuclide inputs return `None`.
///
/// From atomic masses (AME2020, in u):
///
/// ```text
/// Q = [m(Z,A) − m(Z−2,A−4) − m(He4)] · c²   (MeV/u × 931.49410242)
/// ```
///
/// with [`HELIUM4_MASS_U`] = 4.00260325413 u. Using neutral atomic masses,
/// the two released electrons cancel exactly to first order. Returns `None`
/// when the daughter mass is absent from the table, or for metastable or
/// non-nuclide inputs.
pub fn q_value_alpha(nucid: u32) -> Option<f64> {
    let id = NuclideId::from_nucid(nucid);
    if id.state() != 0 || id.z() <= 2 || id.a() <= 4 {
        return None;
    }
    let daughter = ((id.z() - 2) * 1000 + (id.a() - 4)) * 10_000;
    let q_u = atomic_mass(nucid)? - atomic_mass(daughter)? - HELIUM4_MASS_U;
    Some(q_u * MEV_PER_U)
}

/// Q-value of alpha decay on a named nuclide, in MeV.
pub fn q_value_alpha_by_name(name: &str) -> Option<f64> {
    q_value_alpha(NuclideId::from_name(name).ok()?.nucid())
}

/// Zero-sized façade over this module's lookups.
///
/// Standalone today; intended to satisfy the material crate's future
/// `MassProvider` trait once that integration lands (kept free of cross-crate
/// coupling for now).
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct AmeMasses;

impl AmeMasses {
    /// See [`atomic_mass`].
    pub fn atomic_mass(&self, nucid: u32) -> Option<f64> {
        atomic_mass(nucid)
    }

    /// See [`atomic_mass_by_name`].
    pub fn atomic_mass_by_name(&self, name: &str) -> Option<f64> {
        atomic_mass_by_name(name)
    }

    /// See [`natural_abundance`].
    pub fn natural_abundance(&self, nucid: u32) -> Option<f64> {
        natural_abundance(nucid)
    }

    /// See [`natural_abundance_by_name`].
    pub fn natural_abundance_by_name(&self, name: &str) -> Option<f64> {
        natural_abundance_by_name(name)
    }
}

/// Zero-sized façade over the decay-data lookups ([`half_life`],
/// [`decay_constant`], [`decay_energy_mev`], [`decay_branches`],
/// [`branching_fraction`]).
///
/// Mirrors [`AmeMasses`]; standalone today, kept free of cross-crate
/// coupling for future provider-trait integration.
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct DecayData;

impl DecayData {
    /// See [`half_life`].
    pub fn half_life(&self, nucid: u32) -> Option<f64> {
        half_life(nucid)
    }

    /// See [`half_life_by_name`].
    pub fn half_life_by_name(&self, name: &str) -> Option<f64> {
        half_life_by_name(name)
    }

    /// See [`decay_constant`].
    pub fn decay_constant(&self, nucid: u32) -> Option<f64> {
        decay_constant(nucid)
    }

    /// See [`decay_constant_by_name`].
    pub fn decay_constant_by_name(&self, name: &str) -> Option<f64> {
        decay_constant_by_name(name)
    }

    /// See [`decay_energy_mev`].
    pub fn decay_energy_mev(&self, nucid: u32) -> Option<f64> {
        decay_energy_mev(nucid)
    }

    /// See [`decay_energy_mev_by_name`].
    pub fn decay_energy_mev_by_name(&self, name: &str) -> Option<f64> {
        decay_energy_mev_by_name(name)
    }

    /// See [`decay_branches`].
    pub fn decay_branches(&self, nucid: u32) -> Option<Vec<DecayBranch>> {
        decay_branches(nucid)
    }

    /// See [`decay_branches_by_name`].
    pub fn decay_branches_by_name(&self, name: &str) -> Option<Vec<DecayBranch>> {
        decay_branches_by_name(name)
    }

    /// See [`branching_fraction`].
    pub fn branching_fraction(&self, parent: u32, progeny: u32) -> Option<f64> {
        branching_fraction(parent, progeny)
    }

    /// See [`branching_fraction_by_name`].
    pub fn branching_fraction_by_name(&self, parent: &str, progeny: &str) -> Option<f64> {
        branching_fraction_by_name(parent, progeny)
    }
}

/// Zero-sized façade over the dose-factor lookups ([`dose_factor`],
/// [`dose_f1`], [`dose_lung_model`]).
///
/// Mirrors [`DecayData`]; the material crate implements its `DoseProvider`
/// trait for this façade so analytics get dose factors with no circular deps.
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct DoseData;

impl DoseData {
    /// See [`dose_entry`].
    pub fn dose_entry(
        &self,
        nucid: u32,
        pathway: DosePathway,
        source: DoseSource,
    ) -> Option<DoseEntry> {
        dose_entry(nucid, pathway, source)
    }

    /// See [`dose_factor`].
    pub fn dose_factor(&self, nucid: u32, pathway: DosePathway, source: DoseSource) -> Option<f64> {
        dose_factor(nucid, pathway, source)
    }

    /// See [`dose_factor_by_name`].
    pub fn dose_factor_by_name(
        &self,
        name: &str,
        pathway: DosePathway,
        source: DoseSource,
    ) -> Option<f64> {
        dose_factor_by_name(name, pathway, source)
    }

    /// See [`dose_f1`].
    pub fn dose_f1(&self, nucid: u32, source: DoseSource) -> Option<f64> {
        dose_f1(nucid, source)
    }

    /// See [`dose_f1_by_name`].
    pub fn dose_f1_by_name(&self, name: &str, source: DoseSource) -> Option<f64> {
        dose_f1_by_name(name, source)
    }

    /// See [`dose_lung_model`].
    pub fn dose_lung_model(&self, nucid: u32, source: DoseSource) -> Option<char> {
        dose_lung_model(nucid, source)
    }

    /// See [`dose_lung_model_by_name`].
    pub fn dose_lung_model_by_name(&self, name: &str, source: DoseSource) -> Option<char> {
        dose_lung_model_by_name(name, source)
    }
}

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

    const H1: u32 = 10_010_000;
    const O16: u32 = 80_160_000;
    const FE56: u32 = 260_560_000;
    const U235: u32 = 922_350_000;

    #[test]
    fn h1_exact_ame2020_value() {
        assert_eq!(atomic_mass(H1), Some(1.007825031_898));
        assert_eq!(atomic_mass_by_name("H1"), Some(1.007825031_898));
    }

    #[test]
    fn heavy_nuclide_spot_values() {
        assert_eq!(atomic_mass(U235), Some(235.043_928_117));
        assert_eq!(atomic_mass(FE56), Some(55.934_935_537));
        assert_eq!(atomic_mass(O16), Some(15.994_914_619_26));
    }

    #[test]
    fn non_nuclides_and_unknown_ids_return_none() {
        // Free neutron row (Z = 0) is excluded from the vendored table.
        assert_eq!(atomic_mass(10_000), None);
        assert_eq!(atomic_mass(999_999_999), None);
        // Isomers with their own ENDF tape resolve to the isomer mass;
        // isomers without one fall back to the ground-state mass.
        let ba_m1 = NuclideId::from_name("Ba137_m1").unwrap().nucid();
        assert!(atomic_mass(ba_m1).unwrap() > atomic_mass(561_370_000).unwrap());
        assert!(atomic_mass_by_name("U235_m1").unwrap() > atomic_mass(U235).unwrap());
        assert_eq!(
            atomic_mass_by_name("Pm137_m1"),
            atomic_mass_by_name("Pm137")
        );
    }

    #[test]
    fn by_name_agrees_with_nucid_lookup() {
        for name in ["H1", "O16", "Fe56", "U235", "Og294"] {
            let nucid = NuclideId::from_name(name).unwrap().nucid();
            assert_eq!(atomic_mass_by_name(name), atomic_mass(nucid), "{name}");
        }
    }

    #[test]
    fn natural_abundance_spot_values() {
        assert_eq!(natural_abundance(U235), Some(0.007_204));
        assert_eq!(natural_abundance_by_name("U235"), Some(0.007_204));
        assert_eq!(natural_abundance(O16), Some(0.997_620_6));
        assert_eq!(natural_abundance_by_name("H1"), Some(0.999_844_26));
        // The lone naturally occurring isomer.
        assert_eq!(natural_abundance_by_name("Ta180_m1"), Some(0.000_120_1));
    }

    #[test]
    fn natural_abundance_unknown_returns_none() {
        assert_eq!(natural_abundance(999_999_999), None);
        assert_eq!(natural_abundance_by_name("C14"), None);
        assert_eq!(natural_abundance_by_name("Xx999"), None);
    }

    #[test]
    fn abundances_sum_to_one_per_element() {
        let mut totals = [0.0_f64; 119];
        for (nucid, frac) in abundance_table() {
            totals[NuclideId::from_nucid(*nucid).z() as usize] += frac;
        }
        for (z, total) in totals.iter().enumerate() {
            if *total > 0.0 {
                assert!(
                    (total - 1.0).abs() < 1e-6,
                    "Z={z} abundances sum to {total}"
                );
            }
        }
    }

    #[test]
    fn mass_sanity_sweep() {
        let table = mass_table();
        for (nucid, mass) in table {
            let id = NuclideId::from_nucid(*nucid);
            // Ground states plus ENDF-taped isomers (state 1-9); the
            // excitation energy is at most MeV against a GeV-scale mass.
            assert!(id.state() <= 9, "state out of range: {nucid}");
            let (lo, hi) = (0.9 * f64::from(id.a()), 1.2 * f64::from(id.a()));
            assert!(*mass > lo && *mass < hi, "{} mass {mass}", id.to_name());
            assert!(*mass > 0.0);
        }
    }

    #[test]
    fn vendored_row_counts_match_tables() {
        let mass_rows = AME2020_TSV
            .lines()
            .filter(|l| !l.is_empty() && !l.starts_with('#'))
            .count();
        let abundance_rows = NATURAL_ABUNDANCE_TSV
            .lines()
            .filter(|l| !l.is_empty() && !l.starts_with('#'))
            .count();
        assert_eq!(mass_table().len(), mass_rows);
        assert_eq!(mass_rows, 3557 + 738);
        assert_eq!(abundance_table().len(), abundance_rows);
        assert_eq!(abundance_rows, 289);
    }

    #[test]
    fn isomer_masses_follow_ground_plus_excitation() {
        // Ba137_m1: ELIS 661659 eV on a 136.905827207 u ground state.
        let ba = NuclideId::from_name("Ba137").unwrap().nucid();
        let ba_m1 = NuclideId::from_name("Ba137_m1").unwrap().nucid();
        let expected = atomic_mass(ba).unwrap() + 0.661_659 / MEV_PER_U;
        assert!((atomic_mass(ba_m1).unwrap() - expected).abs() < 1e-9);
        // Second isomers resolve too (ENDF m2 tapes, e.g. Cu70_m2).
        let cu_m2 = NuclideId::from_name("Cu70_m2").unwrap().nucid();
        assert!(atomic_mass(cu_m2).unwrap() >= atomic_mass_by_name("Cu70").unwrap());
        // Pm137m has no ENDF/B-VIII.0 isomer tape: absent, no fill-in, so
        // the lookup falls back to the ground-state mass.
        assert_eq!(
            atomic_mass_by_name("Pm137_m1"),
            atomic_mass_by_name("Pm137")
        );
        // Stable Te123 keeps its ground mass; its isomer has its own row.
        assert!(atomic_mass_by_name("Te123_m1").unwrap() > atomic_mass_by_name("Te123").unwrap());
    }

    #[test]
    fn ame_masses_facade_delegates() {
        let provider = AmeMasses;
        assert_eq!(provider.atomic_mass(U235), atomic_mass(U235));
        assert_eq!(provider.atomic_mass_by_name("Fe56"), Some(55.934_935_537));
        assert_eq!(provider.natural_abundance(O16), Some(0.997_620_6));
        assert_eq!(provider.natural_abundance_by_name("Nope1"), None);
    }

    const U238: u32 = 922_380_000;
    const I135: u32 = 531_350_000;
    const CS137: u32 = 551_370_000;

    #[test]
    fn half_life_spot_values() {
        // U-238: 4.468e9 yr × 3.1557e7 s/yr ≈ 1.4100e17 s.
        let t_u238 = half_life(U238).unwrap();
        assert!((t_u238 - 1.409_99e17).abs() / t_u238 < 1e-9);
        // I-135: chain fixture value 2.3652e4 s.
        assert_eq!(half_life(I135), Some(23_652.0));
        // Cs-135: 7.25825e13 s (≈ 2.3 Myr).
        let t_cs135 = half_life_by_name("Cs135").unwrap();
        assert!((t_cs135 - 7.258_25e13).abs() / t_cs135 < 1e-12);
        // Cs-137: 30.08 yr ≈ 9.49e8 s, in the production Julian-year
        // convention (AD-11: 31,557,600 s/yr; the Gregorian divisor it
        // replaces differed by ~6e-4 yr against this 0.01 yr tolerance).
        let t_cs137 = half_life(CS137).unwrap();
        assert!((t_cs137 / 3.155_76e7 - 30.08).abs() < 0.01, "{t_cs137}");
        // Isomers carry their own rows.
        assert_eq!(half_life_by_name("Am242_m1"), Some(4_449_622_000.0));
    }

    #[test]
    fn stable_and_unknown_nuclides_have_no_half_life() {
        // Stable nuclides are absent from the vendored table.
        assert_eq!(half_life(O16), None);
        assert_eq!(half_life_by_name("Fe56"), None);
        assert_eq!(decay_constant(H1), None);
        assert_eq!(half_life(999_999_999), None);
        assert_eq!(half_life_by_name("Xx999"), None);
    }

    #[test]
    fn decay_constant_is_ln2_over_half_life() {
        for nucid in [
            U238,
            I135,
            CS137,
            NuclideId::from_name("Te132").unwrap().nucid(),
        ] {
            let t_half = half_life(nucid).unwrap();
            let lambda = decay_constant(nucid).unwrap();
            let rel = (lambda * t_half - std::f64::consts::LN_2).abs() / std::f64::consts::LN_2;
            assert!(rel < 1e-12, "nucid {nucid}: rel err {rel}");
        }
        let lam = decay_constant_by_name("I135").unwrap();
        assert!((lam - std::f64::consts::LN_2 / 23_652.0).abs() < 1e-18);
    }

    #[test]
    fn shorter_half_life_gives_larger_decay_constant() {
        // Te132 (3.2 d) vs I135 (6.57 h) vs Xe135 (9.14 h) chain ordering.
        let te = NuclideId::from_name("Te132").unwrap().nucid();
        let xe = NuclideId::from_name("Xe135").unwrap().nucid();
        let (t_te, t_i, t_xe) = (
            half_life(te).unwrap(),
            half_life(I135).unwrap(),
            half_life(xe).unwrap(),
        );
        assert!(t_te > t_xe && t_i < t_xe);
        assert!(decay_constant(te).unwrap() < decay_constant(xe).unwrap());
        assert!(decay_constant(xe).unwrap() < decay_constant(I135).unwrap());
    }

    #[test]
    fn half_lives_are_positive_and_finite() {
        for (nucid, t_half) in half_life_table() {
            assert!(*t_half > 0.0 && t_half.is_finite(), "{nucid}: {t_half}");
            let id = NuclideId::from_nucid(*nucid);
            assert!(id.a() >= id.z(), "{}", id.to_name());
        }
    }

    #[test]
    fn vendored_half_life_row_count_matches_table() {
        let rows = HALF_LIFE_TSV
            .lines()
            .filter(|l| !l.is_empty() && !l.starts_with('#'))
            .count();
        assert_eq!(half_life_table().len(), rows);
        assert_eq!(rows, 3561);
    }

    #[test]
    fn neutron_capture_q_value_anchors() {
        // H1(n,γ): textbook 2.224566 MeV.
        let q_h = q_value_neutron_capture(H1).unwrap();
        assert!((q_h - 2.224_566).abs() < 1e-3, "{q_h}");
        // U238(n,γ): ≈ 4.8 MeV.
        let q_u238 = q_value_neutron_capture(U238).unwrap();
        assert!((q_u238 - 4.806_382).abs() < 1e-3, "{q_u238}");
        assert!(q_u238 > 4.79 && q_u238 < 4.81);
        // O16(n,γ): 4.143 MeV.
        let q_o16 = q_value_neutron_capture(O16).unwrap();
        assert!((q_o16 - 4.143_080).abs() < 1e-3, "{q_o16}");
    }

    #[test]
    fn capture_q_value_matches_manual_formula() {
        let expected = (atomic_mass(U238).unwrap() + NEUTRON_MASS_U
            - atomic_mass(922_390_000).unwrap())
            * MEV_PER_U;
        let q = q_value_neutron_capture(U238).unwrap();
        assert!((q - expected).abs() < 1e-9);
        assert_eq!(q_value_neutron_capture_by_name("U238"), Some(q));
        assert_eq!(neutron_mass_u(), NEUTRON_MASS_U);
        assert_eq!(neutron_mass_u(), 1.008_664_915_95);
    }

    #[test]
    fn capture_q_value_missing_or_bad_targets_return_none() {
        // Product nuclide beyond the chart of nuclides (no A+1 mass row).
        let he10 = NuclideId::from_name("He10").unwrap().nucid();
        assert_eq!(atomic_mass(he10 + 10_000), None);
        assert_eq!(q_value_neutron_capture(he10), None);
        // Metastable states and non-nuclides are rejected outright.
        assert_eq!(q_value_neutron_capture(922_350_001), None);
        assert_eq!(q_value_neutron_capture(10_000), None);
        assert_eq!(q_value_neutron_capture_by_name("U235_m1"), None);
        assert_eq!(q_value_neutron_capture_by_name("Nope1"), None);
    }

    #[test]
    fn alpha_q_value_anchors() {
        // Literature: U238 α → 4.2698 MeV, Po210 → 5.4075 MeV,
        // Ra226 → 4.8706 MeV; atomic-mass formula reproduces all three.
        for (name, lit) in [
            ("U238", 4.269_858),
            ("Po210", 5.407_530),
            ("Ra226", 4.870_703),
        ] {
            let q = q_value_alpha_by_name(name).unwrap();
            assert!((q - lit).abs() < 1e-3, "{name}: {q} vs {lit}");
        }
        assert_eq!(q_value_alpha(U238), q_value_alpha_by_name("U238"));
    }

    #[test]
    fn alpha_q_value_rejects_light_and_metastable() {
        // Too light to alpha-decay within the table's Z/A domain.
        assert_eq!(q_value_alpha(H1), None);
        assert_eq!(q_value_alpha_by_name("He4"), None);
        assert_eq!(q_value_alpha(922_350_001), None);
        assert_eq!(q_value_alpha_by_name("Am242_m1"), None);
        // Endothermic "decay" still yields the (negative) Q-value:
        // O16 → C12 + He4 costs ≈ 7.162 MeV.
        let q_o16 = q_value_alpha_by_name("O16").unwrap();
        assert!((q_o16 + 7.162).abs() < 1e-3, "{q_o16}");
    }

    #[test]
    fn decay_data_facade_delegates() {
        let provider = DecayData;
        assert_eq!(provider.half_life(I135), Some(23_652.0));
        assert_eq!(provider.half_life_by_name("I135"), Some(23_652.0));
        assert_eq!(provider.decay_constant(U238), decay_constant(U238));
        assert_eq!(provider.decay_constant_by_name("Fe56"), None);
    }

    #[test]
    fn simple_xs_h1_anchors_within_ten_percent() {
        // Generated values: H free-atom thermal ~20.84 b, ~0.69 b at 14 MeV.
        let (thermal, fast) = simple_xs_by_name("H1").unwrap();
        assert!((thermal - 20.84).abs() / 20.84 < 0.05, "{thermal}");
        assert!((fast - 0.687).abs() / 0.687 < 0.05, "{fast}");
        assert_eq!(simple_xs(H1), Some((thermal, fast)));
    }

    #[test]
    fn simple_xs_absorber_and_actinide_bands() {
        // Strong absorbers dominate at thermal; actinide totals cluster ~6 b
        // at 14 MeV (ENDF/B-VII.1 MF3/MT1).
        let (b10_th, _) = simple_xs_by_name("B10").unwrap();
        assert!(b10_th > 3700.0 && b10_th < 3950.0, "{b10_th}");
        let (u235_th, u235_fast) = simple_xs_by_name("U235").unwrap();
        assert!(u235_th > 680.0 && u235_th < 710.0, "{u235_th}");
        assert!(u235_fast > 5.0 && u235_fast < 7.0, "{u235_fast}");
        let (pu239_th, _) = simple_xs_by_name("Pu239").unwrap();
        assert!(pu239_th > 1000.0 && pu239_th < 1050.0, "{pu239_th}");
        let (o16_th, _) = simple_xs_by_name("O16").unwrap();
        assert!(o16_th > 3.5 && o16_th < 4.0, "{o16_th}");
    }

    #[test]
    fn simple_xs_coverage_gaps_are_documented() {
        // Resonance absorbers without NIST rows, and isomers, are absent.
        for name in ["Cs137", "Co60", "I135", "Xe135", "Am242_m1"] {
            assert_eq!(simple_xs_by_name(name), None, "{name}");
        }
    }

    #[test]
    fn simple_xs_unknown_returns_none() {
        assert_eq!(simple_xs(999_999_999), None);
        assert_eq!(simple_xs_by_name("Og294"), None);
        assert_eq!(simple_xs_by_name("Xx999"), None);
    }

    #[test]
    fn scattering_length_nist_anchors() {
        // NIST NCNR bound lengths (Sears 1992): H1, D, O16.
        let (h_coh, h_inc) = scattering_length_by_name("H1").unwrap();
        assert!((h_coh - -3.7406).abs() < 1e-3, "{h_coh}");
        assert!((h_inc - 25.274).abs() < 1e-3, "{h_inc}");
        let (d_coh, d_inc) = scattering_length_by_name("H2").unwrap();
        assert!((d_coh - 6.671).abs() < 1e-2, "{d_coh}");
        assert!((d_inc - 4.04).abs() < 1e-2, "{d_inc}");
        let (o_coh, o_inc) = scattering_length_by_name("O16").unwrap();
        assert!((o_coh - 5.803).abs() < 1e-2, "{o_coh}");
        assert_eq!(o_inc, 0.0);
    }

    #[test]
    fn scattering_length_unknown_returns_none() {
        assert_eq!(scattering_length(999_999_999), None);
        assert_eq!(scattering_length_by_name("Og294"), None);
        assert_eq!(scattering_length_by_name("Xx999"), None);
    }

    #[test]
    fn decay_energy_anchors_within_tolerance() {
        // ENDF/B-VII.1 MF8/MT457 prompt means: Cs137 excludes the 662 keV
        // daughter gamma (Ba137_m1's row), Co60 ~2.60, H3 ~0.0057 MeV.
        let cs = decay_energy_mev_by_name("Cs137").unwrap();
        assert!((cs - 0.1794).abs() / 0.1794 < 0.05, "{cs}");
        let co = decay_energy_mev_by_name("Co60").unwrap();
        assert!((co - 2.6006).abs() / 2.6006 < 0.05, "{co}");
        let h3 = decay_energy_mev_by_name("H3").unwrap();
        assert!((h3 - 0.00569).abs() / 0.00569 < 0.05, "{h3}");
    }

    #[test]
    fn decay_energy_stable_and_unknown_return_none() {
        assert_eq!(decay_energy_mev(O16), None);
        assert_eq!(decay_energy_mev(FE56), None);
        assert_eq!(decay_energy_mev(999_999_999), None);
        assert_eq!(decay_energy_mev_by_name("Fe56"), None);
        assert_eq!(decay_energy_mev_by_name("Xx999"), None);
    }

    #[test]
    fn decay_energy_isomers_carry_own_rows() {
        // Ba137m prompt 661 keV gamma (×0.9 branching) vs Cs137 prompt only.
        let ba_m1 = decay_energy_mev_by_name("Ba137_m1").unwrap();
        assert!((ba_m1 - 0.6614).abs() / 0.6614 < 0.05, "{ba_m1}");
        assert!(decay_energy_mev_by_name("Ba137").is_none());
        // Second isomers resolve too (ENDF m2 tapes, e.g. Am242_m2).
        let am_m2 = NuclideId::from_name("Am242_m2").unwrap().nucid();
        assert_eq!(
            decay_energy_mev_by_name("Am242_m2"),
            decay_energy_mev(am_m2)
        );
        assert!(decay_energy_mev(am_m2).is_some());
    }

    #[test]
    fn vendored_generated_row_counts_match_tables() {
        for (tsv, table_len, expected) in [
            (SIMPLE_XS_TSV, simple_xs_table().len(), 241),
            (SCATTERING_LENGTHS_TSV, scattering_length_table().len(), 267),
            (DECAY_ENERGY_TSV, decay_energy_table().len(), 3557),
        ] {
            let rows = tsv
                .lines()
                .filter(|l| !l.is_empty() && !l.starts_with('#'))
                .count();
            assert_eq!(table_len, rows);
            assert_eq!(rows, expected);
        }
    }

    #[test]
    fn generated_by_name_agrees_with_nucid_lookup() {
        for name in ["H1", "B10", "O16", "Fe56", "U235", "Pu239"] {
            let nucid = NuclideId::from_name(name).unwrap().nucid();
            assert_eq!(simple_xs_by_name(name), simple_xs(nucid), "{name}");
            assert_eq!(
                scattering_length_by_name(name),
                scattering_length(nucid),
                "{name}"
            );
        }
        for name in ["H3", "Co60", "Cs137"] {
            let nucid = NuclideId::from_name(name).unwrap().nucid();
            assert_eq!(
                decay_energy_mev_by_name(name),
                decay_energy_mev(nucid),
                "{name}"
            );
        }
    }

    #[test]
    fn decay_data_facade_delegates_decay_energy() {
        let provider = DecayData;
        assert_eq!(provider.decay_energy_mev(CS137), decay_energy_mev(CS137));
        let co60 = provider.decay_energy_mev_by_name("Co60").unwrap();
        assert!((co60 - 2.6006).abs() / 2.6006 < 0.05, "{co60}");
        assert_eq!(provider.decay_energy_mev_by_name("Fe56"), None);
    }

    #[test]
    fn decay_branch_k40_two_branches_sum_to_one() {
        // ENDF/B-VIII.0 NDK values, verbatim: beta- 0.8914 to Ca40 and
        // EC/beta+ 0.1086 to Ar40.
        let k40 = NuclideId::from_name("K40").unwrap().nucid();
        let branches = decay_branches(k40).unwrap();
        assert_eq!(branches.len(), 2);
        let ca40 = NuclideId::from_name("Ca40").unwrap().nucid();
        let ar40 = NuclideId::from_name("Ar40").unwrap().nucid();
        assert!(branches.contains(&DecayBranch {
            progeny: ca40,
            branching_fraction: 0.8914,
            mode: DecayBranchMode::BetaMinus,
        }));
        assert!(branches.contains(&DecayBranch {
            progeny: ar40,
            branching_fraction: 0.1086,
            mode: DecayBranchMode::EcBetaPlus,
        }));
        let total: f64 = branches.iter().map(|b| b.branching_fraction).sum();
        assert!((total - 1.0).abs() < 1e-9, "{total}");
        assert_eq!(branching_fraction(k40, ca40), Some(0.8914));
        assert_eq!(branching_fraction_by_name("K40", "Ar40"), Some(0.1086));
        assert_eq!(branching_fraction_by_name("K40", "K40"), None);
    }

    #[test]
    fn decay_branch_spot_modes() {
        use DecayBranchMode as M;
        // Single alpha branch with unit BR.
        let es254 = decay_branches_by_name("Es254").unwrap();
        assert_eq!(es254.len(), 1);
        assert_eq!(es254[0].mode, M::Alpha);
        assert_eq!(es254[0].branching_fraction, 1.0);
        assert_eq!(NuclideId::from_nucid(es254[0].progeny).to_name(), "Bk250");
        // Isomeric transition keeps Z/A.
        let ba = decay_branches_by_name("Ba137_m1").unwrap();
        assert_eq!(ba.len(), 1);
        assert_eq!(ba[0].mode, M::It);
        assert_eq!(NuclideId::from_nucid(ba[0].progeny).to_name(), "Ba137");
        // Delayed-neutron branch collapses to the beta- initial event with
        // the emitted neutron subtracted from the progeny.
        let he8 = decay_branches_by_name("He8").unwrap();
        assert_eq!(he8.len(), 2);
        let names: Vec<String> = he8
            .iter()
            .map(|b| NuclideId::from_nucid(b.progeny).to_name())
            .collect();
        assert!(names.contains(&"Li8".to_string()), "{names:?}");
        assert!(names.contains(&"Li7".to_string()), "{names:?}");
        assert!(he8.iter().all(|b| b.mode == M::BetaMinus));
        // Es254m mixes alpha, beta-, EC, and IT; the SF branch is dropped.
        let es_m1 = decay_branches_by_name("Es254_m1").unwrap();
        assert_eq!(es_m1.len(), 4);
        assert!(es_m1.iter().all(|b| b.mode != M::Sf));
        assert_eq!(branching_fraction_by_name("Es254_m1", "Fm254"), Some(0.98));
    }

    #[test]
    fn decay_branch_stable_and_unknown_have_no_rows() {
        assert_eq!(decay_branches_by_name("Fe56"), None);
        assert_eq!(decay_branches_by_name("O16"), None);
        // Zero-half-life evaluation dummies stay absent (stable-absent rule).
        assert_eq!(decay_branches_by_name("Te123"), None);
        assert_eq!(decay_branches_by_name("Ca46"), None);
        assert_eq!(branching_fraction_by_name("Fe56", "Fe56"), None);
        assert_eq!(decay_branches_by_name("Xx999"), None);
    }

    #[test]
    fn decay_branch_parents_agree_with_half_life_table() {
        // Every branch parent carries a half-life row (same VIII.0 basis);
        // progeny may be stable (absent) — only names must parse.
        for (parent, branches) in decay_branch_table() {
            assert!(
                half_life(*parent).is_some(),
                "parent without half-life: {parent}"
            );
            assert!(!branches.is_empty());
            for b in branches {
                assert!(
                    (0.0..=1.0).contains(&b.branching_fraction),
                    "BF range: {}",
                    b.branching_fraction
                );
                let id = NuclideId::from_nucid(b.progeny);
                assert!(id.z() >= 1 && id.a() >= id.z(), "{}", id.to_name());
            }
        }
    }

    #[test]
    fn decay_branch_row_count_matches_table() {
        let rows = DECAY_BRANCHES_TSV
            .lines()
            .filter(|l| !l.is_empty() && !l.starts_with('#'))
            .count();
        let table_rows: usize = decay_branch_table().values().map(Vec::len).sum();
        assert_eq!(table_rows, rows);
        assert_eq!(rows, 5068);
        assert_eq!(decay_branch_table().len(), 3541);
    }

    #[test]
    fn decay_branch_mode_parsing() {
        use DecayBranchMode as M;
        assert_eq!(M::parse("beta-"), Some(M::BetaMinus));
        assert_eq!(M::parse("ec/beta+"), Some(M::EcBetaPlus));
        assert_eq!(M::parse("EC/BETA+"), Some(M::EcBetaPlus));
        assert_eq!(M::parse("alpha"), Some(M::Alpha));
        assert_eq!(M::parse("IT"), Some(M::It));
        assert_eq!(M::parse("it"), Some(M::It));
        assert_eq!(M::parse("sf"), Some(M::Sf));
        assert_eq!(M::parse("n"), Some(M::Neutron));
        assert_eq!(M::parse("p"), Some(M::Proton));
        assert_eq!(M::parse("nope"), None);
        assert_eq!(M::BetaMinus.as_str(), "beta-");
        assert_eq!(M::EcBetaPlus.as_str(), "ec/beta+");
        assert_eq!(M::It.as_str(), "IT");
    }

    #[test]
    fn decay_data_facade_delegates_branches() {
        let provider = DecayData;
        let k40 = NuclideId::from_name("K40").unwrap().nucid();
        assert_eq!(provider.decay_branches(k40), decay_branches(k40));
        assert_eq!(
            provider.decay_branches_by_name("Es254"),
            decay_branches_by_name("Es254")
        );
        assert_eq!(
            provider.branching_fraction_by_name("K40", "Ca40"),
            Some(0.8914)
        );
        assert_eq!(provider.decay_branches_by_name("Fe56"), None);
    }

    #[test]
    fn dose_factor_row_count_matches_table() {
        let rows = DOSE_FACTORS_TSV
            .lines()
            .filter(|l| !l.is_empty() && !l.starts_with('#'))
            .count();
        assert_eq!(dose_table().len(), rows);
        assert_eq!(rows, 1116);
    }

    #[test]
    fn dose_factor_spot_values() {
        use DosePathway as P;
        use DoseSource as S;
        // Raw table factors from the PyNE CSVs (folded +D parents).
        assert_eq!(
            dose_factor_by_name("Co60", P::Ingest, S::Epa),
            Some(2.69e-05)
        );
        assert_eq!(
            dose_factor_by_name("Cs137", P::Inhale, S::Epa),
            Some(3.19e-05)
        );
        assert_eq!(dose_factor_by_name("H3", P::Air, S::Epa), Some(4.41e-012));
        assert_eq!(dose_factor_by_name("K40", P::Soil, S::Epa), Some(4.33e02));
    }

    #[test]
    fn dose_factor_missing_air_is_minus_one_sentinel() {
        use DosePathway as P;
        use DoseSource as S;
        // Air is EPA-only: GENII/DOE air rows are -1 sentinels (PyNE).
        let h3 = NuclideId::from_name("H3").unwrap().nucid();
        assert_eq!(dose_factor(h3, P::Air, S::Genii), Some(-1.0));
        assert_eq!(dose_factor(h3, P::Air, S::Doe), Some(-1.0));
        // Unknown nuclides have no row at all.
        assert_eq!(dose_factor(999_999_999, P::Ingest, S::Epa), None);
        assert_eq!(dose_factor_by_name("Fe56", P::Ingest, S::Epa), None);
        assert_eq!(dose_factor_by_name("Xx999", P::Ingest, S::Epa), None);
    }

    #[test]
    fn dose_aux_columns() {
        use DoseSource as S;
        // f1 lives on ingest rows; lung model on inhale rows.
        assert_eq!(dose_f1_by_name("Co60", S::Epa), Some(0.3));
        assert_eq!(dose_f1_by_name("H3", S::Epa), Some(1.0));
        assert_eq!(dose_lung_model_by_name("Co60", S::Epa), Some('Y'));
        assert_eq!(dose_lung_model_by_name("H3", S::Epa), Some('V'));
        assert_eq!(dose_lung_model_by_name("C14", S::Epa), Some('O'));
        // Air/soil rows carry neither aux column.
        let h3 = NuclideId::from_name("H3").unwrap().nucid();
        assert_eq!(dose_entry(h3, DosePathway::Air, S::Epa).unwrap().f1, None);
        assert_eq!(
            dose_entry(h3, DosePathway::Air, S::Epa).unwrap().lung_model,
            None
        );
    }

    #[test]
    fn dose_pathway_source_parsing() {
        assert_eq!(DosePathway::parse("air"), Some(DosePathway::Air));
        assert_eq!(DosePathway::parse("ext_air"), Some(DosePathway::Air));
        assert_eq!(DosePathway::parse("ext_soil"), Some(DosePathway::Soil));
        assert_eq!(DosePathway::parse("INGEST"), Some(DosePathway::Ingest));
        assert_eq!(DosePathway::parse("inhale"), Some(DosePathway::Inhale));
        assert_eq!(DosePathway::parse("nope"), None);
        assert_eq!(DoseSource::parse("epa"), Some(DoseSource::Epa));
        assert_eq!(DoseSource::parse("DOE"), Some(DoseSource::Doe));
        assert_eq!(DoseSource::parse("genii"), Some(DoseSource::Genii));
        assert_eq!(DoseSource::from_int(0), Some(DoseSource::Epa));
        assert_eq!(DoseSource::from_int(1), Some(DoseSource::Doe));
        assert_eq!(DoseSource::from_int(2), Some(DoseSource::Genii));
        assert_eq!(DoseSource::from_int(3), None);
        assert_eq!(DoseSource::Epa.to_int(), 0);
    }

    #[test]
    fn dose_data_facade_delegates() {
        use DosePathway as P;
        use DoseSource as S;
        let provider = DoseData;
        let co60 = NuclideId::from_name("Co60").unwrap().nucid();
        assert_eq!(
            provider.dose_factor(co60, P::Ingest, S::Epa),
            dose_factor(co60, P::Ingest, S::Epa)
        );
        assert_eq!(
            provider.dose_factor_by_name("K40", P::Soil, S::Epa),
            Some(4.33e02)
        );
        assert_eq!(provider.dose_f1_by_name("H3", S::Epa), Some(1.0));
        assert_eq!(provider.dose_lung_model_by_name("H3", S::Epa), Some('V'));
        assert_eq!(
            provider.dose_factor_by_name("Fe56", P::Ingest, S::Epa),
            None
        );
    }
}