materials_for_mc 0.2.0

// ...existing code...
use crate::config::CONFIG;
use crate::data::ELEMENT_NAMES;
use crate::nuclide::{get_or_load_nuclide, Nuclide};
use crate::utilities::interpolate_linear;
use std::collections::HashMap;
use std::sync::Arc;

/// Represents a heterogeneous collection of nuclides (or elements expanded to
/// their naturally abundant isotopes) along with the material density and
/// nuclear data needed for transport / analysis.
///
/// A `Material` starts empty; users add nuclides with [`Material::add_nuclide`]
/// or elements with [`Material::add_element`] (which expands to isotopes using
/// natural abundances). Density (with units), temperature (string key matching
/// JSON data), and optional volume may then be specified. On‑demand the
/// structure loads JSON nuclide data (through the global [`crate::config::Config`]) and builds
/// a unified energy grid for neutrons so reaction cross sections for different
/// nuclides can be interpolated on a common axis.
///
/// Key cached members:
/// * `unified_energy_grid_neutron` – lazily constructed common energy grid.
/// * `macroscopic_xs_neutron` – map MT -> Σ(E) on the unified grid.
/// * `macroscopic_xs_neutron_total_by_nuclide` – optional per‑nuclide Σ_t(E) when
///   requested (used for sampling interacting nuclides).
///
/// Typical workflow:
/// 1. Create with [`Material::new`].
/// 2. Populate composition (nuclides or elements) and set density.
/// 3. Optionally set temperature (clears caches if changed).
/// 4. Call cross section building methods (e.g.
///    [`Material::calculate_macroscopic_xs`]) or sampling utilities.
#[derive(Debug, Clone)]
pub struct Material {
    /// Optional name of the material
    pub name: Option<String>,
    /// Composition of the material as a map of nuclide names to their atomic fractions
    pub nuclides: HashMap<String, f64>,
    /// Density of the material in g/cm³
    pub density: Option<f64>,
    /// Density unit (default: g/cm³)
    pub density_units: String,
    /// Volume of the material in cm³
    pub volume: Option<f64>,
    /// Temperature of the material in K
    pub temperature: String,
    /// Loaded nuclide data (name -> `Arc<Nuclide>`) shared for this material instance
    pub nuclide_data: HashMap<String, Arc<Nuclide>>,
    /// Macroscopic cross sections for different MT numbers (neutron only for now)
    /// Map of MT number (i32) -> cross sections
    pub macroscopic_xs_neutron: HashMap<i32, Vec<f64>>,
    /// Unified energy grid for neutrons
    pub unified_energy_grid_neutron: Vec<f64>,
    /// Optional: Per-nuclide macroscopic total cross section (MT=1) on the unified grid
    /// Map: nuclide name -> `Vec<f64>` (same length as unified_energy_grid_neutron)
    pub macroscopic_xs_neutron_total_by_nuclide: Option<HashMap<String, Vec<f64>>>,
}

impl Material {
    pub fn new() -> Self {
        Material {
            name: None,
            nuclides: HashMap::new(),
            density: None,
            density_units: String::from("g/cm3"),
            volume: None,                     // Initialize volume as None
            temperature: String::from("294"), // Default temperature in K (room temperature)
            nuclide_data: HashMap::new(),
            macroscopic_xs_neutron: HashMap::new(),
            unified_energy_grid_neutron: Vec::new(),
            macroscopic_xs_neutron_total_by_nuclide: None,
        }
    }

    /// Set the name of the material
    pub fn set_name(&mut self, name: impl Into<String>) {
        self.name = Some(name.into());
    }

    /// Get the name of the material
    pub fn get_name(&self) -> Option<&str> {
        self.name.as_deref()
    }

    /// Clear all cached cross section data
    fn invalidate_xs_cache(&mut self) {
        self.macroscopic_xs_neutron.clear();
        self.macroscopic_xs_neutron_total_by_nuclide = None;
        self.unified_energy_grid_neutron.clear();
    }

    pub fn add_nuclide(&mut self, nuclide: impl AsRef<str>, fraction: f64) -> Result<(), String> {
        if fraction < 0.0 {
            return Err(String::from("Fraction cannot be negative"));
        }

        self.nuclides
            .insert(String::from(nuclide.as_ref()), fraction);
        
        // Clear cached data since composition changed
        self.invalidate_xs_cache();
        Ok(())
    }

    /// Sample the distance to the next collision for a neutron at the given energy.
    /// Uses the total macroscopic cross section (MT=1).
    /// Returns None if the cross section is zero or not available.
    pub fn sample_distance_to_collision<R: rand::Rng + ?Sized>(
        &self,
        energy: f64,
        rng: &mut R,
    ) -> Option<f64> {
        let xs_vec = self.macroscopic_xs_neutron.get(&1)?;
        let sigma_t =
            crate::utilities::interpolate_linear(&self.unified_energy_grid_neutron, xs_vec, energy);
        if sigma_t <= 0.0 {
            return None;
        }
        let xi: f64 = rng.gen_range(0.0..1.0);
        Some(-xi.ln() / sigma_t)
    }

    pub fn set_density(&mut self, unit: impl AsRef<str>, value: f64) -> Result<(), String> {
        if value <= 0.0 {
            return Err(String::from("Density must be positive"));
        }

        self.density = Some(value);
        self.density_units = String::from(unit.as_ref());
        
        // Clear cached data since density affects macroscopic cross sections
        self.invalidate_xs_cache();
        Ok(())
    }

    pub fn volume(&mut self, value: Option<f64>) -> Result<Option<f64>, String> {
        if let Some(v) = value {
            if v <= 0.0 {
                return Err(String::from("Volume must be positive"));
            }
            self.volume = Some(v);
        }
        Ok(self.volume)
    }

    pub fn set_temperature(&mut self, temperature: impl AsRef<str>) {
        self.temperature = String::from(temperature.as_ref());
        // Clear cached data that depends on temperature
        self.unified_energy_grid_neutron.clear();
        self.macroscopic_xs_neutron.clear();
    }

    pub fn get_nuclides(&self) -> Vec<String> {
        let mut nuclides: Vec<String> = self.nuclides.keys().cloned().collect();
        nuclides.sort(); // Sort alphabetically for consistent output
        nuclides
    }

    /// Read nuclide data from JSON files for this material
    pub fn read_nuclides_from_json(
        &mut self,
        nuclide_json_map: &HashMap<String, String>,
    ) -> Result<(), Box<dyn std::error::Error>> {
        // Collect needed nuclide names
        let mut nuclide_names: Vec<String> = self.nuclides.keys().cloned().collect();
        nuclide_names.sort(); // ensure deterministic alphabetical load order

        // Build merged source map: explicit entries override, missing filled from CONFIG
        let mut merged: HashMap<String, String> = HashMap::new();

        // Start with global config entries for required nuclides, with proper error handling
        let cfg = crate::config::CONFIG
            .lock()
            .unwrap_or_else(|poisoned| poisoned.into_inner());

        for n in &nuclide_names {
            if let Some(p) = cfg.get_cross_section(n) {
                merged.insert(n.clone(), p.clone());
            }
        }
        drop(cfg);

        // Override with any provided mapping entries (even if extra keys not in composition)
        for (k, v) in nuclide_json_map {
            merged.insert(k.clone(), v.clone());
        }
        let source_map: &HashMap<String, String> = &merged;

        // Load nuclides using the centralized function in the nuclide module
        use std::collections::HashSet;
        let mut temp_set: HashSet<String> = HashSet::new();
        temp_set.insert(self.temperature.clone());

        for nuclide_name in nuclide_names {
            let nuclide = get_or_load_nuclide(&nuclide_name, source_map, Some(&temp_set))?;
            self.nuclide_data.insert(nuclide_name, nuclide);
        }

        // Clear cached data since new nuclear data affects cross sections
        self.invalidate_xs_cache();
        Ok(())
    }

    /// Read nuclide data from either a JSON mapping or a keyword string
    pub fn read_nuclides_from_json_or_keyword(
        &mut self,
        source: &str,
    ) -> Result<(), Box<dyn std::error::Error>> {
        if crate::url_cache::is_keyword(source) {
            // It's a keyword - apply to all nuclides in this material
            let mut keyword_map = HashMap::new();
            for nuclide_name in self.nuclides.keys() {
                keyword_map.insert(nuclide_name.clone(), source.to_string());
            }
            self.read_nuclides_from_json(&keyword_map)?;
        } else {
            // Treat as a single nuclide with a path (legacy behavior)
            let mut single_map = HashMap::new();
            single_map.insert(source.to_string(), source.to_string());
            self.read_nuclides_from_json(&single_map)?;
        }
        Ok(())
    }



    /// Read nuclides from a keyword string that will be applied to all nuclides in this material
    pub fn read_nuclides_from_json_keyword(
        &mut self,
        keyword: &str,
    ) -> Result<(), Box<dyn std::error::Error>> {
        self.read_nuclides_from_json_or_keyword(keyword)
    }

    /// Read nuclides from either a HashMap or handle None case
    pub fn read_nuclides_from_optional_map(
        &mut self,
        map: Option<&HashMap<String, String>>,
    ) -> Result<(), Box<dyn std::error::Error>> {
        match map {
            Some(m) => self.read_nuclides_from_json(m),
            None => {
                let empty_map = HashMap::new();
                self.read_nuclides_from_json(&empty_map)
            }
        }
    }

    /// Load nuclear data from extracted input data (pure Rust, no PyO3 dependencies)
    pub fn load_nuclear_data_from_input(
        &mut self,
        dict_data: Option<HashMap<String, String>>,
        keyword_data: Option<String>,
    ) -> Result<(), Box<dyn std::error::Error>> {
        if let Some(map) = dict_data {
            self.read_nuclides_from_json(&map)
        } else if let Some(keyword) = keyword_data {
            self.read_nuclides_from_json_keyword(&keyword)
        } else {
            let empty_map = HashMap::new();
            self.read_nuclides_from_json(&empty_map)
        }
    }

    /// Read nuclides from a string keyword (for Python wrapper)
    pub fn read_nuclides_from_string(
        &mut self,
        keyword: &str,
    ) -> Result<(), Box<dyn std::error::Error>> {
        self.read_nuclides_from_json_keyword(keyword)
    }

    /// Read nuclides from a HashMap (for Python wrapper)
    pub fn read_nuclides_from_map(
        &mut self,
        map: &HashMap<String, String>,
    ) -> Result<(), Box<dyn std::error::Error>> {
        self.read_nuclides_from_json(map)
    }

    /// Read nuclides with no input - use defaults (for Python wrapper)
    pub fn read_nuclides_from_none(
        &mut self,
    ) -> Result<(), Box<dyn std::error::Error>> {
        let empty_map = HashMap::new();
        self.read_nuclides_from_json(&empty_map)
    }

    /// Directly load a nuclide from a JSON string (e.g., for WASM in-memory usage) and insert into nuclide_data.
    /// If the nuclide already exists it will be overwritten.
    pub fn load_nuclide_from_json_str(
        &mut self,
        nuclide_name: &str,
        json_content: &str,
    ) -> Result<(), Box<dyn std::error::Error>> {
        let nuclide = crate::nuclide::read_nuclide_from_json_str(json_content)?;
        self.nuclide_data
            .insert(nuclide_name.to_string(), Arc::new(nuclide));
        Ok(())
    }

    /// Ensure all nuclides are loaded, using the global configuration if needed
    fn ensure_nuclides_loaded(&mut self) -> Result<(), Box<dyn std::error::Error>> {
        let nuclide_names: Vec<String> = self
            .nuclides
            .keys()
            .filter(|name| !self.nuclide_data.contains_key(*name))
            .cloned()
            .collect();

        if nuclide_names.is_empty() {
            return Ok(());
        }

        // Get the global configuration with proper error handling
        let config = CONFIG
            .lock()
            .unwrap_or_else(|poisoned| poisoned.into_inner());

        // Load any missing nuclides
        for nuclide_name in nuclide_names {
            use std::collections::HashSet;
            let mut temps = HashSet::new();
            temps.insert(self.temperature.clone());
            
            // Build a temporary source map with the global default fallback
            let mut source_map = HashMap::new();
            if let Some(path) = config.get_cross_section(&nuclide_name) {
                source_map.insert(nuclide_name.clone(), path);
            }
            
            match get_or_load_nuclide(&nuclide_name, &source_map, Some(&temps)) {
                Ok(nuclide) => {
                    self.nuclide_data.insert(nuclide_name.clone(), nuclide);
                }
                Err(e) => {
                    return Err(format!("Failed to load nuclide '{}': {}", nuclide_name, e).into());
                }
            }
        }

        // Clear cached data since new nuclear data affects cross sections
        self.invalidate_xs_cache();
        Ok(())
    }

    /// Build a unified energy grid for all nuclides for neutrons across all MT reactions
    /// This method also stores the result in the material's unified_energy_grid_neutron property
    pub fn unified_energy_grid_neutron(&mut self) -> Vec<f64> {
        // Ensure nuclides are loaded before proceeding
        if let Err(e) = self.ensure_nuclides_loaded() {
            panic!("Error loading nuclides: {}", e);
        }

        // Check if we already have this grid in the cache
        if !self.unified_energy_grid_neutron.is_empty() {
            return self.unified_energy_grid_neutron.clone();
        }

        // If not cached, build the grid
        let mut all_energies = Vec::new();
        let temperature = &self.temperature;
        let _particle = "neutron"; // This is now specifically for neutrons

        for nuclide in self.nuclides.keys() {
            if let Some(nuclide_data) = self.nuclide_data.get(nuclide) {
                // Check if there's a top-level energy grid
                if let Some(energy_map) = &nuclide_data.energy {
                    if let Some(energy_grid) = energy_map.get(temperature) {
                        all_energies.extend(energy_grid);
                    }
                }
            }
        }

        // Sort and deduplicate
        all_energies.sort_by(|a: &f64, b: &f64| a.partial_cmp(b).unwrap());
        all_energies.dedup_by(|a, b| (*a - *b).abs() < 1e-12);

        // Cache the result
        self.unified_energy_grid_neutron = all_energies.clone();

        all_energies
    }

    /// Calculate microscopic cross sections for neutrons on the unified energy grid
    ///
    /// This method interpolates the microscopic cross sections for each nuclide
    /// onto the unified energy grid for all available MT reactions, or only for the specified MTs if provided.
    /// If mt_filter is Some, only those MTs will be included (by int match).
    /// Returns a nested HashMap: nuclide -> mt -> cross_section values
    pub fn calculate_microscopic_xs_neutron(
        &mut self,
        mt_filter: Option<&Vec<i32>>,
    ) -> HashMap<String, HashMap<i32, Vec<f64>>> {
        // Ensure nuclides are loaded before proceeding
        if let Err(e) = self.ensure_nuclides_loaded() {
            panic!("Error loading nuclides: {}", e);
        }

        let grid = self.unified_energy_grid_neutron();
        let mut micro_xs: HashMap<String, HashMap<i32, Vec<f64>>> = HashMap::new();
        let temperature = &self.temperature;
        // Unified logic: iterate all reactions; if a filter is provided skip non-matching MTs.
        let mt_set_opt: Option<std::collections::HashSet<i32>> = mt_filter.map(|v| v.iter().copied().collect());
        for nuclide_name in self.nuclides.keys() {
            if let Some(nuclide_data) = self.nuclide_data.get(nuclide_name) {
                let mut nuclide_reactions_map: HashMap<i32, Vec<f64>> = HashMap::new();
                if let Some(temp_reactions) = nuclide_data.reactions.get(temperature) {
                    if let Some(energy_map) = &nuclide_data.energy {
                        if let Some(energy_grid) = energy_map.get(temperature) {
                            for (&mt, reaction) in temp_reactions {
                                if let Some(ref set) = mt_set_opt {
                                    if !set.contains(&mt) { continue; }
                                }
                                let threshold_idx = reaction.threshold_idx;
                                if threshold_idx < energy_grid.len() {
                                    let reaction_energy = &energy_grid[threshold_idx..];
                                    if reaction.cross_section.len() == reaction_energy.len() {
                                        let mut xs_values = Vec::with_capacity(grid.len());
                                        for &grid_energy in &grid {
                                            if grid_energy < reaction_energy[0] {
                                                xs_values.push(0.0);
                                            } else {
                                                let xs = interpolate_linear(
                                                    &reaction_energy,
                                                    &reaction.cross_section,
                                                    grid_energy,
                                                );
                                                xs_values.push(xs);
                                            }
                                        }
                                        nuclide_reactions_map.insert(mt, xs_values);
                                    }
                                }
                            }
                        }
                    }
                }
                if !nuclide_reactions_map.is_empty() {
                    micro_xs.insert(nuclide_name.clone(), nuclide_reactions_map);
                }
            }
        }
        micro_xs
    }

    /// Calculate macroscopic cross sections for neutrons on the unified energy grid
    ///
    /// This method calculates the total macroscopic cross section by:
    /// 1. Interpolating the microscopic cross sections onto the unified grid
    /// 2. Multiplying by atom density for each nuclide
    /// 3. Summing over all nuclides
    /// If mt_filter is Some, only those MTs will be included (by string match).
    /// If by_nuclide is true, populates the struct field with per-nuclide macroscopic total xs (MT=1) on the unified grid.
    pub fn calculate_macroscopic_xs(
        &mut self,
        mt_filter: &Vec<i32>,
        by_nuclide: bool,
    ) -> (Vec<f64>, HashMap<i32, Vec<f64>>) {
        // Ensure nuclides are loaded before proceeding
        if let Err(e) = self.ensure_nuclides_loaded() {
            panic!("Error loading nuclides: {}", e);
        }
        // Get the energy grid
        let energy_grid = self.unified_energy_grid_neutron();
        // ...existing code...

        // If by_nuclide is true, ensure MT=1 is in the filter
        if by_nuclide && !mt_filter.contains(&1) {
            panic!("If by_nuclide is true, mt_filter must contain 1 (total). Otherwise, per-nuclide total cross section makes no sense.");
        }

        // Use the filter directly, as all hierarchical MTs are now present in the JSON files.
        let micro_xs = self.calculate_microscopic_xs_neutron(Some(mt_filter));

        // Create a map to hold macroscopic cross sections for each MT (i32)
        let mut macro_xs: HashMap<i32, Vec<f64>> = HashMap::new();
        // Find all unique MT numbers across all nuclides (as i32)
        let mut all_mts = std::collections::HashSet::new();
        for nuclide_data in micro_xs.values() {
            for &mt in nuclide_data.keys() {
                if mt_filter.contains(&mt) {
                    all_mts.insert(mt);
                }
            }
        }
        // Get the grid length (from any MT reaction of any nuclide, all should have same length)
        let grid_length = micro_xs
            .values()
            .next()
            .and_then(|xs| xs.values().next())
            .map_or(0, |v| v.len());
        // Initialize macro_xs with zeros for each MT
        for &mt in &all_mts {
            macro_xs.insert(mt, vec![0.0; grid_length]);
        }
        // Calculate macroscopic cross section for each MT
        // Get atoms per barn-cm for all nuclides
        let atoms_per_bcm_map = self.get_atoms_per_barn_cm();
        // ...existing code...
        // Optionally: collect per-nuclide macroscopic total xs (MT=1) if requested
        let mut by_nuclide_map: Option<HashMap<String, Vec<f64>>> = if by_nuclide {
            Some(HashMap::new())
        } else {
            None
        };

        for (nuclide, _) in &self.nuclides {
            let atoms_per_bcm = atoms_per_bcm_map.get(nuclide);
            let nuclide_data = micro_xs.get(nuclide);
            // Always try to store per-nuclide MT=1 if by_nuclide is true
            if by_nuclide_map.is_some() {
                if let (Some(nuclide_data), Some(atoms_per_bcm)) = (nuclide_data, atoms_per_bcm) {
                    if let Some(xs_values) = nuclide_data.get(&1) {
                        let macro_vec: Vec<f64> =
                            xs_values.iter().map(|&xs| atoms_per_bcm * xs).collect();
                        by_nuclide_map
                            .as_mut()
                            .unwrap()
                            .insert(nuclide.clone(), macro_vec);
                    } else {
                        by_nuclide_map
                            .as_mut()
                            .unwrap()
                            .insert(nuclide.clone(), vec![0.0; energy_grid.len()]);
                    }
                } else {
                    by_nuclide_map
                        .as_mut()
                        .unwrap()
                        .insert(nuclide.clone(), vec![0.0; energy_grid.len()]);
                }
            }
            if let (Some(nuclide_data), Some(atoms_per_bcm)) = (nuclide_data, atoms_per_bcm) {
                for (&mt, xs_values) in nuclide_data {
                    if let Some(macro_values) = macro_xs.get_mut(&mt) {
                        for (i, &xs) in xs_values.iter().enumerate() {
                            macro_values[i] += atoms_per_bcm * xs;
                        }
                    }
                }
            }
        }

        // If by_nuclide was requested, update struct field
        if by_nuclide {
            self.macroscopic_xs_neutron_total_by_nuclide = by_nuclide_map;
        } else {
            self.macroscopic_xs_neutron_total_by_nuclide = None;
        }
        // Cache the results in the material
        self.macroscopic_xs_neutron = macro_xs.clone();
        // All hierarchical MTs are now constructed in Python and present in the JSON files.
        // No need to generate or copy hierarchical MTs here.
        (energy_grid, macro_xs)
    }

    /// Calculate macroscopic neutron cross sections with flexible reaction parameter.
    ///
    /// This method calculates the macroscopic cross section by accepting either
    /// integer MT numbers or string reaction names (like "(n,gamma)", "fission").
    ///
    /// # Arguments
    /// * `reaction` - Either an integer MT number or a string reaction name
    ///
    /// # Returns
    /// * A tuple of (energy_grid, cross_section_values)
    pub fn macroscopic_cross_section<R>(&mut self, reaction: R) -> (Vec<f64>, Vec<f64>)
    where
        R: Into<crate::nuclide::ReactionIdentifier>,
    {
        // Convert reaction identifier to MT number
        let reaction_id: crate::nuclide::ReactionIdentifier = reaction.into();
        let mt = match reaction_id {
            crate::nuclide::ReactionIdentifier::Mt(mt_num) => mt_num,
            crate::nuclide::ReactionIdentifier::Name(name) => {
                crate::data::REACTION_MT.get(name.as_str())
                    .copied()
                    .unwrap_or_else(|| panic!("Unknown reaction name '{}'", name))
            }
        };

        // Calculate macroscopic cross sections for this MT
        let mt_filter = vec![mt];
        let (energy_grid, xs_map) = self.calculate_macroscopic_xs(&mt_filter, false);
        
        // Extract the cross section for the requested MT
        let xs_values = xs_map.get(&mt)
            .unwrap_or_else(|| panic!("No cross section data found for MT {}", mt))
            .clone();

        (xs_values, energy_grid)
    }

    /// Calculate the neutron mean free path at a given energy
    ///
    /// This method calculates the mean free path of a neutron at a specific energy
    /// by interpolating the total macroscopic cross section and then taking 1/Σ.
    ///
    /// If the total macroscopic cross section hasn't been calculated yet, it will
    /// automatically call calculate_total_xs_neutron() first.
    ///
    /// # Arguments
    /// * `energy` - The energy of the neutron in eV
    ///
    /// # Returns
    /// * The mean free path in cm, or None if there's no cross section data
    pub fn mean_free_path_neutron(&mut self, energy: f64) -> Option<f64> {
        // Ensure we have a total cross section
        if !self.macroscopic_xs_neutron.contains_key(&1) {
            let mt_filter = vec![1];
            self.calculate_macroscopic_xs(&mt_filter, false);
        }
        // If we still don't have a total cross section, return None
        if !self.macroscopic_xs_neutron.contains_key(&1) {
            return None;
        }
        // Get the total cross section and energy grid
        let total_xs = &self.macroscopic_xs_neutron[&1];

        // If we have an empty cross section array, return None
        if total_xs.is_empty() || self.unified_energy_grid_neutron.is_empty() {
            return None;
        }

        // Make sure the energy grid and cross section have the same length
        if total_xs.len() != self.unified_energy_grid_neutron.len() {
            eprintln!("Error: Energy grid and cross section lengths don't match");
            return None;
        }

        // Interpolate to get the cross section at the requested energy
        // Using linear-linear interpolation
        let cross_section = interpolate_linear(&self.unified_energy_grid_neutron, total_xs, energy);

        // Mean free path = 1/Σ
        // Check for zero to avoid division by zero
        if cross_section <= 0.0 {
            None
        } else {
            Some(1.0 / cross_section)
        }
    }

    /// Calculate atoms per barn-centimeter for each nuclide in the material
    ///
    /// This method calculates the number density of atoms for each nuclide,
    /// using the atomic fractions and material density.
    ///
    /// Returns a HashMap mapping nuclide symbols to their atom density in atoms/b-cm,
    /// which is the unit used by OpenMC (atoms per barn-centimeter).
    /// Returns an empty HashMap if the material density is not set.
    pub fn get_atoms_per_barn_cm(&self) -> HashMap<String, f64> {
        let mut atoms_per_bcm = HashMap::new();

        // Return empty HashMap if density is not set
        if self.density.is_none() {
            panic!("Cannot calculate atoms per barn-cm: Material has no density defined");
        }

        // Return empty HashMap if no nuclides are defined
        if self.nuclides.is_empty() {
            panic!("Cannot calculate atoms per barn-cm: Material has no nuclides defined");
        }

        // Convert density to g/cm³ if necessary
        let mut density = self.density.unwrap();

        // Handle different density units
        match self.density_units.as_str() {
            "g/cm3" => (),                         // Already in the right units
            "kg/m3" => density = density / 1000.0, // Convert kg/m³ to g/cm³
            _ => {
                // For any other units, just use the value as is, but it may give incorrect results
            }
        }

        // Use canonical atomic masses from crate::data::ATOMIC_MASSES
        let atomic_masses = &crate::data::ATOMIC_MASSES;
        // First, get the atomic masses for all nuclides
        let mut nuclide_masses = HashMap::new();
        for (nuclide, _) in &self.nuclides {
            let mass = if let Some(mass_value) = atomic_masses.get(nuclide.as_str()) {
                *mass_value
            } else {
                panic!(
                    "Atomic mass for nuclide '{}' not found in the database",
                    nuclide
                );
            };
            nuclide_masses.insert(nuclide.clone(), mass);
        }

        // Normalize the fractions to sum to 1.0 for all cases
        let total_fraction: f64 = self.nuclides.values().sum();

        // Calculate the average molar mass (weighted)
        let mut weighted_mass_sum = 0.0;
        for (nuclide, &fraction) in &self.nuclides {
            let mass = nuclide_masses.get(nuclide).unwrap();
            weighted_mass_sum += fraction * mass;
        }
        let average_molar_mass = weighted_mass_sum / total_fraction;

        // Calculate atom densities using OpenMC's approach
        for (nuclide, &fraction) in &self.nuclides {
            let normalized_fraction = fraction / total_fraction;
            let _mass = nuclide_masses.get(nuclide).unwrap();

            // For a mixture, use the formula:
            // atom_density = density * N_A / avg_molar_mass * normalized_fraction * 1e-24
            const AVOGADRO: f64 = 6.02214076e23;
            let atom_density =
                density * AVOGADRO / average_molar_mass * normalized_fraction * 1.0e-24;

            atoms_per_bcm.insert(nuclide.clone(), atom_density);
        }

        atoms_per_bcm
    }

    pub fn add_element(&mut self, element: impl AsRef<str>, fraction: f64) -> Result<(), String> {
        if fraction <= 0.0 {
            return Err(String::from("Fraction must be positive"));
        }

        // Canonicalize input: trim only (do not lowercase or otherwise change user input)
        let input = element.as_ref().trim();

        // Try to match as symbol (case-sensitive, exact match)
        let mut found_symbol: Option<String> = None;
        for (symbol, _name) in ELEMENT_NAMES.iter() {
            if *symbol == input {
                found_symbol = Some(symbol.to_string());
                break;
            }
        }
        // If not found as symbol, try to match as name (case-sensitive, exact match)
        if found_symbol.is_none() {
            for (symbol, name) in ELEMENT_NAMES.iter() {
                if *name == input {
                    found_symbol = Some(symbol.to_string());
                    break;
                }
            }
        }
        let element_sym = match found_symbol {
            Some(sym) => sym,
            None => {
                return Err(format!(
                    "Element '{}' is not a recognized element symbol or name (case-sensitive, must match exactly)",
                    element.as_ref()
                ));
            }
        };

        // Get the isotopes for this element
        // Use the static ELEMENT_NUCLIDES map directly to avoid rebuilding a full map each call.
        let isotopes_vec = crate::data::ELEMENT_NUCLIDES
            .get(element_sym.as_str())
            .ok_or_else(|| {
                format!(
                    "Element '{}' not found in the natural abundance database",
                    element_sym
                )
            })?;

        // Add each isotope with its natural abundance
        for &isotope in isotopes_vec.iter() {
            if let Some(abundance) = crate::data::NATURAL_ABUNDANCE.get(isotope) {
                let isotope_fraction = fraction * abundance;
                if isotope_fraction > 0.0 {
                    self.add_nuclide(isotope, isotope_fraction)?;
                }
            }
        }
        Ok(())
    }

    /// Returns a sorted list of all unique MT numbers available in this material (across all nuclides).
    /// Ensures all nuclide JSON data is loaded.
    pub fn reaction_mts(&mut self) -> Result<Vec<i32>, Box<dyn std::error::Error>> {
        // Ensure all nuclides are loaded using the global config
        self.ensure_nuclides_loaded()?;
        let mut mt_set = std::collections::HashSet::new();
        for nuclide in self.nuclide_data.values() {
            if let Some(mts) = nuclide.reaction_mts() {
                for mt in mts {
                    mt_set.insert(mt);
                }
            }
        }
        let mut mt_vec: Vec<i32> = mt_set.into_iter().collect();
        mt_vec.sort();
        Ok(mt_vec)
    }

    /// Sample which nuclide a neutron interacts with at a given energy, using per-nuclide macroscopic total xs
    /// Returns the nuclide name as a String, or None if not possible
    pub fn sample_interacting_nuclide<R: rand::Rng + ?Sized>(
        &self,
        energy: f64,
        rng: &mut R,
    ) -> String {
        let by_nuclide = self.macroscopic_xs_neutron_total_by_nuclide.as_ref().expect("macroscopic_xs_neutron_total_by_nuclide is None: call calculate_macroscopic_xs with by_nuclide=true first");
        let mut xs_by_nuclide = Vec::new();
        let mut total = 0.0;
        let mut debug_info = String::new();
        for (nuclide, xs_vec) in by_nuclide.iter() {
            if xs_vec.is_empty() || self.unified_energy_grid_neutron.is_empty() {
                debug_info.push_str(&format!("{}: EMPTY\n", nuclide));
                continue;
            }
            let xs = crate::utilities::interpolate_linear(
                &self.unified_energy_grid_neutron,
                xs_vec,
                energy,
            );
            debug_info.push_str(&format!("{}: xs = {}\n", nuclide, xs));
            if xs > 0.0 {
                xs_by_nuclide.push((nuclide, xs));
                total += xs;
            }
        }
        if xs_by_nuclide.is_empty() || total <= 0.0 {
            panic!(
                "No nuclide has nonzero macroscopic total cross section at energy {}. Details:\n{}",
                energy, debug_info
            );
        }
        let xi = rng.gen_range(0.0..total);
        let mut accum = 0.0;
        for (nuclide, xs) in xs_by_nuclide {
            accum += xs;
            if xi < accum {
                return nuclide.clone();
            }
        }
        panic!("Failed to sample nuclide: numerical error in sampling loop");
    }
}


    #[cfg(test)]
    mod tests {
        use super::*;
        #[test]
        fn test_set_and_get_name() {
            let mut mat = Material::new();
            assert_eq!(mat.get_name(), None);
            mat.set_name("TestMaterial");
            assert_eq!(mat.get_name(), Some("TestMaterial"));
            mat.set_name("AnotherName");
            assert_eq!(mat.get_name(), Some("AnotherName"));
        }
        #[test]
        fn test_sample_distance_to_collision() {
            use rand::rngs::StdRng;
            use rand::SeedableRng;
            let mut material = Material::new();
            // Set up a mock total cross section and energy grid
            material.unified_energy_grid_neutron = vec![1.0, 10.0, 100.0];
            material
                .macroscopic_xs_neutron
                .insert(1, vec![2.0, 2.0, 2.0]);
            let mut rng = StdRng::seed_from_u64(42);
            let energy = 5.0;
            // Sample 200 times and check the average is close to expected mean
            let mut samples = Vec::with_capacity(200);
            for _ in 0..200 {
                let distance = material.sample_distance_to_collision(energy, &mut rng);
                assert!(distance.is_some());
                samples.push(distance.unwrap());
            }
            // For sigma_t = 2.0, mean = 1/sigma_t = 0.5
            let avg: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
            let expected_mean = 0.5;
            let tolerance = 0.05; // 10% tolerance
            assert!(
                (avg - expected_mean).abs() < tolerance,
                "Average sampled distance incorrect: got {}, expected {}",
                avg,
                expected_mean
            );
        }
    #[allow(unused_imports)]
    use super::Material;
    #[test]
    fn test_macroscopic_xs_neutron_total_by_nuclide_li6_li7() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 0.5).unwrap();
        material.add_nuclide("Li7", 0.5).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // Call with by_nuclide = true
        let mt_filter = vec![1];
        let (_grid, _macro_xs) = material.calculate_macroscopic_xs(&mt_filter, true);
        // Check that macroscopic_xs_neutron_total_by_nuclide is Some and contains both nuclides
        let by_nuclide = material
            .macroscopic_xs_neutron_total_by_nuclide
            .as_ref()
            .expect("macroscopic_xs_neutron_total_by_nuclide should be Some");
        assert!(
            by_nuclide.contains_key("Li6"),
            "macroscopic_xs_neutron_total_by_nuclide should contain Li6"
        );
        assert!(
            by_nuclide.contains_key("Li7"),
            "macroscopic_xs_neutron_total_by_nuclide should contain Li7"
        );
        // Optionally, check that the vectors are non-empty and same length as energy grid
        let grid_len = material.unified_energy_grid_neutron.len();
        assert_eq!(
            by_nuclide["Li6"].len(),
            grid_len,
            "Li6 xs vector should match grid length"
        );
        assert_eq!(
            by_nuclide["Li7"].len(),
            grid_len,
            "Li7 xs vector should match grid length"
        );
    }

    #[test]
    fn test_macroscopic_xs_mt3_does_not_generate_mt1() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap();
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        let mt_filter = vec![3];
        let (_grid, macro_xs) = material.calculate_macroscopic_xs(&mt_filter, false);
        assert!(
            !macro_xs.contains_key(&1),
            "MT=1 should NOT be present when only MT=3 is requested"
        );
    }

    #[test]
    fn test_macroscopic_xs_mt24_does_not_generate_mt1() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap();
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        let mt_filter = vec![24];
        let (_grid, macro_xs) = material.calculate_macroscopic_xs(&mt_filter, false);
        assert!(
            !macro_xs.contains_key(&1),
            "MT=1 should NOT be present when only MT=24 is requested"
        );
    }
    #[test]
    fn test_hierarchical_mt3_generated_for_li6() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap();
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // Request only MT=3 (now present in JSON)
        let mt_filter = vec![3];
        let (_grid, macro_xs) = material.calculate_macroscopic_xs(&mt_filter, false);
        assert!(
            macro_xs.contains_key(&3),
            "MT=3 should be present in macro_xs for Li6"
        );
        // Just check that the cross section is nonzero and matches the JSON
        let xs = &macro_xs[&3];
        assert!(
            xs.iter().any(|&v| v > 0.0),
            "MT=3 cross section should have nonzero values"
        );
    }
    // Removed unused `use super::*;` (all required items referenced explicitly)

    #[test]
    fn test_new_material() {
        let material = Material::new();
        assert!(material.nuclides.is_empty());
        assert_eq!(material.density, None);
        assert_eq!(material.density_units, "g/cm3");
    }

    #[test]
    fn test_add_nuclide() {
        let mut material = Material::new();

        // Test adding a valid nuclide
        let result = material.add_nuclide("U235", 0.05);
        assert!(result.is_ok());
        assert_eq!(material.nuclides.get("U235"), Some(&0.05));

        // Test adding another nuclide
        let result = material.add_nuclide("U238", 0.95);
        assert!(result.is_ok());
        assert_eq!(material.nuclides.get("U238"), Some(&0.95));
        assert_eq!(material.nuclides.len(), 2);

        // Test overwriting an existing nuclide
        let result = material.add_nuclide("U235", 0.1);
        assert!(result.is_ok());
        assert_eq!(material.nuclides.get("U235"), Some(&0.1));
    }

    #[test]
    fn test_add_nuclide_negative_fraction() {
        let mut material = Material::new();

        // Test adding a nuclide with negative fraction
        let result = material.add_nuclide("U235", -0.05);
        assert!(result.is_err());
        assert_eq!(result.unwrap_err(), "Fraction cannot be negative");
        assert!(material.nuclides.is_empty());
    }

    #[test]
    fn test_set_density() {
        let mut material = Material::new();

        // Test setting a valid density
        let result = material.set_density("g/cm3", 10.5);
        assert!(result.is_ok());
        assert_eq!(material.density, Some(10.5));
        assert_eq!(material.density_units, "g/cm3");

        // Test setting a different unit
        let result = material.set_density("kg/m3", 10500.0);
        assert!(result.is_ok());
        assert_eq!(material.density, Some(10500.0));
        assert_eq!(material.density_units, "kg/m3");
    }

    #[test]
    fn test_set_density_negative_value() {
        let mut material = Material::new();

        // Test setting a negative density
        let result = material.set_density("g/cm3", -10.5);
        assert!(result.is_err());
        assert_eq!(result.unwrap_err(), "Density must be positive");
        assert_eq!(material.density, None);
    }

    #[test]
    fn test_set_density_zero_value() {
        let mut material = Material::new();

        // Test setting a zero density
        let result = material.set_density("g/cm3", 0.0);
        assert!(result.is_err());
        assert_eq!(result.unwrap_err(), "Density must be positive");
        assert_eq!(material.density, None);
    }

    #[test]
    fn test_material_clone() {
        let mut material = Material::new();
        material.add_nuclide("U235", 0.05).unwrap();
        material.add_nuclide("U238", 0.95).unwrap();
        material.set_density("g/cm3", 19.1).unwrap();

        let cloned = material.clone();

        assert_eq!(cloned.nuclides.get("U235"), Some(&0.05));
        assert_eq!(cloned.nuclides.get("U238"), Some(&0.95));
        assert_eq!(cloned.density, Some(19.1));
        assert_eq!(cloned.density_units, "g/cm3");
    }

    #[test]
    fn test_material_debug() {
        let mut material = Material::new();
        material.add_nuclide("U235", 0.05).unwrap();
        material.set_density("g/cm3", 19.1).unwrap();

        // This test merely ensures that the Debug implementation doesn't panic
        let _debug_str = format!("{:?}", material);
        assert!(true);
    }

    #[test]
    fn test_volume_get_and_set() {
        let mut material = Material::new();

        // Test setting a valid volume
        let result = material.volume(Some(100.0));
        assert!(result.is_ok());
        assert_eq!(material.volume, Some(100.0));

        // Test getting the current volume
        let current_volume = material.volume(None).unwrap();
        assert_eq!(current_volume, Some(100.0));

        // Test setting an invalid (negative) volume
        let result = material.volume(Some(-50.0));
        assert!(result.is_err());
        assert_eq!(result.unwrap_err(), "Volume must be positive");
        assert_eq!(material.volume, Some(100.0)); // Ensure the volume wasn't changed
    }

    #[test]
    fn test_get_nuclides() {
        let mut material = Material::new();

        // Empty material should return empty vector
        assert!(material.get_nuclides().is_empty());

        // Add some nuclides
        material.add_nuclide("U235", 0.05).unwrap();
        material.add_nuclide("U238", 0.95).unwrap();
        material.add_nuclide("O16", 2.0).unwrap();

        // Check the result is sorted
        let nuclides = material.get_nuclides();
        assert_eq!(
            nuclides,
            vec!["O16".to_string(), "U235".to_string(), "U238".to_string()]
        );
    }

    #[test]
    fn test_get_atoms_per_barn_cm() {
        let material = Material::new();

        // Test with no density set - should panic
        let result = std::panic::catch_unwind(|| material.get_atoms_per_barn_cm());
        assert!(result.is_err(), "Should panic when density is not set");

        // Test with single nuclide case
        let mut material_single = Material::new();
        material_single.add_nuclide("Li6", 2.5).unwrap();
        material_single.set_density("g/cm3", 1.0).unwrap();

        let atoms_single = material_single.get_atoms_per_barn_cm();
        assert_eq!(
            atoms_single.len(),
            1,
            "Should have 1 nuclide in the HashMap"
        );

        // For a single nuclide, we normalize the fraction to 1.0
        let avogadro = 6.02214076e23;
        let li6_mass = 6.01512288742;
        let li6_expected = 1.0 * avogadro / li6_mass * 1.0e-24; // Fraction is normalized to 1.0

        let li6_actual = atoms_single.get("Li6").unwrap();
        let tolerance = 0.01; // 1%

        assert!(
            (li6_actual - li6_expected).abs() / li6_expected < tolerance,
            "Li6 atoms/cc calculation (single nuclide) is incorrect: got {}, expected {}",
            li6_actual,
            li6_expected
        );

        // Test with multiple nuclides
        let mut material_multi = Material::new();
        material_multi.add_nuclide("Li6", 0.5).unwrap();
        material_multi.add_nuclide("Li7", 0.5).unwrap();
        material_multi.set_density("g/cm3", 1.0).unwrap();

        let atoms_multi = material_multi.get_atoms_per_barn_cm();
        assert_eq!(
            atoms_multi.len(),
            2,
            "Should have 2 nuclides in the HashMap"
        );

        // For multiple nuclides, the fractions are normalized and used with average molar mass
        let li7_mass = 7.016004;
        let avg_mass = (0.5 * li6_mass + 0.5 * li7_mass) / 1.0; // weighted average

        let li6_expected_multi = 1.0 * avogadro / avg_mass * (0.5 / 1.0) * 1.0e-24;
        let li7_expected_multi = 1.0 * avogadro / avg_mass * (0.5 / 1.0) * 1.0e-24;

        let li6_actual_multi = atoms_multi.get("Li6").unwrap();
        let li7_actual_multi = atoms_multi.get("Li7").unwrap();

        assert!(
            (li6_actual_multi - li6_expected_multi).abs() / li6_expected_multi < tolerance,
            "Li6 atoms/cc calculation (multiple nuclides) is incorrect: got {}, expected {}",
            li6_actual_multi,
            li6_expected_multi
        );
        assert!(
            (li7_actual_multi - li7_expected_multi).abs() / li7_expected_multi < tolerance,
            "Li7 atoms/cc calculation (multiple nuclides) is incorrect: got {}, expected {}",
            li7_actual_multi,
            li7_expected_multi
        );

        // Test with non-normalized fractions
        let mut material_non_norm = Material::new();
        material_non_norm.add_nuclide("Li6", 1.0).unwrap();
        material_non_norm.add_nuclide("Li7", 1.0).unwrap(); // Total fractions = 2.0
        material_non_norm.set_density("g/cm3", 1.0).unwrap();

        let atoms_non_norm = material_non_norm.get_atoms_per_barn_cm();

        // Fractions should be normalized to 0.5 each (1.0/2.0)
        let avg_mass_non_norm = (1.0 * li6_mass + 1.0 * li7_mass) / 2.0;
        let li6_expected_non_norm = 1.0 * avogadro / avg_mass_non_norm * (1.0 / 2.0) * 1.0e-24;
        let li7_expected_non_norm = 1.0 * avogadro / avg_mass_non_norm * (1.0 / 2.0) * 1.0e-24;

        let li6_actual_non_norm = atoms_non_norm.get("Li6").unwrap();
        let li7_actual_non_norm = atoms_non_norm.get("Li7").unwrap();

        assert!(
            (li6_actual_non_norm - li6_expected_non_norm).abs() / li6_expected_non_norm < tolerance,
            "Li6 normalized atoms/cc calculation is incorrect: got {}, expected {}",
            li6_actual_non_norm,
            li6_expected_non_norm
        );
        assert!(
            (li7_actual_non_norm - li7_expected_non_norm).abs() / li7_expected_non_norm < tolerance,
            "Li7 normalized atoms/cc calculation is incorrect: got {}, expected {}",
            li7_actual_non_norm,
            li7_expected_non_norm
        );
    }

    #[test]
    fn test_get_atoms_per_barn_cm_no_density() {
        let material = Material::new();

        // Should panic when density is not set
        let result = std::panic::catch_unwind(|| material.get_atoms_per_barn_cm());
        assert!(result.is_err(), "Should panic when density is not set");

        // Should also panic when nuclides are not added
        let mut material_with_density = Material::new();
        material_with_density.set_density("g/cm3", 1.0).unwrap();

        let result = std::panic::catch_unwind(|| material_with_density.get_atoms_per_barn_cm());
        assert!(result.is_err(), "Should panic when no nuclides are defined");
    }

    #[test]
    fn test_mean_free_path_neutron() {
        // Create a properly set up material
        let mut material = Material::new();

        // Test with empty material - should return None
        // We can't use catch_unwind easily with &mut material, so we'll bypass the part that panics
        // by directly setting up the test with mock data

        // Skip the initial test with empty material that would cause a panic

        // Add a nuclide and set density
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();

        // Create mock cross sections directly (bypassing the normal calculation path)
        let energy_grid = vec![
            1.0,
            10.0,
            100.0,
            1000.0,
            10000.0,
            100000.0,
            1000000.0,
            10000000.0,
            100000000.0,
        ];
        material.unified_energy_grid_neutron = energy_grid.clone();

        // Set total cross section (in barns * atoms/cm³, which gives cm⁻¹)
        // Intentionally using a simple pattern that's easy to verify
        let total_xs = vec![
            1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625,
        ];
        material.macroscopic_xs_neutron.insert(1, total_xs.clone());

        // Test exact values from our mock data
        assert_eq!(material.mean_free_path_neutron(1.0), Some(1.0));
        assert_eq!(material.mean_free_path_neutron(10.0), Some(2.0));
        assert_eq!(material.mean_free_path_neutron(100.0), Some(4.0));

        // Test interpolated value
        // At energy = 3.0, we're using linear interpolation between 1.0 and 10.0
        // Cross section should be about 0.889 (linearly interpolated between 1.0 and 0.5)
        // Mean free path should be about 1.125
        let mfp_3ev = material.mean_free_path_neutron(3.0).unwrap();
        assert!(
            (mfp_3ev - 1.125).abs() < 0.01,
            "Expected ~1.125, got {}",
            mfp_3ev
        );

        // Test outside of range (should use endpoint value)
        assert_eq!(material.mean_free_path_neutron(0.1), Some(1.0)); // Below range
        assert_eq!(material.mean_free_path_neutron(1e9), Some(256.0)); // Above range
    }

    #[test]
    fn test_mean_free_path_lithium_14mev() {
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap(); // lithium density
                                                       // Set up mock cross section data for 14 MeV (1.4e7 eV)
                                                       // We'll use a simple grid and cross section for demonstration
        let energy_grid = vec![1e6, 1.4e7, 1e8]; // eV
        let total_xs = vec![1.0, 0.5, 0.2]; // barns * atoms/cm³, so cm⁻¹
        material.unified_energy_grid_neutron = energy_grid.clone();
        material.macroscopic_xs_neutron.insert(1, total_xs.clone());
        // 14 MeV = 1.4e7 eV
        let mfp = material.mean_free_path_neutron(1.4e7);
        assert!(mfp.is_some());
        let mfp_val = mfp.unwrap();
        // At 14 MeV, total_xs = 0.5, so mean free path = 1/0.5 = 2.0 cm
        assert!(
            (mfp_val - 2.0).abs() < 1e-6,
            "Expected 2.0 cm, got {}",
            mfp_val
        );
    }

    #[test]
    fn test_add_element() {
        let mut material = Material::new();
        // Test adding natural lithium
        let result = material.add_element("Li", 1.0);
        assert!(result.is_ok());
        // Verify the isotopes were added correctly
        assert!(material.nuclides.contains_key("Li6"));
        assert!(material.nuclides.contains_key("Li7"));
        // Check the fractions are correct
        assert_eq!(*material.nuclides.get("Li6").unwrap(), 0.07589);
        assert_eq!(*material.nuclides.get("Li7").unwrap(), 0.92411);
        // Test adding an element with many isotopes
        let mut material2 = Material::new();
        let result = material2.add_element("Sn", 1.0); // Tin has 10 isotopes
        assert!(result.is_ok());
        assert_eq!(material2.nuclides.len(), 10);
    }

    #[test]
    fn test_add_element_invalid() {
        let mut material = Material::new();
        // Test with negative fraction
        let result = material.add_element("Li", -1.0);
        assert!(result.is_err());
        // Test with invalid element
        let result = material.add_element("Xx", 1.0);
        assert!(result.is_err());
    }

    #[test]
    fn test_add_element_by_symbol_and_name() {
        let mut material = Material::new();
        // By symbol (case-sensitive, exact match)
        assert!(material.add_element("Li", 1.0).is_ok());
        assert!(material.nuclides.contains_key("Li6"));
        assert!(material.nuclides.contains_key("Li7"));
        // By full name (case-sensitive, exact match)
        let mut material2 = Material::new();
        assert!(material2.add_element("gold", 1.0).is_ok());
        assert!(material2.nuclides.contains_key("Au197"));
        // By full name (lowercase) - should fail
        let mut material3 = Material::new();
        assert!(material3.add_element("Lithium", 1.0).is_err());
        // By symbol (lowercase) - should fail
        let mut material4 = Material::new();
        assert!(material4.add_element("li", 1.0).is_err());
        // By symbol (uppercase) - should fail
        let mut material5 = Material::new();
        assert!(material5.add_element("LI", 1.0).is_err());
        // Invalid name
        let mut material6 = Material::new();
        let result = material6.add_element("notanelement", 1.0);
        assert!(result.is_err());
    }

    #[test]
    fn test_add_element_beryllium_and_iron() {
        let mut mat_be = Material::new();
        assert!(mat_be.add_element("Be", 1.0).is_ok());
        // Beryllium has only one stable isotope
        assert_eq!(mat_be.nuclides.len(), 1);
        assert!(mat_be.nuclides.contains_key("Be9"));
        // Check the fraction is 1.0 for Be9
        assert_eq!(*mat_be.nuclides.get("Be9").unwrap(), 1.0);

        let mut mat_fe = Material::new();
        assert!(mat_fe.add_element("Fe", 1.0).is_ok());
        // Iron has four stable isotopes
        assert!(mat_fe.nuclides.contains_key("Fe54"));
        assert!(mat_fe.nuclides.contains_key("Fe56"));
        assert!(mat_fe.nuclides.contains_key("Fe57"));
        assert!(mat_fe.nuclides.contains_key("Fe58"));
        // Check that the sum of fractions is 1.0 (within tolerance)
        let sum: f64 = mat_fe.nuclides.values().sum();
        assert!((sum - 1.0).abs() < 1e-6);
    }
    #[test]
    fn test_material_reaction_mts_lithium() {
        use crate::material::Material;
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        // Prepare the nuclide JSON map for Li6 and Li7
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        // Read in the nuclear data
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // This will load Li6 and Li7, so the MTs should be the union of both, including hierarchical MTs
        let mts = material.reaction_mts().expect("Failed to get reaction MTs");
        // The expected list should match the actual list from the JSON, including hierarchical MTs
        let expected = vec![
            1, 2, 3, 4, 5, 16, 24, 25, 27, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
            65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 101, 102, 103,
            104, 105, 203, 204, 205, 206, 207, 301, 444,
        ];
        assert_eq!(
            mts, expected,
            "Material lithium MT list does not match expected. Got {:?}",
            mts
        );
    }

    #[test]
    fn test_macroscopic_cross_section_flexible_reactions() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 0.5).unwrap();
        material.add_nuclide("Li7", 0.5).unwrap();
        material.set_density("g/cm3", 2.0).unwrap();
        
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");

        // Test with integer MT number
        let (xs1, energy1) = material.macroscopic_cross_section(1);
        assert!(!energy1.is_empty(), "Energy grid should not be empty");
        assert!(!xs1.is_empty(), "Cross section should not be empty");
        assert_eq!(energy1.len(), xs1.len(), "Energy and cross section arrays should have same length");

        // Test with string reaction name - same reaction
        let (xs2, energy2) = material.macroscopic_cross_section("(n,total)".to_string());
        assert_eq!(energy1.len(), energy2.len(), "Energy grids should be same length");
        assert_eq!(xs1.len(), xs2.len(), "Cross sections should be same length");
        
        // Values should be identical (or very close due to floating point)
        for (i, (&val1, &val2)) in xs1.iter().zip(xs2.iter()).enumerate() {
            assert!(
                (val1 - val2).abs() < 1e-10,
                "Cross section values should be identical at index {}: {} vs {}",
                i, val1, val2
            );
        }

        // Test with gamma capture reaction
        let (xs3, energy3) = material.macroscopic_cross_section("(n,gamma)");
        assert!(!energy3.is_empty(), "Energy grid should not be empty for (n,gamma)");
        assert!(!xs3.is_empty(), "Cross section should not be empty for (n,gamma)");
        assert_eq!(energy3.len(), xs3.len(), "Energy and cross section arrays should have same length for (n,gamma)");
    }

    #[test]
    fn test_selective_temperature_load_be9_300() {
        crate::nuclide::clear_nuclide_cache();
        let mut mat = Material::new();
        mat.add_nuclide("Be9", 1.0).unwrap();
        mat.set_temperature("300");
        let mut map = std::collections::HashMap::new();
        map.insert("Be9".to_string(), "tests/Be9.json".to_string());
        mat.read_nuclides_from_json(&map).unwrap();
        let be9 = mat.nuclide_data.get("Be9").expect("Be9 not loaded");
        assert_eq!(
            be9.available_temperatures,
            vec!["294".to_string(), "300".to_string()]
        );
        assert_eq!(
            be9.loaded_temperatures,
            vec!["300".to_string()],
            "Should only load 300K data"
        );
    }

    #[test]
    fn test_selective_temperature_load_be9_294() {
        crate::nuclide::clear_nuclide_cache();
        let mut mat = Material::new();
        mat.add_nuclide("Be9", 1.0).unwrap();
        mat.set_temperature("294");
        let mut map = std::collections::HashMap::new();
        map.insert("Be9".to_string(), "tests/Be9.json".to_string());
        mat.read_nuclides_from_json(&map).unwrap();
        let be9 = mat.nuclide_data.get("Be9").expect("Be9 not loaded");
        assert_eq!(
            be9.available_temperatures,
            vec!["294".to_string(), "300".to_string()]
        );
        assert_eq!(
            be9.loaded_temperatures,
            vec!["294".to_string()],
            "Should only load 294K data"
        );
    }

    #[test]
    fn test_calculate_microscopic_xs_neutron_lithium() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        // Prepare the nuclide JSON map for Li6 and Li7
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        // Read in the nuclear data
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // Build the unified energy grid
        let grid = material.unified_energy_grid_neutron();
        // Calculate microscopic cross sections
        let micro_xs = material.calculate_microscopic_xs_neutron(None);
        // Check that both Li6 and Li7 are present
        assert!(micro_xs.contains_key("Li6"));
        assert!(micro_xs.contains_key("Li7"));
        // Check that for a known MT (e.g., 2), both nuclides have cross section data
        let mt = 2;
        assert!(micro_xs["Li6"].contains_key(&mt), "Li6 missing MT=2");
        assert!(micro_xs["Li7"].contains_key(&mt), "Li7 missing MT=2");
        // Check that the cross section arrays are the same length as the grid
        assert_eq!(micro_xs["Li6"][&mt].len(), grid.len());
        assert_eq!(micro_xs["Li7"][&mt].len(), grid.len());
    }

    #[test]
    fn test_material_vs_nuclide_microscopic_xs_li6() {
        use crate::nuclide::get_or_load_nuclide;
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        // Prepare the nuclide JSON map for Li6
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        // Read in the nuclear data
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // Build the unified energy grid
        let grid = material.unified_energy_grid_neutron();
        // Calculate microscopic cross sections for the material
        let micro_xs_mat = material.calculate_microscopic_xs_neutron(None);
        // Get the nuclide directly
        let nuclide =
            get_or_load_nuclide("Li6", &nuclide_json_map, None).expect("Failed to load Li6");
        let temperature = &material.temperature;
        // Get reactions and energy grid for nuclide
        let reactions = nuclide
            .reactions
            .get(temperature)
            .expect("No reactions for Li6");
        let energy_map = nuclide.energy.as_ref().expect("No energy map for Li6");
        let energy_grid = energy_map.get(temperature).expect("No energy grid for Li6");
        // For each MT in the material, compare the cross sections
        for (mt, xs_mat) in micro_xs_mat["Li6"].iter() {
            // Only compare if MT exists in nuclide
            if let Some(reaction) = reactions.get(&mt) {
                let threshold_idx = reaction.threshold_idx;
                let nuclide_energy = if threshold_idx < energy_grid.len() {
                    &energy_grid[threshold_idx..]
                } else {
                    continue;
                };
                let xs_nuclide = &reaction.cross_section;
                // Interpolate nuclide xs onto the material grid
                let mut xs_nuclide_interp = Vec::with_capacity(grid.len());
                for &g in &grid {
                    if g < nuclide_energy[0] {
                        xs_nuclide_interp.push(0.0);
                    } else {
                        let xs =
                            crate::utilities::interpolate_linear(nuclide_energy, xs_nuclide, g);
                        xs_nuclide_interp.push(xs);
                    }
                }
                // Compare arrays (allow small tolerance)
                let tol = 1e-10;
                for (a, b) in xs_mat.iter().zip(xs_nuclide_interp.iter()) {
                    assert!(
                        (a - b).abs() < tol,
                        "Mismatch for MT {}: {} vs {}",
                        mt,
                        a,
                        b
                    );
                }
            }
        }
    }

    #[test]
    fn test_calculate_microscopic_xs_neutron_mt_filter() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        // Prepare the nuclide JSON map for Li6 and Li7
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        // Read in the nuclear data
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        // Build the unified energy grid
        let _grid = material.unified_energy_grid_neutron();
        // Calculate microscopic cross sections for all MTs
        let micro_xs_all = material.calculate_microscopic_xs_neutron(None);
        // Calculate microscopic cross sections for only MT=2
        let mt_filter = vec![2];
        let micro_xs_mt2 = material.calculate_microscopic_xs_neutron(Some(&mt_filter));
        // For each nuclide, only MT=2 should be present
        for nuclide in &["Li6", "Li7"] {
            assert!(
                micro_xs_mt2.contains_key(*nuclide),
                "{} missing in filtered result",
                nuclide
            );
            let xs_map = &micro_xs_mt2[*nuclide];
            // Assert that only the requested MT is present
            assert!(
                xs_map.keys().all(|k| *k == 2),
                "Filtered result for {} contains non-filtered MTs: {:?}",
                nuclide,
                xs_map.keys()
            );
            assert_eq!(
                xs_map.len(),
                1,
                "Filtered result for {} should have only one MT",
                nuclide
            );
            // The cross section array for MT=2 should match the unfiltered result
            let xs_all = &micro_xs_all[*nuclide][&2];
            let xs_filtered = &xs_map[&2];
            assert_eq!(
                xs_all, xs_filtered,
                "Filtered and unfiltered MT=2 xs do not match for {}",
                nuclide
            );
        }
    }

    #[test]
    fn test_calculate_macroscopic_xs_mt_filter() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap(); // lithium density
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        let mt_filter = vec![
            2, 16, 24, 25, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
            69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 101, 102, 103, 104, 105, 203,
            204, 205, 206, 207, 301, 444,
        ];
        let (_grid, macro_xs) = material.calculate_macroscopic_xs(&mt_filter, false);
        for mt in &mt_filter {
            match macro_xs.get(mt) {
                Some(xs) => assert_eq!(xs.len(), material.unified_energy_grid_neutron.len()),
                None => println!("MT {} not present in macro_xs", mt),
            }
        }
    }

    #[test]
    fn test_panic_if_no_density_for_macroscopic_xs_and_mean_free_path() {
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();

        // Prepare the nuclide JSON map for Li6
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        let _grid = material.unified_energy_grid_neutron();

        // Should panic when calculating macroscopic cross sections with no density
        let result = std::panic::catch_unwind(move || {
            let mut material = material;
            material.calculate_macroscopic_xs(&vec![1], false);
        });
        assert!(
            result.is_err(),
            "Should panic if density is not set for calculate_macroscopic_xs"
        );

        // Re-create material for the next test
        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");
        material.unified_energy_grid_neutron();

        // Should panic when calculating mean free path with no density
        let result = std::panic::catch_unwind(move || {
            let mut material = material;
            material.mean_free_path_neutron(1e6);
        });
        assert!(
            result.is_err(),
            "Should panic if density is not set for mean_free_path_neutron"
        );
    }

    #[test]
    fn test_mean_free_path_lithium_real_data() {
        use std::collections::HashMap;
        let mut material = Material::new();
        material.add_element("Li", 1.0).unwrap();
        material.set_density("g/cm3", 0.534).unwrap(); // lithium density

        // Prepare the nuclide JSON map for Li6 and Li7
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());

        // Read in the nuclear data
        material
            .read_nuclides_from_json(&nuclide_json_map)
            .expect("Failed to read nuclide JSON");

        // Calculate the mean free path at 14 MeV (1.4e7 eV)
        let mfp = material.mean_free_path_neutron(1.4e7);

        assert!(mfp.is_some(), "Mean free path should be Some for real data");
        let mfp_val = mfp.unwrap();
        // Print the value for inspection
        println!("Mean free path for lithium at 14 MeV: {} cm", mfp_val);
        // Check that the value is positive
        assert!(mfp_val > 0.0, "Mean free path should be positive");
        let expected = 14.963768069986559;
        let rel_tol = 1e-5;
        assert!(
            (mfp_val - expected).abs() / expected < rel_tol,
            "Expected ~{:.8} cm, got {:.8} cm",
            expected,
            mfp_val
        );
    }

    #[test]
    fn test_sample_distance_to_collision_li6() {
        use std::collections::HashMap;
        // Create Li6 material
        let mut mat = Material::new();
        mat.add_nuclide("Li6", 1.0).unwrap();
        // Load nuclide data from JSON
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        mat.read_nuclides_from_json(&nuclide_json_map).unwrap();
        mat.set_density("g/cm3", 1.0).unwrap();
        mat.set_temperature("294");

        // Check that the total cross section is present and nonzero at 14 MeV
        mat.calculate_macroscopic_xs(&vec![1], false);

        // Sample 1000 distances
        let mut sum = 0.0;
        let n_samples = 1000;
        use rand::rngs::StdRng;
        use rand::SeedableRng; // Needed for seed_from_u64
        for seed in 0..n_samples {
            let mut rng = StdRng::seed_from_u64(seed as u64);
            let dist = mat
                .sample_distance_to_collision(14_000_000.0, &mut rng)
                .unwrap_or_else(|| panic!("sample_distance_to_collision returned None at 14 MeV!"));
            sum += dist;
        }
        let avg = sum / n_samples as f64;
        println!("Average distance: {}", avg);
        assert!(
            (avg - 6.9).abs() < 0.1,
            "Average {} not within 0.2 of 6.9",
            avg
        );
    }

    #[test]
    fn test_sample_interacting_nuclide_li6_li7() {
        use rand::rngs::StdRng;
        use rand::SeedableRng;
        use std::collections::HashMap;

        let mut material = Material::new();
        material.add_nuclide("Li6", 0.1).unwrap();
        material.add_nuclide("Li7", 0.9).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclide data from JSON
        let mut nuclide_json_map = HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Calculate total xs to ensure everything is set up
        // For this test, calculate per-nuclide macroscopic total xs as well
        material.calculate_macroscopic_xs(&vec![1], true);

        // Check that macroscopic_xs_neutron_total_by_nuclide is present and not empty
        let by_nuclide = material.macroscopic_xs_neutron_total_by_nuclide.as_ref();
        assert!(
            by_nuclide.is_some(),
            "macroscopic_xs_neutron_total_by_nuclide should be Some after calculation"
        );
        let by_nuclide = by_nuclide.unwrap();
        assert!(
            !by_nuclide.is_empty(),
            "macroscopic_xs_neutron_total_by_nuclide should not be empty"
        );

        // Sample the interacting nuclide many times at 14 MeV
        let energy = 14_000_000.0;
        let n_samples = 10000;
        let mut counts = HashMap::new();
        for seed in 0..n_samples {
            let mut rng = StdRng::seed_from_u64(seed as u64);
            let nuclide = material.sample_interacting_nuclide(energy, &mut rng);
            *counts.entry(nuclide).or_insert(0) += 1;
        }

        let count_li6 = *counts.get("Li6").unwrap_or(&0) as f64;
        let count_li7 = *counts.get("Li7").unwrap_or(&0) as f64;
        let total = count_li6 + count_li7;
        let frac_li6 = count_li6 / total;
        let frac_li7 = count_li7 / total;

        println!("Li6 fraction: {}, Li7 fraction: {}", frac_li6, frac_li7);

        // The sampled fractions should be close to the expected probability
        // (proportional to macroscopic total xs for each nuclide at this energy)
        // For a rough test, just check both are nonzero and sum to 1
        assert!(
            frac_li6 > 0.0 && frac_li7 > 0.0,
            "Both nuclides should be sampled"
        );
        assert!(
            (frac_li6 + frac_li7 - 1.0).abs() < 1e-6,
            "Fractions should sum to 1"
        );

        // Optionally, check that Li7 is sampled much more often than Li6 (since its fraction is higher)
        assert!(
            frac_li7 > frac_li6,
            "Li7 should be sampled more often than Li6"
        );
    }

    #[test]
    fn test_cache_invalidation_add_nuclide() {
        crate::nuclide::clear_nuclide_cache();

        let mut material = Material::new();
        material.add_nuclide("Li6", 0.5).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclear data
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        nuclide_json_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Pre-populate cache by calculating cross sections
        material.calculate_macroscopic_xs(&vec![1], false);
        assert!(!material.macroscopic_xs_neutron.is_empty(), "Cache should be populated");

        // Adding a nuclide should clear the cache
        material.add_nuclide("Li7", 0.5).unwrap();
        assert!(material.macroscopic_xs_neutron.is_empty(), "Cache should be cleared after adding nuclide");
    }

    #[test]
    fn test_cache_invalidation_set_density() {
        crate::nuclide::clear_nuclide_cache();

        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclear data
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Pre-populate cache
        material.calculate_macroscopic_xs(&vec![1], false);
        assert!(!material.macroscopic_xs_neutron.is_empty(), "Cache should be populated");

        // Changing density should clear the cache
        material.set_density("g/cm3", 2.0).unwrap();
        assert!(material.macroscopic_xs_neutron.is_empty(), "Cache should be cleared after density change");
    }

    #[test]
    fn test_cache_invalidation_set_temperature() {
        crate::nuclide::clear_nuclide_cache();

        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclear data
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Pre-populate cache
        material.calculate_macroscopic_xs(&vec![1], false);
        assert!(!material.macroscopic_xs_neutron.is_empty(), "Cache should be populated");

        // Changing temperature should clear the cache
        material.set_temperature("300");
        assert!(material.macroscopic_xs_neutron.is_empty(), "Cache should be cleared after temperature change");
    }

    #[test]
    fn test_cache_invalidation_read_nuclides() {
        crate::nuclide::clear_nuclide_cache();

        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclear data first time
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Pre-populate cache
        material.calculate_macroscopic_xs(&vec![1], false);
        assert!(!material.macroscopic_xs_neutron.is_empty(), "Cache should be populated");

        // Loading nuclear data again should clear the cache (use same map since empty map would fail)
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();
        assert!(material.macroscopic_xs_neutron.is_empty(), "Cache should be cleared after loading nuclear data");
    }

    #[test]
    fn test_cache_behavior_after_calculation() {
        crate::nuclide::clear_nuclide_cache();

        let mut material = Material::new();
        material.add_nuclide("Li6", 1.0).unwrap();
        material.set_density("g/cm3", 1.0).unwrap();
        material.set_temperature("294");

        // Load nuclear data
        let mut nuclide_json_map = std::collections::HashMap::new();
        nuclide_json_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        material.read_nuclides_from_json(&nuclide_json_map).unwrap();

        // Calculate cross sections multiple times - each call replaces the cache
        material.calculate_macroscopic_xs(&vec![1], false);
        let first_size = material.macroscopic_xs_neutron.len();
        assert_eq!(first_size, 1, "Cache should contain only MT 1");

        material.calculate_macroscopic_xs(&vec![1, 2], false);
        let second_size = material.macroscopic_xs_neutron.len();
        assert_eq!(second_size, 2, "Cache should contain MT 1 and 2");

        // Verify that calling with just MT 1 replaces cache with only MT 1
        material.calculate_macroscopic_xs(&vec![1], false);
        assert_eq!(material.macroscopic_xs_neutron.len(), 1, "Cache should contain only MT 1 again");
        
        // Verify that calling with all MTs gives us all MTs in cache
        material.calculate_macroscopic_xs(&vec![1, 2], false);
        assert_eq!(material.macroscopic_xs_neutron.len(), 2, "Cache should contain MT 1 and 2 again");
    }

    #[test]
    fn test_material_different_data_sources() {
        // Test that material loading respects different data sources
        
        crate::nuclide::clear_nuclide_cache();

        // Material 1: Li6 from file
        let mut mat_li6 = Material::new();
        mat_li6.add_nuclide("Li6", 1.0).unwrap();
        let mut li6_map = std::collections::HashMap::new();
        li6_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        mat_li6.read_nuclides_from_json(&li6_map).unwrap();
        mat_li6.set_density("g/cm3", 0.534).unwrap();

        // Material 2: Li7 from file
        let mut mat_li7 = Material::new();
        mat_li7.add_nuclide("Li7", 1.0).unwrap();
        let mut li7_map = std::collections::HashMap::new();
        li7_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        mat_li7.read_nuclides_from_json(&li7_map).unwrap();
        mat_li7.set_density("g/cm3", 0.534).unwrap();

        // Get macroscopic cross sections
        let (xs_li6, _) = mat_li6.macroscopic_cross_section("(n,gamma)");
        let (xs_li7, _) = mat_li7.macroscopic_cross_section("(n,gamma)");

        // Should have different data (different nuclides)
        let data_different = xs_li6.len() != xs_li7.len() || 
                            xs_li6.iter().zip(&xs_li7).any(|(a, b)| (a - b).abs() > 1e-10);
        
        assert!(data_different, "Li6 and Li7 materials should have different cross sections");
        println!("Material Li6: {} points, Li7: {} points", xs_li6.len(), xs_li7.len());
    }

    #[test]
    fn test_material_file_and_keyword_sources() {
        // Test that materials can use both file paths and keywords
        
        crate::nuclide::clear_nuclide_cache();

        // Material 1: Li6 from file
        let mut mat_file = Material::new();
        mat_file.add_nuclide("Li6", 1.0).unwrap();
        let mut nuclide_map = std::collections::HashMap::new();
        nuclide_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        mat_file.read_nuclides_from_json(&nuclide_map).unwrap();
        mat_file.set_density("g/cm3", 1.0).unwrap();

        // Material 2: Li7 from file (different nuclide, different file)
        let mut mat_other = Material::new();
        mat_other.add_nuclide("Li7", 1.0).unwrap();
        let mut other_map = std::collections::HashMap::new();
        other_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        mat_other.read_nuclides_from_json(&other_map).unwrap();
        mat_other.set_density("g/cm3", 1.0).unwrap();

        // Both should work and give different results (different nuclides)
        let (xs_file, _) = mat_file.macroscopic_cross_section("(n,gamma)");
        let (xs_other, _) = mat_other.macroscopic_cross_section("(n,gamma)");

        assert!(!xs_file.is_empty() && !xs_other.is_empty(), "Both materials should have cross section data");
        
        // Should be different since Li6 vs Li7
        let data_different = xs_file.len() != xs_other.len() || 
                            xs_file.iter().zip(&xs_other).any(|(a, b)| (a - b).abs() > 1e-10);
        assert!(data_different, "Li6 file vs Li7 file should give different results");
    }

    #[test]
    fn test_material_cache_respects_data_source_boundaries() {
        // Test that material cache properly separates different data sources
        
        crate::nuclide::clear_nuclide_cache();

        // Material 1: Li6 from file first time
        let mut mat_li6_1 = Material::new();
        mat_li6_1.add_nuclide("Li6", 1.0).unwrap();
        let mut li6_map = std::collections::HashMap::new();
        li6_map.insert("Li6".to_string(), "tests/Li6.json".to_string());
        mat_li6_1.read_nuclides_from_json(&li6_map).unwrap();
        mat_li6_1.set_density("g/cm3", 0.534).unwrap();

        // Material 2: Li7 from file
        let mut mat_li7 = Material::new();
        mat_li7.add_nuclide("Li7", 1.0).unwrap();
        let mut li7_map = std::collections::HashMap::new();
        li7_map.insert("Li7".to_string(), "tests/Li7.json".to_string());
        mat_li7.read_nuclides_from_json(&li7_map).unwrap();
        mat_li7.set_density("g/cm3", 0.534).unwrap();

        // Material 3: Li6 from file again (should use cache)
        let mut mat_li6_2 = Material::new();
        mat_li6_2.add_nuclide("Li6", 1.0).unwrap();
        mat_li6_2.read_nuclides_from_json(&li6_map).unwrap(); // reuse same map
        mat_li6_2.set_density("g/cm3", 0.534).unwrap();

        // Get cross sections
        let (xs_li6_1, _) = mat_li6_1.macroscopic_cross_section("(n,gamma)");
        let (xs_li7, _) = mat_li7.macroscopic_cross_section("(n,gamma)");
        let (xs_li6_2, _) = mat_li6_2.macroscopic_cross_section("(n,gamma)");

        // Li6 materials should be identical (cache working)
        assert_eq!(xs_li6_1, xs_li6_2, "Li6 materials should be identical (cache working)");

        // Li6 vs Li7 should be different (different nuclides)
        let li6_vs_li7_different = xs_li6_1.len() != xs_li7.len() || 
                                  xs_li6_1.iter().zip(&xs_li7).any(|(a, b)| (a - b).abs() > 1e-10);
        assert!(li6_vs_li7_different, "Li6 and Li7 materials should have different data");
    }

    #[test]
    fn test_material_path_normalization_in_cache() {
        // Test that different path formats for same file use same cache entry
        
        crate::nuclide::clear_nuclide_cache();

        // Material 1: relative path
        let mut mat_rel = Material::new();
        mat_rel.add_nuclide("Li6", 1.0).unwrap();
        let mut map_rel = std::collections::HashMap::new();
        map_rel.insert("Li6".to_string(), "tests/Li6.json".to_string());
        mat_rel.read_nuclides_from_json(&map_rel).unwrap();
        mat_rel.set_density("g/cm3", 1.0).unwrap();

        // Material 2: absolute path to same file
        let mut mat_abs = Material::new();
        mat_abs.add_nuclide("Li6", 1.0).unwrap();
        let mut map_abs = std::collections::HashMap::new();
        let abs_path = std::env::current_dir().unwrap().join("tests/Li6.json");
        map_abs.insert("Li6".to_string(), abs_path.to_string_lossy().to_string());
        mat_abs.read_nuclides_from_json(&map_abs).unwrap();
        mat_abs.set_density("g/cm3", 1.0).unwrap();

        // Should give identical results (same file, cache working)
        let (xs_rel, _) = mat_rel.macroscopic_cross_section("(n,gamma)");
        let (xs_abs, _) = mat_abs.macroscopic_cross_section("(n,gamma)");

        assert_eq!(xs_rel, xs_abs, "Relative and absolute paths to same file should give identical results");
    }
} // close mod tests