oxiphysics-materials 0.1.0

// Copyright 2026 COOLJAPAN OU (Team KitaSan)
// SPDX-License-Identifier: Apache-2.0

//! Advanced nuclear reactor materials science.
//!
//! This module covers:
//!
//! - **Zircaloy cladding** — oxidation kinetics (cubic/linear breakaway),
//!   hydrogen pickup fraction, hydride embrittlement
//! - **UO₂ fuel pellet** — thermal conductivity vs burnup (Lucuta model),
//!   fission gas release (Booth diffusion), pellet cracking and swelling
//! - **Radiation void swelling** — Garner model for austenitic steels
//! - **Irradiation creep** — in-reactor creep for fuel cladding
//! - **Radiation-induced segregation (RIS)** — inverse Kirkendall mechanism
//! - **Radiation embrittlement** — ductile-to-brittle transition temperature
//!   shift (Charpy / USE models for RPV steel)
//! - **Helium embrittlement** — grain boundary helium bubble model
//! - **Reactor pressure vessel (RPV) steel** — Charpy USE correlation,
//!   master curve approach (T₀ reference temperature)
//! - **Graphite moderator** — dimensional change, thermal conductivity
//!   degradation under irradiation
//! - **Sodium fast reactor structural materials** — 316 SS and HT-9 ferritic-
//!   martensitic steel properties, creep-fatigue interaction

#![allow(dead_code)]
#![allow(clippy::too_many_arguments)]

use std::f64::consts::PI;

// ---------------------------------------------------------------------------
// Section 1 — Zircaloy cladding
// ---------------------------------------------------------------------------

/// Zircaloy alloy grade.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum ZircaloyGrade {
    /// Zircaloy-2 (0.1% Sn, used in BWR).
    Zircaloy2,
    /// Zircaloy-4 (1.45% Sn, used in PWR).
    Zircaloy4,
    /// ZIRLO — enhanced corrosion resistance.
    Zirlo,
    /// M5 — Nb-based alloy.
    M5,
}

/// Zircaloy cladding tube.
#[derive(Debug, Clone)]
pub struct ZircaloyClad {
    /// Alloy grade.
    pub grade: ZircaloyGrade,
    /// Outer diameter (m).
    pub outer_diameter: f64,
    /// Wall thickness (m).
    pub wall_thickness: f64,
    /// Current oxide layer thickness (m).
    pub oxide_thickness: f64,
    /// Cumulative fast neutron fluence (n/m², E > 1 MeV).
    pub fluence: f64,
    /// Operating temperature (K).
    pub temperature: f64,
}

impl ZircaloyClad {
    /// Construct a new Zircaloy cladding tube.
    pub fn new(
        grade: ZircaloyGrade,
        outer_diameter: f64,
        wall_thickness: f64,
        temperature: f64,
    ) -> Self {
        Self {
            grade,
            outer_diameter,
            wall_thickness,
            oxide_thickness: 0.0,
            fluence: 0.0,
            temperature,
        }
    }

    /// Corrosion rate constant A (kg²/m⁴/s) for the pre-breakaway cubic law.
    ///
    /// Empirical fit to Urbanic & Heidrick (1978) data.
    fn cubic_rate_constant(&self) -> f64 {
        let ea_j = match self.grade {
            ZircaloyGrade::Zircaloy2 => 1.42e4_f64,
            ZircaloyGrade::Zircaloy4 => 1.35e4_f64,
            ZircaloyGrade::Zirlo => 1.20e4_f64,
            ZircaloyGrade::M5 => 1.10e4_f64,
        };
        // A = A0 * exp(-Ea/(R*T))
        let a0 = 2.88e-3_f64; // kg²/(m⁴·s)
        let r = 8.314_f64;
        a0 * (-(ea_j) / (r * self.temperature)).exp()
    }

    /// Oxide weight gain (kg/m²) after time `dt` seconds using cubic parabolic law.
    ///
    /// Δw³ = A · t  (pre-transition regime, w < transition weight gain w_t)
    pub fn oxide_weight_gain_cubic(&self, time_s: f64) -> f64 {
        let a = self.cubic_rate_constant();
        (a * time_s).cbrt()
    }

    /// Post-transition linear oxidation rate constant (kg/m²/s).
    ///
    /// After breakaway, a linear rate governs: Δw = B · (t − t_t).
    pub fn linear_rate_constant(&self) -> f64 {
        // Empirical B values
        let ea_j = 1.10e4_f64;
        let b0 = 6.0e-6_f64;
        let r = 8.314_f64;
        b0 * (-(ea_j) / (r * self.temperature)).exp()
    }

    /// Breakaway transition time (s) — approximate.
    ///
    /// Returns the time at which oxide weight gain transitions from cubic to linear.
    pub fn breakaway_time(&self) -> f64 {
        // Transition weight gain ~ 30 mg/dm² = 0.30 kg/m² for Zircaloy-4
        let wt = match self.grade {
            ZircaloyGrade::Zircaloy2 => 0.30_f64,
            ZircaloyGrade::Zircaloy4 => 0.30_f64,
            ZircaloyGrade::Zirlo => 0.40_f64,
            ZircaloyGrade::M5 => 0.50_f64,
        };
        wt * wt * wt / self.cubic_rate_constant()
    }

    /// Hydrogen pickup fraction (dimensionless, 0–1).
    ///
    /// Fraction of hydrogen generated by oxidation that is absorbed by the Zr matrix.
    /// Based on empirical correlation with temperature and fluence.
    pub fn hydrogen_pickup_fraction(&self) -> f64 {
        // f = f0 * exp(Ea/RT) * (1 + k_phi * fluence)
        let f0 = match self.grade {
            ZircaloyGrade::Zircaloy2 => 0.15_f64,
            ZircaloyGrade::Zircaloy4 => 0.10_f64,
            ZircaloyGrade::Zirlo => 0.08_f64,
            ZircaloyGrade::M5 => 0.07_f64,
        };
        let ea = 2500.0_f64; // J/mol
        let r = 8.314_f64;
        let fluence_factor = 1.0 + 2.0e-26 * self.fluence;
        (f0 * (ea / (r * self.temperature)).exp() * fluence_factor).min(1.0)
    }

    /// Hydrogen concentration in cladding (wppm).
    ///
    /// Computed from cumulative oxide weight gain and pickup fraction.
    /// `oxide_wg` is oxide weight gain (kg/m²).
    pub fn hydrogen_concentration_wppm(&self, oxide_wg: f64) -> f64 {
        // Every kg/m² of ZrO2 generates ~1.11e4 wppm·m² if all absorbed
        // Scaling for wall thickness
        let rho_zr = 6520.0_f64; // kg/m³
        let wt = self.wall_thickness;
        let h_gen = oxide_wg * 1.11e4 / (rho_zr * wt);
        h_gen * self.hydrogen_pickup_fraction()
    }

    /// Yield strength (Pa) of irradiated Zircaloy.
    ///
    /// Uses the Northwood & Gilbert correlation including irradiation hardening.
    pub fn yield_strength(&self) -> f64 {
        // Base yield strength (Pa) at temperature
        let sigma_0 = match self.grade {
            ZircaloyGrade::Zircaloy2 => 450e6_f64,
            ZircaloyGrade::Zircaloy4 => 480e6_f64,
            ZircaloyGrade::Zirlo => 500e6_f64,
            ZircaloyGrade::M5 => 510e6_f64,
        };
        // Temperature softening
        let t0 = 293.0_f64;
        let softening = 1.0 - 5.0e-4 * (self.temperature - t0);
        // Irradiation hardening: Δσ = k * sqrt(Φ)
        let k_irr = 1.2e-12_f64;
        let irr_hardening = k_irr * self.fluence.sqrt();
        (sigma_0 * softening + irr_hardening).max(0.0)
    }

    /// Irradiation creep rate (s⁻¹) under hoop stress `sigma` (Pa).
    ///
    /// In-reactor creep: ε̇ = B_0 · σ · Φ̇  (irradiation creep component).
    pub fn irradiation_creep_rate(&self, sigma_pa: f64, flux_n_m2_s: f64) -> f64 {
        // B_0 ~ 3e-27 m²/(n·Pa)  for Zircaloy
        let b0 = 3.0e-27_f64;
        b0 * sigma_pa * flux_n_m2_s
    }

    /// Ductility reduction factor due to hydrides (at test temperature).
    ///
    /// Returns a multiplier (1.0 = no degradation, → 0 at very high H).
    pub fn hydride_ductility_factor(&self, h_wppm: f64) -> f64 {
        // Simple exponential decay — EPRI correlation
        (1.0 - h_wppm / 1000.0).max(0.05)
    }
}

// ---------------------------------------------------------------------------
// Section 2 — UO₂ fuel pellet
// ---------------------------------------------------------------------------

/// UO₂ fuel pellet.
#[derive(Debug, Clone)]
pub struct Uo2Pellet {
    /// Pellet diameter (m).
    pub diameter: f64,
    /// Pellet height (m).
    pub height: f64,
    /// Burnup (MWd/kgU).
    pub burnup: f64,
    /// Centre-line temperature (K).
    pub temperature_centre: f64,
    /// Pellet-cladding gap width (m).
    pub gap_width: f64,
    /// Enrichment (weight fraction U-235, 0–1).
    pub enrichment: f64,
}

impl Uo2Pellet {
    /// Construct a standard UO₂ pellet.
    pub fn new(diameter: f64, height: f64, enrichment: f64) -> Self {
        Self {
            diameter,
            height,
            burnup: 0.0,
            temperature_centre: 900.0,
            gap_width: 1.0e-4,
            enrichment,
        }
    }

    /// Thermal conductivity of UO₂ (W/m/K) as a function of burnup.
    ///
    /// Uses the Lucuta et al. (1996) model:
    /// k = k_0(T) · F_D · F_P · F_R · F_Por
    pub fn thermal_conductivity(&self, temperature_k: f64) -> f64 {
        // Base conductivity (unirradiated, theoretical density)
        let t = temperature_k;
        let k0 = 1.0 / (0.0452 + 2.46e-4 * t) + 88.0e9 / (t * t) * (-18531.7 / t).exp();

        // Dissolved fission product factor F_D (Lucuta)
        let bu = self.burnup.min(100.0);
        let f_d = 1.0 / (1.0 + 0.019 * bu / (3.0 - 0.019 * bu));

        // Precipitate factor F_P
        let f_p = 1.0 - 0.2 / (1.0 + (t - 900.0).powi(2) / 40000.0) * (bu / 100.0);

        // Radiation damage factor F_R
        let f_r = 1.0 - 0.18 / (1.0 + (t - 900.0).powi(2) / 40000.0) * (bu / 100.0);

        // Porosity factor (assume 95% TD → ~5% porosity)
        let porosity = 0.05_f64;
        let f_por = (1.0 - porosity).powf(2.5);

        k0 * f_d * f_p * f_r * f_por
    }

    /// Fission gas release fraction (dimensionless, 0–1).
    ///
    /// Booth single-sphere diffusion model.  `temperature_k` is the representative
    /// temperature.  `time_s` is irradiation time.
    pub fn fission_gas_release_booth(&self, temperature_k: f64, time_s: f64) -> f64 {
        // Diffusion coefficient of Xe in UO₂
        let d0 = 5.0e-8_f64; // m²/s (pre-exponential)
        let ea = 40200.0_f64; // J/mol
        let r = 8.314_f64;
        let d = d0 * (-(ea) / (r * temperature_k)).exp();
        // Booth's sphere radius (equivalent)
        let radius = self.diameter / 2.0;
        let tau = d * time_s / (radius * radius);
        // Fractional release (Booth 1957)
        if tau < 0.02 {
            6.0 * (tau / PI).sqrt() - 3.0 * tau
        } else {
            1.0 - (6.0 / (PI * PI))
                * (1..=20)
                    .map(|n| (-(n as f64).powi(2) * PI * PI * tau).exp() / (n as f64).powi(2))
                    .sum::<f64>()
        }
        .clamp(0.0, 1.0)
    }

    /// Radial cracking factor — ratio of actual effective conductivity to uncracked.
    ///
    /// Approximate power-law fit to Oguma (1983) data.
    pub fn cracking_conductivity_factor(&self) -> f64 {
        // Decreases with burnup due to pellet cracking/relocation
        let bu = self.burnup;
        1.0 / (1.0 + 0.005 * bu)
    }

    /// Volumetric swelling (ΔV/V₀, dimensionless) due to fission products.
    ///
    /// Empirical correlation: ~0.98 vol% per 10 GWd/tU burnup.
    pub fn fission_product_swelling(&self) -> f64 {
        0.98e-2 * self.burnup / 10.0
    }

    /// Gap conductance (W/m²/K).
    ///
    /// Simple model: k_gap / gap_width + radiation term.
    pub fn gap_conductance(&self, fill_gas_conductivity: f64) -> f64 {
        let k_gas = fill_gas_conductivity.max(0.01);
        let gap = self.gap_width.max(1e-7);
        // Conduction through gas + radiation (Stefan-Boltzmann, simplified)
        let t_clad = self.temperature_centre - 400.0; // approximate
        let sigma = 5.67e-8_f64;
        let emissivity = 0.85_f64;
        let radiation = 4.0 * sigma * emissivity * t_clad.powi(3);
        k_gas / gap + radiation
    }

    /// Linear heat generation rate (W/m) for a given specific power (W/kgU).
    pub fn linear_heat_rate(&self, specific_power_w_per_kgu: f64) -> f64 {
        let rho_uo2 = 10_400.0_f64; // kg/m³ at 95%TD
        let area = PI * (self.diameter / 2.0).powi(2);
        // Mass of uranium per unit length
        let u_mass_fraction = 238.03 / (238.03 + 2.0 * 16.0); // UO₂ molar mass ratio
        let m_u_per_m = rho_uo2 * area * u_mass_fraction;
        specific_power_w_per_kgu * m_u_per_m / 1000.0
    }
}

// ---------------------------------------------------------------------------
// Section 3 — Radiation void swelling (Garner model)
// ---------------------------------------------------------------------------

/// Void swelling parameters for austenitic stainless steel (Garner model).
#[derive(Debug, Clone, Copy)]
pub struct VoidSwellingParams {
    /// Transient swelling incubation dose (dpa).
    pub incubation_dose_dpa: f64,
    /// Steady-state swelling rate (%/dpa) after incubation.
    pub swelling_rate_pct_per_dpa: f64,
    /// Temperature at peak swelling (K).
    pub peak_temperature_k: f64,
    /// Width of temperature bell curve (K).
    pub temperature_width_k: f64,
}

impl Default for VoidSwellingParams {
    fn default() -> Self {
        // Typical 316 SS (Garner 1994)
        Self {
            incubation_dose_dpa: 10.0,
            swelling_rate_pct_per_dpa: 1.0,
            peak_temperature_k: 773.0,
            temperature_width_k: 80.0,
        }
    }
}

/// Compute void swelling fraction (ΔV/V) using the Garner ramp-rate model.
///
/// S(Φ) = 0  for Φ < Φ_inc
/// S(Φ) = ṡ · (Φ − Φ_inc) · f_T(T)  for Φ ≥ Φ_inc
pub fn garner_void_swelling(dose_dpa: f64, temperature_k: f64, params: &VoidSwellingParams) -> f64 {
    if dose_dpa <= params.incubation_dose_dpa {
        return 0.0;
    }
    let excess_dpa = dose_dpa - params.incubation_dose_dpa;
    // Temperature factor: Gaussian bell
    let dt = temperature_k - params.peak_temperature_k;
    let f_t = (-(dt * dt) / (2.0 * params.temperature_width_k * params.temperature_width_k)).exp();
    let swelling_pct = params.swelling_rate_pct_per_dpa * excess_dpa * f_t;
    swelling_pct / 100.0
}

/// Void number density (m⁻³) — approximate from swelling and mean void radius.
pub fn void_number_density(swelling: f64, mean_void_radius_m: f64) -> f64 {
    if mean_void_radius_m < 1e-15 {
        return 0.0;
    }
    swelling / ((4.0 / 3.0) * PI * mean_void_radius_m.powi(3))
}

// ---------------------------------------------------------------------------
// Section 4 — Irradiation creep
// ---------------------------------------------------------------------------

/// Irradiation creep rate (s⁻¹) for austenitic stainless steel.
///
/// Total creep rate = thermal creep + irradiation-enhanced creep.
///
/// Uses the simplified model: ε̇_irr = B_irr · σ · Φ̇
/// and ε̇_th = A · σⁿ · exp(−Q / RT).
pub fn irradiation_creep_rate(
    stress_pa: f64,
    temperature_k: f64,
    flux_dpa_per_s: f64,
    b_irr: f64,
    a_thermal: f64,
    n_creep: f64,
    q_creep: f64,
) -> f64 {
    let r = 8.314_f64;
    let eps_irr = b_irr * stress_pa * flux_dpa_per_s;
    let eps_th = a_thermal * stress_pa.powf(n_creep) * (-(q_creep) / (r * temperature_k)).exp();
    eps_irr + eps_th
}

/// Stress relaxation factor due to irradiation creep after dose `dose_dpa`.
///
/// σ(Φ) = σ₀ · exp(−B_irr · Φ)
pub fn irradiation_stress_relaxation(sigma_0: f64, b_irr: f64, dose_dpa: f64) -> f64 {
    sigma_0 * (-(b_irr * dose_dpa)).exp()
}

// ---------------------------------------------------------------------------
// Section 5 — Radiation-induced segregation (RIS)
// ---------------------------------------------------------------------------

/// Radiation-induced segregation model (inverse Kirkendall / thermodynamic).
///
/// Computes the segregated solute concentration at a grain boundary sink.
#[derive(Debug, Clone, Copy)]
pub struct RisParams {
    /// Bulk solute concentration (atomic fraction).
    pub bulk_concentration: f64,
    /// Solute-interstitial binding energy (eV).
    pub binding_energy_ev: f64,
    /// Temperature (K).
    pub temperature_k: f64,
    /// Dose rate (dpa/s).
    pub dose_rate_dpa_s: f64,
    /// Recombination coefficient (dimensionless).
    pub recombination_coeff: f64,
}

impl Default for RisParams {
    fn default() -> Self {
        Self {
            bulk_concentration: 0.18, // 18 at% Cr in 316 SS
            binding_energy_ev: 0.15,
            temperature_k: 573.0,
            dose_rate_dpa_s: 1e-8,
            recombination_coeff: 1e17,
        }
    }
}

/// Steady-state solute depletion at a grain boundary (Cr depletion in stainless steel).
///
/// Returns the grain boundary concentration (atomic fraction) using the
/// simplified Perks (1986) two-sink model.
pub fn ris_grain_boundary_concentration(params: &RisParams) -> f64 {
    let kb = 8.617e-5_f64; // eV/K
    let kbt = kb * params.temperature_k;
    // Partition factor for solute flux: positive → depletion
    let binding_factor = (params.binding_energy_ev / kbt).exp();
    // Steady-state supersaturation of defects
    let defect_ratio = (params.dose_rate_dpa_s / params.recombination_coeff).sqrt();
    // Concentration change at GB
    let delta_c = params.bulk_concentration * (1.0 - binding_factor * defect_ratio);
    (params.bulk_concentration + delta_c).clamp(0.0, 1.0)
}

// ---------------------------------------------------------------------------
// Section 6 — Radiation embrittlement — RPV steel
// ---------------------------------------------------------------------------

/// Reactor pressure vessel steel grade.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum RpvSteelGrade {
    /// A508 Class 3 (forging — typical US PWR vessel).
    A508Class3,
    /// A533 Grade B (plate — typical US BWR vessel).
    A533GradeB,
    /// 15Kh2MFA (Russian VVER vessel steel).
    Vver15Kh2Mfa,
}

/// Charpy impact properties of RPV steel.
#[derive(Debug, Clone)]
pub struct RpvCharpy {
    /// Steel grade.
    pub grade: RpvSteelGrade,
    /// Unirradiated upper-shelf energy (J).
    pub use_unirradiated: f64,
    /// Unirradiated reference temperature (°C).
    pub rt_ndt_0: f64,
    /// Fast neutron fluence at vessel wall (n/m², E > 1 MeV).
    pub fluence: f64,
    /// Copper content (weight%).
    pub copper_wt: f64,
    /// Nickel content (weight%).
    pub nickel_wt: f64,
    /// Phosphorus content (weight%).
    pub phosphorus_wt: f64,
}

impl RpvCharpy {
    /// Construct a new RPV Charpy record.
    pub fn new(
        grade: RpvSteelGrade,
        use_unirradiated: f64,
        rt_ndt_0: f64,
        fluence: f64,
        copper_wt: f64,
        nickel_wt: f64,
        phosphorus_wt: f64,
    ) -> Self {
        Self {
            grade,
            use_unirradiated,
            rt_ndt_0,
            fluence,
            copper_wt,
            nickel_wt,
            phosphorus_wt,
        }
    }

    /// Charpy upper-shelf energy drop (ΔUE, J) due to irradiation.
    ///
    /// Regulatory Guide 1.99 Rev.2 correlation (NRC, 1988):
    /// ΔUE/USE₀ = A · Φ^(0.28 − 0.10 log Φ)
    pub fn use_drop(&self) -> f64 {
        if self.fluence < 1e18 {
            return 0.0;
        }
        let phi = self.fluence / 1e19_f64; // normalise to 10¹⁹ n/cm²
        // A factor depends on copper and nickel
        let a =
            0.11 * (1.0 + 1.58 * self.copper_wt.max(0.0)) * (1.0 + 3.77 * self.nickel_wt.max(0.0));
        let exponent = 0.28 - 0.10 * phi.log10();
        let fraction = a * phi.powf(exponent);
        fraction * self.use_unirradiated
    }

    /// Irradiated upper-shelf energy (J).
    pub fn use_irradiated(&self) -> f64 {
        (self.use_unirradiated - self.use_drop()).max(0.0)
    }

    /// RTNDT shift (°C) — Reg Guide 1.99 Rev.2 CF method.
    ///
    /// ΔRTNDT = CF · Φ^(0.28 − 0.10 log Φ)
    pub fn rtndt_shift(&self) -> f64 {
        if self.fluence < 1e18 {
            return 0.0;
        }
        let phi = self.fluence / 1e19_f64;
        let cf = chemistry_factor(self.copper_wt, self.nickel_wt);
        let exponent = 0.28 - 0.10 * phi.log10();
        cf * phi.powf(exponent)
    }

    /// Irradiated reference temperature RTNDT (°C).
    pub fn rtndt_irradiated(&self) -> f64 {
        self.rt_ndt_0 + self.rtndt_shift()
    }

    /// Fracture toughness K_Ic (MPa√m) from ASME Code Case N-629 master curve.
    ///
    /// K_Ic(T) = 30 + 70 · exp(0.019 · (T − T₀))  where T₀ = RTNDT − 33 °C.
    pub fn fracture_toughness_kic(&self, test_temperature_c: f64) -> f64 {
        let t0 = self.rtndt_irradiated() - 33.0;
        30.0 + 70.0 * (0.019 * (test_temperature_c - t0)).exp()
    }
}

/// Chemistry factor CF (°C) from Reg Guide 1.99 Rev.2 Table 2.
///
/// Approximate bilinear interpolation over copper content.
pub fn chemistry_factor(copper_wt: f64, nickel_wt: f64) -> f64 {
    let cu = copper_wt.clamp(0.0, 0.40);
    let ni = nickel_wt.clamp(0.0, 1.2);
    // Simplified analytical fit
    let cf_cu = if cu < 0.10 { 0.0 } else { (cu - 0.10) * 250.0 };
    let cf_ni = ni * 20.0;
    (cf_cu + cf_ni).max(0.0)
}

// ---------------------------------------------------------------------------
// Section 7 — Helium embrittlement
// ---------------------------------------------------------------------------

/// Helium embrittlement model for high-fluence structural alloys.
///
/// He accumulates at grain boundaries via (n,α) transmutation reactions.
#[derive(Debug, Clone, Copy)]
pub struct HeliumEmbrittlement {
    /// Grain boundary helium bubble density (m⁻²).
    pub bubble_density: f64,
    /// Mean helium bubble radius (m).
    pub bubble_radius: f64,
    /// Grain size (m).
    pub grain_size: f64,
}

impl HeliumEmbrittlement {
    /// Construct a new helium embrittlement record.
    pub fn new(bubble_density: f64, bubble_radius: f64, grain_size: f64) -> Self {
        Self {
            bubble_density,
            bubble_radius,
            grain_size,
        }
    }

    /// Grain boundary helium coverage fraction (dimensionless, 0–1).
    ///
    /// f_He = N_b · π · r_b²  (projected area fraction on boundary)
    pub fn coverage_fraction(&self) -> f64 {
        (self.bubble_density * PI * self.bubble_radius * self.bubble_radius).min(1.0)
    }

    /// Grain boundary fracture strength reduction factor.
    ///
    /// σ_gb / σ_gb,0 = 1 − A_He · f_He
    /// where A_He ≈ 0.8 (empirical).
    pub fn strength_reduction_factor(&self) -> f64 {
        let a_he = 0.8_f64;
        (1.0 - a_he * self.coverage_fraction()).max(0.0)
    }

    /// Helium concentration (appm) from transmutation.
    ///
    /// C_He = σ_nα · Φ_th · t · N  (simplified, ignoring burnout)
    pub fn helium_appm(
        &self,
        cross_section_m2: f64,
        thermal_flux: f64,
        time_s: f64,
        atom_density: f64,
    ) -> f64 {
        cross_section_m2 * thermal_flux * time_s * atom_density * 1.0e6
    }
}

// ---------------------------------------------------------------------------
// Section 8 — Graphite moderator
// ---------------------------------------------------------------------------

/// Nuclear graphite grade.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum GraphiteGrade {
    /// Pile Grade A (UK Magnox reactors).
    PileGradeA,
    /// IG-110 (Japanese HTGR).
    Ig110,
    /// H-327 (US HTGR — Fort St. Vrain).
    H327,
    /// NBG-18 (European PBMR).
    Nbg18,
}

/// Graphite moderator block.
#[derive(Debug, Clone)]
pub struct GraphiteBlock {
    /// Graphite grade.
    pub grade: GraphiteGrade,
    /// Fast neutron fluence (n/m²).
    pub fluence: f64,
    /// Irradiation temperature (K).
    pub temperature: f64,
    /// Initial bulk density (kg/m³).
    pub initial_density: f64,
}

impl GraphiteBlock {
    /// Construct a new graphite block.
    pub fn new(grade: GraphiteGrade, temperature: f64) -> Self {
        let density = match grade {
            GraphiteGrade::PileGradeA => 1750.0,
            GraphiteGrade::Ig110 => 1770.0,
            GraphiteGrade::H327 => 1740.0,
            GraphiteGrade::Nbg18 => 1850.0,
        };
        Self {
            grade,
            fluence: 0.0,
            temperature,
            initial_density: density,
        }
    }

    /// Unirradiated thermal conductivity (W/m/K).
    pub fn thermal_conductivity_unirradiated(&self) -> f64 {
        match self.grade {
            GraphiteGrade::PileGradeA => 170.0,
            GraphiteGrade::Ig110 => 130.0,
            GraphiteGrade::H327 => 80.0,
            GraphiteGrade::Nbg18 => 140.0,
        }
    }

    /// Irradiated thermal conductivity (W/m/K).
    ///
    /// Rapid initial degradation followed by recovery at very high fluence.
    /// Empirical fit to IAEA-TECDOC-1154 data.
    pub fn thermal_conductivity_irradiated(&self) -> f64 {
        let k0 = self.thermal_conductivity_unirradiated();
        let phi = self.fluence / 1.0e25_f64; // normalise to 10²⁵ n/m²
        // Rapid initial drop: k_irr/k0 = 1/(1 + A·phi) * recovery
        let a = match self.grade {
            GraphiteGrade::PileGradeA => 0.10,
            GraphiteGrade::Ig110 => 0.12,
            GraphiteGrade::H327 => 0.15,
            GraphiteGrade::Nbg18 => 0.11,
        };
        let degradation = 1.0 / (1.0 + a * phi);
        k0 * degradation
    }

    /// Dimensional change (ΔL/L₀) in the with-grain direction.
    ///
    /// Initially contracts (negative), then expands after turnaround.
    pub fn dimensional_change_axial(&self) -> f64 {
        let phi = self.fluence / 1.0e25_f64;
        // Turnaround fluence depends on grade
        let phi_t = match self.grade {
            GraphiteGrade::PileGradeA => 3.0,
            GraphiteGrade::Ig110 => 5.0,
            GraphiteGrade::H327 => 4.0,
            GraphiteGrade::Nbg18 => 4.5,
        };
        let contraction_rate = -0.03_f64;
        let expansion_rate = 0.02_f64;
        if phi < phi_t {
            contraction_rate * phi
        } else {
            contraction_rate * phi_t + expansion_rate * (phi - phi_t)
        }
    }

    /// Young's modulus (Pa) vs irradiation dose.
    ///
    /// Initially increases (~20% at 1 dpa), then decreases at high dose.
    pub fn youngs_modulus(&self) -> f64 {
        let e0 = match self.grade {
            GraphiteGrade::PileGradeA => 8.0e9_f64,
            GraphiteGrade::Ig110 => 9.8e9_f64,
            GraphiteGrade::H327 => 11.0e9_f64,
            GraphiteGrade::Nbg18 => 12.0e9_f64,
        };
        let phi = self.fluence / 1.0e25_f64;
        let factor = 1.0 + 0.20 * phi * (-phi / 2.0).exp();
        e0 * factor
    }

    /// Wigner energy stored in graphite (J/kg) — relevant for Windscale-type incidents.
    ///
    /// Based on Kennedy (1950) empirical correlation.
    pub fn wigner_energy_j_per_kg(&self) -> f64 {
        // Saturation value ~2700 J/kg at low temperature, released above ~200°C
        let e_sat = 2700.0_f64;
        let phi = self.fluence / 1.0e25_f64;
        let t_factor = if self.temperature < 473.0 {
            1.0 - (self.temperature - 300.0) / 173.0 * 0.5
        } else {
            0.5
        };
        e_sat * (1.0 - (-phi * 0.5).exp()) * t_factor.clamp(0.0, 1.0)
    }
}

// ---------------------------------------------------------------------------
// Section 9 — Sodium fast reactor structural materials
// ---------------------------------------------------------------------------

/// Sodium fast reactor structural steel grade.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum SfrSteelGrade {
    /// 316 austenitic stainless steel (AISI 316).
    Ss316,
    /// HT-9 ferritic-martensitic steel (12% Cr).
    Ht9,
    /// D9 austenitic (Ti-modified 316 SS).
    D9,
    /// EP-450 Russian F/M steel.
    Ep450,
}

/// Sodium fast reactor structural material.
#[derive(Debug, Clone)]
pub struct SfrMaterial {
    /// Steel grade.
    pub grade: SfrSteelGrade,
    /// Operating temperature (K).
    pub temperature: f64,
    /// Cumulative dose (dpa).
    pub dose_dpa: f64,
}

impl SfrMaterial {
    /// Construct a new SFR structural material record.
    pub fn new(grade: SfrSteelGrade, temperature: f64) -> Self {
        Self {
            grade,
            temperature,
            dose_dpa: 0.0,
        }
    }

    /// Unirradiated tensile strength (MPa).
    pub fn tensile_strength_unirradiated(&self) -> f64 {
        match self.grade {
            SfrSteelGrade::Ss316 => 515.0,
            SfrSteelGrade::Ht9 => 690.0,
            SfrSteelGrade::D9 => 550.0,
            SfrSteelGrade::Ep450 => 680.0,
        }
    }

    /// Irradiation hardening — yield strength increase (MPa).
    ///
    /// Δσ_y = A · (dpa)^0.5  (dispersed barrier hardening model).
    pub fn irradiation_hardening_mpa(&self) -> f64 {
        let a = match self.grade {
            SfrSteelGrade::Ss316 => 50.0,
            SfrSteelGrade::Ht9 => 30.0, // F/M steels saturate faster
            SfrSteelGrade::D9 => 45.0,
            SfrSteelGrade::Ep450 => 28.0,
        };
        a * self.dose_dpa.sqrt()
    }

    /// Void swelling (%ΔV/V₀) for austenitic steels — Garner correlation.
    pub fn void_swelling_pct(&self) -> f64 {
        match self.grade {
            SfrSteelGrade::Ss316 | SfrSteelGrade::D9 => {
                let params = VoidSwellingParams::default();
                garner_void_swelling(self.dose_dpa, self.temperature, &params) * 100.0
            }
            // F/M steels have negligible void swelling
            SfrSteelGrade::Ht9 | SfrSteelGrade::Ep450 => {
                0.002 * self.dose_dpa // < 0.2% even at 100 dpa
            }
        }
    }

    /// Thermal creep rate (s⁻¹) at given stress (MPa) using Norton power law.
    pub fn thermal_creep_rate(&self, stress_mpa: f64) -> f64 {
        // ε̇ = A · σ^n · exp(−Q/RT)
        let (a, n, q) = match self.grade {
            SfrSteelGrade::Ss316 => (1.5e-32_f64, 5.0_f64, 280_000.0_f64),
            SfrSteelGrade::Ht9 => (3.0e-28_f64, 4.5_f64, 260_000.0_f64),
            SfrSteelGrade::D9 => (1.2e-32_f64, 5.0_f64, 275_000.0_f64),
            SfrSteelGrade::Ep450 => (2.5e-28_f64, 4.5_f64, 255_000.0_f64),
        };
        let r = 8.314_f64;
        a * stress_mpa.powf(n) * (-(q) / (r * self.temperature)).exp()
    }

    /// Fatigue life (cycles) at given strain amplitude using Coffin-Manson.
    ///
    /// 2N_f = (Δε / ε_f')^(1/c)  where c ≈ −0.6 (typical austenitic SS).
    pub fn fatigue_life_cycles(&self, strain_amplitude: f64) -> f64 {
        let (eps_f, c) = match self.grade {
            SfrSteelGrade::Ss316 | SfrSteelGrade::D9 => (0.3_f64, -0.60_f64),
            SfrSteelGrade::Ht9 | SfrSteelGrade::Ep450 => (0.25_f64, -0.57_f64),
        };
        if strain_amplitude < 1e-9 {
            return f64::MAX;
        }
        let two_nf = (strain_amplitude / eps_f).powf(1.0 / c);
        (two_nf / 2.0).max(1.0)
    }

    /// Creep-fatigue damage parameter D_cf = D_c + D_f.
    ///
    /// D_c = Σ(Δt/t_r)  and D_f = Σ(n/N_f).
    /// Returns true if D_cf ≥ 1.0 (failure predicted).
    pub fn creep_fatigue_interaction(
        &self,
        delta_t_s: f64,
        rupture_time_s: f64,
        cycles: f64,
        n_f: f64,
    ) -> bool {
        let d_c = delta_t_s / rupture_time_s.max(1e-9);
        let d_f = cycles / n_f.max(1.0);
        (d_c + d_f) >= 1.0
    }
}

// ---------------------------------------------------------------------------
// Section 10 — Nuclear fuel cycle and transmutation helpers
// ---------------------------------------------------------------------------

/// Simple Bateman-style two-nuclide decay chain.
///
/// Parent nuclide A decays to stable daughter B with decay constant λ.
#[derive(Debug, Clone, Copy)]
pub struct DecayChain {
    /// Decay constant λ (s⁻¹).
    pub lambda: f64,
    /// Initial number of parent atoms.
    pub n0_parent: f64,
}

impl DecayChain {
    /// Construct from half-life (s).
    pub fn from_half_life(t_half_s: f64, n0_parent: f64) -> Self {
        Self {
            lambda: 2.0_f64.ln() / t_half_s,
            n0_parent,
        }
    }

    /// Parent atom count at time `t` (s).
    pub fn parent_count(&self, t: f64) -> f64 {
        self.n0_parent * (-self.lambda * t).exp()
    }

    /// Daughter atom count at time `t` (s) assuming zero initial daughters.
    pub fn daughter_count(&self, t: f64) -> f64 {
        self.n0_parent * (1.0 - (-self.lambda * t).exp())
    }

    /// Radioactivity (Bq) at time `t`.
    pub fn activity_bq(&self, t: f64) -> f64 {
        self.lambda * self.parent_count(t)
    }

    /// Specific activity (Bq/kg) given molar mass (g/mol).
    pub fn specific_activity(&self, molar_mass_g_per_mol: f64, n0_mass_kg: f64) -> f64 {
        let avogadro = 6.022e23_f64;
        let n0 = n0_mass_kg * 1000.0 / molar_mass_g_per_mol * avogadro;
        let chain = Self {
            lambda: self.lambda,
            n0_parent: n0,
        };
        chain.activity_bq(0.0)
    }
}

/// Compute reactivity worth of a control rod absorber (simplified 1-group).
///
/// Uses Hurst-type formula: ρ = −Σ_a_rod / (ν · Σ_f).
pub fn control_rod_worth(sigma_a_rod_m2: f64, nu_sigma_f_m2: f64) -> f64 {
    if nu_sigma_f_m2 < 1e-40 {
        return 0.0;
    }
    -sigma_a_rod_m2 / nu_sigma_f_m2
}

// ---------------------------------------------------------------------------
// Section 11 — Pellet-Cladding Interaction (PCI)
// ---------------------------------------------------------------------------

/// Pellet-cladding interaction parameters.
#[derive(Debug, Clone, Copy)]
pub struct PciParams {
    /// Pellet outer radius at operating conditions (m).
    pub pellet_radius: f64,
    /// Cladding inner radius (m).
    pub clad_inner_radius: f64,
    /// Cladding elastic modulus (Pa).
    pub clad_youngs_modulus: f64,
    /// Cladding Poisson's ratio.
    pub clad_poisson: f64,
    /// Cladding wall thickness (m).
    pub clad_thickness: f64,
}

impl Default for PciParams {
    fn default() -> Self {
        Self {
            pellet_radius: 4.1e-3,
            clad_inner_radius: 4.18e-3,
            clad_youngs_modulus: 75e9,
            clad_poisson: 0.37,
            clad_thickness: 0.57e-3,
        }
    }
}

/// Contact pressure (Pa) between pellet and cladding due to swelling.
///
/// Uses Lamé cylinder theory for the cladding.
pub fn pci_contact_pressure(params: &PciParams, pellet_swelling_fraction: f64) -> f64 {
    let delta = params.pellet_radius * pellet_swelling_fraction
        - (params.clad_inner_radius - params.pellet_radius);
    if delta <= 0.0 {
        return 0.0;
    }
    // Thin-shell approximation
    let e = params.clad_youngs_modulus;
    let nu = params.clad_poisson;
    let ri = params.clad_inner_radius;
    let t = params.clad_thickness;
    // p = E · δ · t / ((1 − ν) · ri²)
    e * delta * t / ((1.0 - nu) * ri * ri)
}

// ---------------------------------------------------------------------------
// Section 12 — Burnup conversion helpers
// ---------------------------------------------------------------------------

/// Convert specific power (W/kgU) and irradiation time (days) to burnup (MWd/kgU).
pub fn power_to_burnup(specific_power_w_per_kgu: f64, time_days: f64) -> f64 {
    specific_power_w_per_kgu * time_days / 1e6
}

/// Convert burnup (MWd/kgU) to dpa for fuel cladding (approximate).
///
/// For Zircaloy in PWR: ~1 dpa per 10 MWd/kgU of local burnup.
pub fn burnup_to_dpa_clad(burnup_mwd_per_kgu: f64) -> f64 {
    burnup_mwd_per_kgu / 10.0
}

/// Fission density (fissions/m³) from burnup.
pub fn fission_density(burnup_mwd_per_kgu: f64, uo2_density_kg_per_m3: f64) -> f64 {
    // 1 MWd = 86400 MJ; 1 fission ~ 3.2e-11 J
    let u_fraction = 238.03 / (238.03 + 32.0); // UO₂ heavy metal fraction
    let mhm_per_m3 = uo2_density_kg_per_m3 * u_fraction;
    // burnup in J/kgHM
    let burnup_j = burnup_mwd_per_kgu * 1e6 * 86400.0;
    mhm_per_m3 * burnup_j / 3.2e-11
}

// ---------------------------------------------------------------------------
// Section 13 — Unit tests
// ---------------------------------------------------------------------------

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

    #[test]
    fn test_zircaloy_clad_construction() {
        let clad = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        assert_eq!(clad.grade, ZircaloyGrade::Zircaloy4);
        assert_eq!(clad.oxide_thickness, 0.0);
    }

    #[test]
    fn test_cubic_oxide_weight_gain_zero_time() {
        let clad = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        let wg = clad.oxide_weight_gain_cubic(0.0);
        assert!(wg.abs() < 1e-15, "zero time → zero weight gain");
    }

    #[test]
    fn test_cubic_oxide_weight_gain_increases_with_time() {
        let clad = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        let wg1 = clad.oxide_weight_gain_cubic(1e6);
        let wg2 = clad.oxide_weight_gain_cubic(1e7);
        assert!(wg2 > wg1, "longer time should give greater weight gain");
    }

    #[test]
    fn test_hydrogen_pickup_fraction_range() {
        let clad = ZircaloyClad::new(ZircaloyGrade::M5, 9.5e-3, 0.57e-3, 580.0);
        let f = clad.hydrogen_pickup_fraction();
        assert!(
            (0.0..=1.0).contains(&f),
            "pickup fraction must be [0,1]: {}",
            f
        );
    }

    #[test]
    fn test_hydrogen_pickup_m5_lower_than_zry4() {
        let clad4 = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        let cladm5 = ZircaloyClad::new(ZircaloyGrade::M5, 9.5e-3, 0.57e-3, 600.0);
        assert!(
            cladm5.hydrogen_pickup_fraction() < clad4.hydrogen_pickup_fraction(),
            "M5 should have lower pickup than Zry-4"
        );
    }

    #[test]
    fn test_yield_strength_increases_with_fluence() {
        let mut clad = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        let ys0 = clad.yield_strength();
        clad.fluence = 1e25;
        let ys1 = clad.yield_strength();
        assert!(ys1 > ys0, "irradiation should harden Zircaloy");
    }

    #[test]
    fn test_irradiation_creep_rate_proportional_to_flux() {
        let clad = ZircaloyClad::new(ZircaloyGrade::Zircaloy4, 9.5e-3, 0.57e-3, 600.0);
        let r1 = clad.irradiation_creep_rate(100e6, 1e13);
        let r2 = clad.irradiation_creep_rate(100e6, 2e13);
        assert!(
            (r2 / r1 - 2.0).abs() < 1e-9,
            "creep rate should scale linearly with flux"
        );
    }

    #[test]
    fn test_uo2_thermal_conductivity_positive() {
        let pellet = Uo2Pellet::new(8.19e-3, 10e-3, 0.0315);
        let k = pellet.thermal_conductivity(900.0);
        assert!(k > 0.0, "thermal conductivity must be positive, got {}", k);
    }

    #[test]
    fn test_uo2_thermal_conductivity_decreases_with_burnup() {
        let mut pellet = Uo2Pellet::new(8.19e-3, 10e-3, 0.0315);
        let k0 = pellet.thermal_conductivity(900.0);
        pellet.burnup = 50.0;
        let k50 = pellet.thermal_conductivity(900.0);
        assert!(k50 < k0, "conductivity should decrease with burnup");
    }

    #[test]
    fn test_fission_gas_release_zero_at_t0() {
        let pellet = Uo2Pellet::new(8.19e-3, 10e-3, 0.0315);
        let fgr = pellet.fission_gas_release_booth(1200.0, 0.0);
        assert!(fgr.abs() < 1e-9, "FGR at t=0 must be 0");
    }

    #[test]
    fn test_fission_gas_release_bounded() {
        let pellet = Uo2Pellet::new(8.19e-3, 10e-3, 0.0315);
        let fgr = pellet.fission_gas_release_booth(1400.0, 1e10);
        assert!((0.0..=1.0).contains(&fgr), "FGR must be in [0,1]: {}", fgr);
    }

    #[test]
    fn test_pellet_swelling_increases_with_burnup() {
        let mut p = Uo2Pellet::new(8.19e-3, 10e-3, 0.0315);
        p.burnup = 0.0;
        let s0 = p.fission_product_swelling();
        p.burnup = 50.0;
        let s50 = p.fission_product_swelling();
        assert!(s50 > s0, "swelling increases with burnup");
    }

    #[test]
    fn test_garner_void_swelling_below_incubation_zero() {
        let params = VoidSwellingParams::default();
        let s = garner_void_swelling(5.0, 773.0, &params);
        assert_eq!(s, 0.0, "no swelling before incubation dose");
    }

    #[test]
    fn test_garner_void_swelling_above_incubation_positive() {
        let params = VoidSwellingParams::default();
        let s = garner_void_swelling(50.0, 773.0, &params);
        assert!(s > 0.0, "swelling should be positive above incubation dose");
    }

    #[test]
    fn test_garner_void_swelling_peak_temperature() {
        let params = VoidSwellingParams::default();
        let s_peak = garner_void_swelling(50.0, params.peak_temperature_k, &params);
        let s_off = garner_void_swelling(50.0, params.peak_temperature_k + 100.0, &params);
        assert!(
            s_peak > s_off,
            "peak temperature should give maximum swelling"
        );
    }

    #[test]
    fn test_void_number_density_zero_for_zero_swelling() {
        let n = void_number_density(0.0, 5e-9);
        assert_eq!(n, 0.0);
    }

    #[test]
    fn test_irradiation_creep_additive() {
        let r = irradiation_creep_rate(200e6, 773.0, 1e-8, 3e-27, 1e-33, 5.0, 280000.0);
        assert!(r > 0.0, "creep rate must be positive");
    }

    #[test]
    fn test_stress_relaxation_decreases_with_dose() {
        let sigma = irradiation_stress_relaxation(200e6, 0.01, 100.0);
        assert!(sigma < 200e6, "stress should relax with dose");
        assert!(sigma > 0.0, "stress must remain positive");
    }

    #[test]
    fn test_ris_concentration_bounded() {
        let params = RisParams::default();
        let c = ris_grain_boundary_concentration(&params);
        assert!(
            (0.0..=1.0).contains(&c),
            "RIS concentration out of [0,1]: {}",
            c
        );
    }

    #[test]
    fn test_rpv_charpy_use_drop_low_fluence() {
        let charpy = RpvCharpy::new(
            RpvSteelGrade::A508Class3,
            100.0,
            -20.0,
            1e17,
            0.05,
            0.7,
            0.01,
        );
        let drop = charpy.use_drop();
        assert_eq!(drop, 0.0, "no USE drop below fluence threshold");
    }

    #[test]
    fn test_rpv_charpy_use_irradiated_non_negative() {
        let charpy = RpvCharpy::new(
            RpvSteelGrade::A533GradeB,
            100.0,
            -30.0,
            1e20,
            0.20,
            0.8,
            0.02,
        );
        assert!(
            charpy.use_irradiated() >= 0.0,
            "irradiated USE must be non-negative"
        );
    }

    #[test]
    fn test_rpv_rtndt_shift_positive_copper() {
        let charpy = RpvCharpy::new(
            RpvSteelGrade::A508Class3,
            100.0,
            -20.0,
            5e19,
            0.20,
            0.7,
            0.01,
        );
        let shift = charpy.rtndt_shift();
        assert!(shift > 0.0, "RTNDT shift must be positive with high copper");
    }

    #[test]
    fn test_fracture_toughness_increases_with_temperature() {
        let charpy = RpvCharpy::new(
            RpvSteelGrade::A508Class3,
            100.0,
            -20.0,
            1e19,
            0.1,
            0.7,
            0.01,
        );
        let kic_cold = charpy.fracture_toughness_kic(-100.0);
        let kic_warm = charpy.fracture_toughness_kic(50.0);
        assert!(
            kic_warm > kic_cold,
            "fracture toughness increases above RTNDT"
        );
    }

    #[test]
    fn test_helium_coverage_fraction_bounded() {
        let he = HeliumEmbrittlement::new(1e16, 2e-9, 20e-6);
        let f = he.coverage_fraction();
        assert!(
            (0.0..=1.0).contains(&f),
            "coverage fraction out of [0,1]: {}",
            f
        );
    }

    #[test]
    fn test_helium_strength_reduction_non_negative() {
        let he = HeliumEmbrittlement::new(1e16, 5e-9, 20e-6);
        let sr = he.strength_reduction_factor();
        assert!(
            sr >= 0.0,
            "strength reduction factor must be non-negative: {}",
            sr
        );
    }

    #[test]
    fn test_graphite_conductivity_decreases_with_fluence() {
        let mut g = GraphiteBlock::new(GraphiteGrade::Ig110, 673.0);
        let k0 = g.thermal_conductivity_irradiated();
        g.fluence = 1e26;
        let k_hi = g.thermal_conductivity_irradiated();
        assert!(
            k_hi < k0,
            "graphite conductivity should decrease with fluence"
        );
    }

    #[test]
    fn test_graphite_wigner_zero_at_zero_fluence() {
        let g = GraphiteBlock::new(GraphiteGrade::PileGradeA, 400.0);
        let w = g.wigner_energy_j_per_kg();
        assert!(w.abs() < 1e-9, "zero fluence → zero Wigner energy");
    }

    #[test]
    fn test_sfr_void_swelling_316_positive() {
        let mut mat = SfrMaterial::new(SfrSteelGrade::Ss316, 773.0);
        mat.dose_dpa = 50.0;
        let s = mat.void_swelling_pct();
        assert!(s >= 0.0, "void swelling must be non-negative");
    }

    #[test]
    fn test_sfr_ht9_swelling_lower_than_316() {
        let mut mat_316 = SfrMaterial::new(SfrSteelGrade::Ss316, 773.0);
        let mut mat_ht9 = SfrMaterial::new(SfrSteelGrade::Ht9, 773.0);
        mat_316.dose_dpa = 100.0;
        mat_ht9.dose_dpa = 100.0;
        assert!(
            mat_ht9.void_swelling_pct() < mat_316.void_swelling_pct(),
            "HT-9 should swell less than 316 SS"
        );
    }

    #[test]
    fn test_thermal_creep_rate_increases_with_temperature() {
        let mat_lo = SfrMaterial::new(SfrSteelGrade::Ss316, 700.0);
        let mat_hi = SfrMaterial::new(SfrSteelGrade::Ss316, 900.0);
        let r_lo = mat_lo.thermal_creep_rate(100.0);
        let r_hi = mat_hi.thermal_creep_rate(100.0);
        assert!(r_hi > r_lo, "creep rate should increase with temperature");
    }

    #[test]
    fn test_fatigue_life_decreases_with_strain() {
        let mat = SfrMaterial::new(SfrSteelGrade::Ss316, 773.0);
        let n1 = mat.fatigue_life_cycles(0.001);
        let n2 = mat.fatigue_life_cycles(0.01);
        assert!(
            n1 > n2,
            "higher strain amplitude gives shorter fatigue life"
        );
    }

    #[test]
    fn test_decay_chain_parent_count_at_zero() {
        let dc = DecayChain::from_half_life(1e9, 1e20);
        assert!((dc.parent_count(0.0) - 1e20).abs() < 1.0);
    }

    #[test]
    fn test_decay_chain_conservation() {
        let dc = DecayChain::from_half_life(3600.0, 1e20);
        let t = 3600.0 * 3.0;
        let total = dc.parent_count(t) + dc.daughter_count(t);
        assert!(
            (total - 1e20).abs() / 1e20 < 1e-9,
            "atom count must be conserved"
        );
    }

    #[test]
    fn test_pci_contact_pressure_zero_gap() {
        let params = PciParams::default();
        let p = pci_contact_pressure(&params, 0.0);
        assert_eq!(p, 0.0, "no contact pressure when gap is open");
    }

    #[test]
    fn test_pci_contact_pressure_positive_on_swelling() {
        let params = PciParams::default();
        let p = pci_contact_pressure(&params, 0.01);
        assert!(p >= 0.0, "contact pressure must be non-negative");
    }

    #[test]
    fn test_burnup_to_dpa_scaling() {
        let dpa = burnup_to_dpa_clad(50.0);
        assert!((dpa - 5.0).abs() < 1e-9, "50 MWd/kgU → 5 dpa");
    }

    #[test]
    fn test_fission_density_positive() {
        let fd = fission_density(50.0, 10400.0);
        assert!(fd > 0.0, "fission density must be positive");
    }

    #[test]
    fn test_chemistry_factor_zero_low_cu() {
        let cf = chemistry_factor(0.05, 0.0);
        assert_eq!(cf, 0.0, "no CF for very low copper");
    }
}