KiThe 0.3.6

//! # User Substances Module
//!
//! This module provides comprehensive management of thermodynamic and transport properties
//! for chemical substances across multiple databases. It implements a priority-based search
//! system and supports multiple data formats including numerical values, functions, and
//! symbolic expressions.
//!
//! ## Key Features
//!
//! - **Multi-Database Search**: Searches across NASA, NIST, CEA, Aramco and other databases
//! - **Priority System**: Hierarchical search with priority and permitted libraries
//! - **Multiple Data Formats**: Values, functions, and symbolic expressions
//! - **Comprehensive Properties**: Heat capacity, enthalpy, entropy, viscosity, thermal conductivity
//! - **Error Management**: Detailed logging and graceful error handling
//! - **Batch Processing**: Efficient handling of multiple substances
//!
//! ## Usage Example
//!
//! ```rust, ignore
//! use crate::Thermodynamics::User_substances::*;
//!
//! let mut subs_data = SubsData::new();
//! subs_data.set_substances(vec!["CO".to_string(), "CO2".to_string()]);
//! subs_data.set_multiple_library_priorities(
//!     vec!["NASA_gas".to_string()],
//!     LibraryPriority::Priority
//! );
//! subs_data.search_substances().unwrap();
//! subs_data.calculate_therm_map_of_properties(298.15).unwrap();
//! ```
/// Enumeration of all supported thermodynamic and transport property types
///
/// This enum covers three formats for each property:
/// - Direct values (e.g., `Cp`): Numerical results at specific conditions
/// - Functions (e.g., `Cp_fun`): Closures for temperature-dependent calculations  
/// - Symbolic (e.g., `Cp_sym`): Mathematical expressions for analytical work
use crate::Thermodynamics::DBhandlers::Diffusion::MultiSubstanceDiffusion;
/// agregator for all the thermo and transport calculations for user defined substances
use crate::Thermodynamics::DBhandlers::thermo_api::{
    ThermoCalculator, ThermoEnum, ThermoError, create_thermal_by_name,
};
use crate::Thermodynamics::DBhandlers::transport_api::{
    TransportCalculator, TransportEnum, create_transport_calculator_by_name,
};
use crate::Thermodynamics::User_substances_error::SimpleExceptionLogger;
use crate::Thermodynamics::User_substances_error::{SubsDataError, SubsDataResult};
pub use crate::Thermodynamics::physical_state::PhysicalState as Phases;
use crate::Thermodynamics::thermo_lib_api::{
    LibraryCapability, LibraryId, ResolvedThermoRecord, ThermoData, ThermoRepository,
};
use std::fmt;

use RustedSciThe::symbolic::symbolic_engine::Expr;
use nalgebra::DMatrix;
use serde_json::Value;
use std::collections::{HashMap, HashSet};
use std::sync::Arc;
/*
 * User Substances Module - Comprehensive Thermodynamic and Transport Property Management
 *
 * This module provides a unified interface for managing thermodynamic and transport properties
 * of chemical substances across multiple databases and libraries. It implements a sophisticated
 * search and priority system for substance data retrieval.
 *
 * Core Architecture:
 * ================
 *
 * Data Types:
 * - DataType: Enumeration of all supported property types (Cp, dH, dS, Lambda, Visc, etc.)
 *   with variants for values, functions, and symbolic expressions
 * - CalculatorType: Wrapper for thermodynamic or transport calculators
 * - WhatIsFound: Classification of search results (Thermo, Transport, NotFound)
 * - SearchResult: Complete search result with library info, priority, and calculator
 * - LibraryPriority: Priority classification (Priority, Permitted, Explicit)
 * - Phases: Physical phase enumeration (Liquid, Gas, Solid, Condensed)
 *
 * Main Structure - SubsData:
 * - Manages substance lists and library priorities
 * - Performs hierarchical searches across multiple databases
 * - Stores calculated properties in multiple formats (values, functions, symbolic)
 * - Provides comprehensive error logging and handling
 * - Supports both thermodynamic and transport property calculations
 *
 * Search Strategy:
 * ===============
 * 1. Priority Libraries: Search high-priority databases first
 * 2. Permitted Libraries: Fallback to secondary databases
 * 3. Explicit Instructions: Direct substance-to-library mappings
 * 4. Result Storage: Maintains search provenance and calculator instances
 *
 * Property Management:
 * ===================
 * - Values: Numerical results at specific conditions
 * - Functions: Closures for temperature-dependent calculations
 * - Symbolic: Mathematical expressions for analytical work
 * - Multi-substance: Batch processing capabilities
 *
 * Supported Libraries:
 * ===================
 * Thermodynamic: NASA, NIST, NASA7, NUIG_thermo, Cantera databases
 * Transport: CEA, Aramco_transport for viscosity, thermal conductivity, diffusion
 *
 * Error Handling:
 * ===============
 * - Comprehensive logging system with substance and function tracking
 * - Graceful degradation for missing data
 * - Detailed error propagation and reporting
 */

#[derive(Debug, Eq, PartialEq, Hash, Clone, Copy)]
#[allow(non_camel_case_types)]
pub enum DataType {
    /// Heat capacity at constant pressure [J/mol·K]
    Cp,
    /// Heat capacity function Cp(T)
    Cp_fun,
    /// Heat capacity symbolic expression
    Cp_sym,
    /// Enthalpy change [J/mol]
    dH,
    /// Enthalpy function dH(T)
    dH_fun,
    /// Enthalpy symbolic expression
    dH_sym,
    /// Entropy change [J/mol·K]
    dS,
    /// Entropy function dS(T)
    dS_fun,
    /// Entropy symbolic expression
    dS_sym,
    /// Chemical potential change [J/mol]
    dmu,
    /// Chemical potential function dmu(T)
    dmu_fun,
    /// Chemical potential symbolic expression
    dmu_sym,
    /// Thermal conductivity [W/m·K]
    Lambda,
    /// Thermal conductivity function Lambda(T)
    Lambda_fun,
    /// Thermal conductivity symbolic expression
    Lambda_sym,
    /// Dynamic viscosity [Pa·s]
    Visc,
    /// Viscosity function Visc(T)
    Visc_fun,
    /// Viscosity symbolic expression
    Visc_sym,
}
#[allow(non_camel_case_types)]
/// Canonical property identity without representation-specific suffixes.
#[derive(Debug, Eq, PartialEq, Hash, Clone, Copy)]
pub enum PropertyKind {
    /// Heat capacity at constant pressure.
    Cp,
    /// Enthalpy.
    dH,
    /// Entropy.
    dS,
    /// Chemical potential.
    dmu,
    /// Thermal conductivity.
    Lambda,
    /// Dynamic viscosity.
    Visc,
}

impl PropertyKind {
    /// Returns the numeric representation key for this property kind.
    pub fn numeric_type(self) -> DataType {
        match self {
            PropertyKind::Cp => DataType::Cp,
            PropertyKind::dH => DataType::dH,
            PropertyKind::dS => DataType::dS,
            PropertyKind::dmu => DataType::dmu,
            PropertyKind::Lambda => DataType::Lambda,
            PropertyKind::Visc => DataType::Visc,
        }
    }

    /// Returns the closure representation key for this property kind.
    pub fn function_type(self) -> DataType {
        match self {
            PropertyKind::Cp => DataType::Cp_fun,
            PropertyKind::dH => DataType::dH_fun,
            PropertyKind::dS => DataType::dS_fun,
            PropertyKind::dmu => DataType::dmu_fun,
            PropertyKind::Lambda => DataType::Lambda_fun,
            PropertyKind::Visc => DataType::Visc_fun,
        }
    }

    /// Returns the symbolic representation key for this property kind.
    pub fn symbolic_type(self) -> DataType {
        match self {
            PropertyKind::Cp => DataType::Cp_sym,
            PropertyKind::dH => DataType::dH_sym,
            PropertyKind::dS => DataType::dS_sym,
            PropertyKind::dmu => DataType::dmu_sym,
            PropertyKind::Lambda => DataType::Lambda_sym,
            PropertyKind::Visc => DataType::Visc_sym,
        }
    }

    /// Returns `true` when the property belongs to the transport branch.
    pub fn is_transport(self) -> bool {
        matches!(self, PropertyKind::Lambda | PropertyKind::Visc)
    }

    /// Returns the canonical identity from a representation-specific key.
    pub fn from_data_type(data_type: DataType) -> Option<Self> {
        match data_type {
            DataType::Cp | DataType::Cp_fun | DataType::Cp_sym => Some(PropertyKind::Cp),
            DataType::dH | DataType::dH_fun | DataType::dH_sym => Some(PropertyKind::dH),
            DataType::dS | DataType::dS_fun | DataType::dS_sym => Some(PropertyKind::dS),
            DataType::dmu | DataType::dmu_fun | DataType::dmu_sym => Some(PropertyKind::dmu),
            DataType::Lambda | DataType::Lambda_fun | DataType::Lambda_sym => {
                Some(PropertyKind::Lambda)
            }
            DataType::Visc | DataType::Visc_fun | DataType::Visc_sym => Some(PropertyKind::Visc),
        }
    }
}

/// Wrapper enum for different types of property calculators
///
/// Encapsulates either thermodynamic calculators (for Cp, dH, dS) or
/// transport calculators (for viscosity, thermal conductivity, diffusion)
#[derive(Debug, Clone)]
pub enum CalculatorType {
    /// Thermodynamic property calculator (NASA, NIST, etc.)
    Thermo(ThermoEnum),
    /// Transport property calculator (CEA, Aramco, etc.)
    Transport(TransportEnum),
}
/// Classification of search results for substance data
///
/// Indicates what type of data was found during library searches
#[derive(Debug, Copy, Clone, Eq, Hash, PartialEq)]
pub enum WhatIsFound {
    /// Thermodynamic data found (Cp, dH, dS)
    Thermo,
    /// Transport data found (viscosity, thermal conductivity)
    Transport,
    /// No data found in any searched library
    NotFound,
}

/// Canonical per-property search state for one substance.
#[derive(Debug, Clone)]
pub enum PropertySearchState {
    NotSearched,
    Found(SearchResult),
    Missing,
    Failed(String),
}

impl PropertySearchState {
    /// Returns `true` when the property has not been searched yet.
    pub fn is_pending(&self) -> bool {
        matches!(self, PropertySearchState::NotSearched)
    }

    /// Returns `true` when the property has been resolved successfully.
    pub fn is_resolved(&self) -> bool {
        matches!(self, PropertySearchState::Found(_))
    }

    /// Returns `true` when the property lookup ended in a missing-record state.
    pub fn is_missing(&self) -> bool {
        matches!(self, PropertySearchState::Missing)
    }

    /// Returns `true` when the property lookup failed with an error message.
    pub fn is_failed(&self) -> bool {
        matches!(self, PropertySearchState::Failed(_))
    }
}

/// Explicit readiness state for a derived numeric cache entry.
///
/// This keeps the public read-only API from confusing a missing cache entry
/// with a valid numeric zero.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum DerivedValueState<'a, T: ?Sized> {
    /// No calculated value is available yet.
    NotCalculated,
    /// A concrete calculated value exists.
    Ready(&'a T),
}

impl<'a, T: ?Sized> DerivedValueState<'a, T> {
    /// Returns `true` when the entry has not been calculated yet.
    pub fn is_not_calculated(&self) -> bool {
        matches!(self, DerivedValueState::NotCalculated)
    }

    /// Returns `true` when the entry contains a calculated value.
    pub fn is_ready(&self) -> bool {
        matches!(self, DerivedValueState::Ready(_))
    }

    /// Returns the calculated value when it exists.
    pub fn value(&self) -> Option<&'a T> {
        match self {
            DerivedValueState::Ready(value) => Some(*value),
            DerivedValueState::NotCalculated => None,
        }
    }
}

/// Transport calculation input mode.
///
/// The CEA branch does not require pressure or molar mass inputs, so we keep
/// that branch explicit instead of smuggling dummy values through the solver.
#[derive(Debug, Clone)]
pub(crate) enum TransportCalculationMode {
    Cea,
    Standard {
        pressure: f64,
        molar_mass: f64,
        pressure_unit: Option<String>,
        molar_mass_unit: Option<String>,
    },
}

/// Canonical search state for a single substance.
#[derive(Debug, Clone)]
pub struct SubstanceSearchState {
    pub thermo: PropertySearchState,
    pub transport: PropertySearchState,
}

impl SubstanceSearchState {
    /// Returns the canonical thermo search state.
    pub fn thermo(&self) -> &PropertySearchState {
        &self.thermo
    }

    /// Returns the canonical transport search state.
    pub fn transport(&self) -> &PropertySearchState {
        &self.transport
    }

    /// Returns the canonical state for one property, if the property is known.
    pub fn property_state(&self, found_kind: WhatIsFound) -> Option<&PropertySearchState> {
        match found_kind {
            WhatIsFound::Thermo => Some(&self.thermo),
            WhatIsFound::Transport => Some(&self.transport),
            WhatIsFound::NotFound => None,
        }
    }

    /// Returns the canonical found record for one property, if present.
    pub fn result(&self, found_kind: WhatIsFound) -> Option<&SearchResult> {
        match self.property_state(found_kind)? {
            PropertySearchState::Found(result) => Some(result),
            PropertySearchState::NotSearched
            | PropertySearchState::Missing
            | PropertySearchState::Failed(_) => None,
        }
    }

    /// Returns a failure message for one property, if present.
    pub fn failure(&self, found_kind: WhatIsFound) -> Option<&str> {
        match self.property_state(found_kind)? {
            PropertySearchState::Failed(message) => Some(message.as_str()),
            PropertySearchState::NotSearched
            | PropertySearchState::Found(_)
            | PropertySearchState::Missing => None,
        }
    }

    pub(crate) fn new() -> Self {
        Self {
            thermo: PropertySearchState::NotSearched,
            transport: PropertySearchState::NotSearched,
        }
    }

    pub(crate) fn set_found(
        &mut self,
        found_kind: WhatIsFound,
        result: SearchResult,
        overwrite: bool,
    ) {
        let target = match found_kind {
            WhatIsFound::Thermo => &mut self.thermo,
            WhatIsFound::Transport => &mut self.transport,
            WhatIsFound::NotFound => return,
        };

        if overwrite
            || matches!(
                target,
                PropertySearchState::NotSearched | PropertySearchState::Missing
            )
        {
            *target = PropertySearchState::Found(result);
        }
    }

    pub(crate) fn mark_missing(&mut self) {
        if matches!(self.thermo, PropertySearchState::NotSearched)
            && matches!(self.transport, PropertySearchState::NotSearched)
        {
            self.thermo = PropertySearchState::Missing;
            self.transport = PropertySearchState::Missing;
        }
    }

    pub(crate) fn to_compat_map(&self) -> HashMap<WhatIsFound, Option<SearchResult>> {
        let mut map = HashMap::new();
        match self.thermo() {
            PropertySearchState::Found(result) => {
                map.insert(WhatIsFound::Thermo, Some(result.clone()));
            }
            PropertySearchState::Missing | PropertySearchState::Failed(_) => {
                map.insert(WhatIsFound::NotFound, None);
            }
            PropertySearchState::NotSearched => {}
        }
        match self.transport() {
            PropertySearchState::Found(result) => {
                map.insert(WhatIsFound::Transport, Some(result.clone()));
            }
            PropertySearchState::Missing | PropertySearchState::Failed(_) => {
                map.entry(WhatIsFound::NotFound).or_insert(None);
            }
            PropertySearchState::NotSearched => {}
        }
        map
    }

    /// Returns `true` when both properties were explicitly marked as missing.
    pub(crate) fn is_missing(&self) -> bool {
        matches!(self.thermo, PropertySearchState::Missing)
            && matches!(self.transport, PropertySearchState::Missing)
    }

    /// Returns `true` when at least one property was found in a priority library.
    pub(crate) fn has_priority_found(&self) -> bool {
        matches!(
            &self.thermo,
            PropertySearchState::Found(SearchResult {
                priority_type: LibraryPriority::Priority,
                ..
            })
        ) || matches!(
            &self.transport,
            PropertySearchState::Found(SearchResult {
                priority_type: LibraryPriority::Priority,
                ..
            })
        )
    }

    /// Returns `true` when at least one property was found in a permitted library.
    pub(crate) fn has_permitted_found(&self) -> bool {
        matches!(
            &self.thermo,
            PropertySearchState::Found(SearchResult {
                priority_type: LibraryPriority::Permitted,
                ..
            })
        ) || matches!(
            &self.transport,
            PropertySearchState::Found(SearchResult {
                priority_type: LibraryPriority::Permitted,
                ..
            })
        )
    }
}
/// Complete search result for a substance in a specific library
///
/// Contains all information needed to work with found substance data,
/// including the source library, priority level, raw data, and initialized calculator
#[derive(Debug, Clone)]
pub struct SearchResult {
    /// Name of the library where the substance was found.
    pub(crate) library: String,
    /// Exact key selected from the source library.
    pub(crate) record_key: String,
    /// Priority level of the library (Priority, Permitted, Explicit).
    pub(crate) priority_type: LibraryPriority,
    /// Raw JSON data from the library.
    pub(crate) data: Value,
    /// Initialized calculator instance for property calculations.
    pub(crate) calculator: Option<CalculatorType>,
}

impl SearchResult {
    /// Returns the source library name for the canonical record.
    pub fn library(&self) -> &str {
        &self.library
    }

    /// Returns the exact source-record key, which may include a state tag.
    pub fn record_key(&self) -> &str {
        &self.record_key
    }

    /// Returns the library priority captured when the record was found.
    pub fn priority_type(&self) -> LibraryPriority {
        self.priority_type
    }

    /// Returns the raw JSON payload captured from the library.
    pub fn data(&self) -> &Value {
        &self.data
    }

    /// Returns the attached calculator, if present.
    pub fn calculator(&self) -> Option<&CalculatorType> {
        self.calculator.as_ref()
    }

    /// Returns the attached calculator mutably.
    pub(crate) fn calculator_mut(&mut self) -> Option<&mut CalculatorType> {
        self.calculator.as_mut()
    }
}
/// Library priority levels for hierarchical searching
///
/// Defines the search order and preference for different databases
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum LibraryPriority {
    /// High-priority libraries searched first (most trusted/accurate data)
    Priority,
    /// Secondary libraries used as fallback options
    Permitted,
    /// Direct substance-to-library mapping (bypasses normal search)
    Explicit,
}

/// Main structure for comprehensive substance property management
///
/// `SubsData` is the central hub for managing thermodynamic and transport properties
/// of chemical substances. It provides:
///
/// - **Multi-database integration**: Searches across NASA, NIST, CEA, Aramco databases
/// - **Priority-based search**: Hierarchical search with configurable library priorities
/// - **Multiple data formats**: Stores properties as values, functions, and symbolic expressions
/// - **Comprehensive error handling**: Detailed logging and graceful error recovery
/// - **Batch processing**: Efficient handling of multiple substances simultaneously
///
/// ## Key Components
///
/// ### Search System
/// - `library_priorities`: Maps libraries to priority levels
/// - `search_results`: Stores found data with provenance information
/// - `explicit_search_instructions`: Direct substance-to-library mappings
///
/// ### Property Storage
/// - `therm_map_of_properties_values`: Numerical thermodynamic properties
/// - `therm_map_of_fun`: Function closures for temperature-dependent calculations
/// - `therm_map_of_sym`: Symbolic expressions for analytical work
/// - `transport_map_*`: Equivalent storage for transport properties
///
/// ### Physical Properties
/// - `P`, `T`: Pressure and temperature conditions
/// - `map_of_phases`: Physical phase information for each substance
/// - `molar_mass_by_substance`: Molecular weights
/// - `elem_composition_matrix`: Elemental composition data
///
/// ## Usage Pattern
///
/// 1. Create instance and set substances
/// 2. Configure library priorities
/// 3. Perform search across databases
/// 4. Calculate properties in desired formats
/// 5. Access results through getter methods
pub struct SubsData {
    /// System pressure [Pa] for property calculations
    pub(crate) P: Option<f64>,
    /// Pressure unit (e.g., "Pa", "atm", "bar")
    pub(crate) P_unit: Option<String>,
    /// System temperature [K] for property calculations
    pub(crate) T: Option<f64>,
    /// Molar mass unit (e.g., "g/mol", "kg/mol")
    pub(crate) Molar_mass_unit: Option<String>,
    /// List of chemical substances to search for and analyze
    pub(crate) substances: Vec<String>,
    /// Maps library names to their search priority levels
    pub(crate) library_priorities: HashMap<String, LibraryPriority>,
    /// Direct substance-to-library mappings that bypass normal search hierarchy
    pub(crate) explicit_search_instructions: HashMap<String, String>,
    /// Density values for each substance [kg/m³]
    pub(crate) ro_map: Option<HashMap<String, f64>>,
    /// Symbolic expressions for density calculations
    pub(crate) ro_map_sym: Option<HashMap<String, Box<Expr>>>,
    /// Physical phase context for each substance, used by phase models and
    /// live NIST fallback. It is not a local-library lookup filter.
    pub(crate) map_of_phases: HashMap<String, Option<Phases>>,
    /// Optional explicit physical-state filters for local thermo-record lookup.
    /// An absent entry means "search all records by the requested substance
    /// name", preserving the normal unconstrained lookup contract.
    physical_state_requirements: HashMap<String, Phases>,
    /// Canonical search state for each substance.
    pub(crate) search_states: HashMap<String, SubstanceSearchState>,
    /// Complete search results with library provenance and calculators
    pub(crate) search_results: HashMap<String, HashMap<WhatIsFound, Option<SearchResult>>>,
    /// Interface to all thermodynamic and transport databases
    pub(crate) thermo_data: ThermoData,
    /// Calculated thermodynamic property values at specific conditions
    pub(crate) therm_map_of_properties_values: HashMap<String, HashMap<DataType, Option<f64>>>,
    /// Function closures for temperature-dependent thermodynamic calculations
    pub(crate) therm_map_of_fun:
        HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>>,
    /// Symbolic expressions for analytical thermodynamic calculations
    pub(crate) therm_map_of_sym: HashMap<String, HashMap<DataType, Option<Box<Expr>>>>,
    /// Calculated transport property values at specific conditions
    pub(crate) transport_map_of_properties_values: HashMap<String, HashMap<DataType, Option<f64>>>,
    /// Function closures for temperature-dependent transport calculations
    pub(crate) transport_map_of_fun:
        HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>>,
    /// Symbolic expressions for analytical transport calculations
    pub(crate) transport_map_of_sym: HashMap<String, HashMap<DataType, Option<Box<Expr>>>>,
    /// Matrix of elemental composition for all substances
    pub(crate) elem_composition_matrix: Option<DMatrix<f64>>,
    /// Molar masses for each substance [g/mol or specified unit]
    pub(crate) molar_mass_by_substance: HashMap<String, f64>,
    /// Multi-component diffusion coefficient data and calculations
    pub(crate) diffusion_data: Option<MultiSubstanceDiffusion>,
    /// List of unique chemical elements present in all substances
    pub(crate) unique_elements: Vec<String>,
    /// Error logging system for tracking calculation failures
    pub(crate) logger: SimpleExceptionLogger,
}
impl fmt::Debug for SubsData {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("SubsData")
            .field("substances", &self.substances)
            .field("library_priorities", &self.library_priorities)
            .field("search_states", &self.search_states)
            .field("search_results", &self.search_results)
            // Skip thermo_data as it might be large
            .field(
                "therm_map_of_properties_values",
                &self.therm_map_of_properties_values,
            )
            // Skip therm_map_of_fun as it contains closures that don't implement Debug
            .field("therm_map_of_sym", &self.therm_map_of_sym)
            .finish()
    }
}

/// Clone the shape of a cached function map without attempting to duplicate closures.
///
/// The cache remains structurally intact so callers can rebuild it explicitly after cloning.
fn clone_function_cache(
    source: &HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>>,
) -> HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>> {
    let mut cloned = HashMap::with_capacity(source.len());

    for (substance, type_map) in source {
        let mut cloned_type_map = HashMap::with_capacity(type_map.len());
        for (data_type, _) in type_map {
            cloned_type_map.insert(*data_type, None);
        }
        cloned.insert(substance.clone(), cloned_type_map);
    }

    cloned
}

impl Clone for SubsData {
    fn clone(&self) -> Self {
        // Closure caches are intentionally reset during clone.
        // Callers can rebuild them explicitly when they actually need executable functions.
        let new_therm_map_of_fun = clone_function_cache(&self.therm_map_of_fun);
        let new_transport_map_of_fun = clone_function_cache(&self.transport_map_of_fun);

        let sd = SubsData {
            P: self.P,
            P_unit: self.P_unit.clone(),
            T: self.T.clone(),
            Molar_mass_unit: self.Molar_mass_unit.clone(),
            substances: self.substances.clone(),
            library_priorities: self.library_priorities.clone(),
            explicit_search_instructions: self.explicit_search_instructions.clone(),
            ro_map: self.ro_map.clone(),
            ro_map_sym: self.ro_map_sym.clone(),
            map_of_phases: self.map_of_phases.clone(),
            physical_state_requirements: self.physical_state_requirements.clone(),
            search_states: self.search_states.clone(),
            search_results: self.get_all_results(),
            thermo_data: self.thermo_data.clone(),
            therm_map_of_properties_values: self.therm_map_of_properties_values.clone(),
            therm_map_of_fun: new_therm_map_of_fun,
            therm_map_of_sym: self.therm_map_of_sym.clone(),
            transport_map_of_properties_values: self.transport_map_of_properties_values.clone(),
            transport_map_of_fun: new_transport_map_of_fun,
            transport_map_of_sym: self.transport_map_of_sym.clone(),
            elem_composition_matrix: self.elem_composition_matrix.clone(),
            molar_mass_by_substance: self.molar_mass_by_substance.clone(),
            diffusion_data: self.diffusion_data.clone(),
            unique_elements: self.unique_elements.clone(),
            logger: self.logger.clone(),
        };
        sd
    }
}
impl SubsData {
    /// Creates an empty query bound to an explicitly supplied immutable catalog.
    ///
    /// This is the construction path for higher-level systems that resolve
    /// multiple phase payloads under one repository lifecycle.
    pub fn from_thermo_repository(repository: Arc<ThermoRepository>) -> Self {
        let mut data = Self::empty();
        data.thermo_data = ThermoData::from_repository(repository);
        data
    }

    /// Returns the current query substances as a read-only slice.
    pub fn substances(&self) -> &[String] {
        &self.substances
    }

    /// Returns the configured library priority map.
    pub fn library_priorities(&self) -> &HashMap<String, LibraryPriority> {
        &self.library_priorities
    }

    /// Returns the direct substance-to-library instructions.
    pub fn explicit_search_instructions(&self) -> &HashMap<String, String> {
        &self.explicit_search_instructions
    }

    /// Returns the direct substance-to-library instructions under a shorter
    /// canonical name.
    pub fn explicit_search_map(&self) -> &HashMap<String, String> {
        &self.explicit_search_instructions
    }

    /// Returns the configured pressure value, if present.
    pub fn pressure(&self) -> Option<f64> {
        self.P
    }

    /// Returns the configured pressure unit, if present.
    pub fn pressure_unit(&self) -> Option<&str> {
        self.P_unit.as_deref()
    }

    /// Returns the configured temperature value, if present.
    pub fn temperature(&self) -> Option<f64> {
        self.T
    }

    /// Returns the configured molar-mass unit, if present.
    pub fn molar_mass_unit(&self) -> Option<&str> {
        self.Molar_mass_unit.as_deref()
    }

    /// Returns the configured phase hints for substances.
    pub fn phases(&self) -> &HashMap<String, Option<Phases>> {
        &self.map_of_phases
    }

    /// Returns the canonical lookup state map.
    pub fn search_states(&self) -> &HashMap<String, SubstanceSearchState> {
        &self.search_states
    }

    /// Returns the canonical thermodynamic value cache.
    pub fn therm_values(&self) -> &HashMap<String, HashMap<DataType, Option<f64>>> {
        &self.therm_map_of_properties_values
    }

    /// Returns the canonical thermodynamic closure cache.
    pub fn therm_functions(
        &self,
    ) -> &HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>> {
        &self.therm_map_of_fun
    }

    /// Returns the readiness state of a thermodynamic closure cache entry.
    pub fn therm_function_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, dyn Fn(f64) -> f64 + Send + Sync> {
        self.derived_function_state(&self.therm_map_of_fun, substance, data_type)
    }

    /// Returns the canonical thermodynamic symbolic cache.
    pub fn therm_symbolic(&self) -> &HashMap<String, HashMap<DataType, Option<Box<Expr>>>> {
        &self.therm_map_of_sym
    }

    /// Returns the readiness state for a canonical property kind.
    pub fn property_value_state(
        &self,
        substance: &str,
        kind: PropertyKind,
    ) -> DerivedValueState<'_, f64> {
        if kind.is_transport() {
            self.transport_value_state(substance, kind.numeric_type())
        } else {
            self.therm_value_state(substance, kind.numeric_type())
        }
    }

    /// Returns the readiness state for a canonical property closure.
    pub fn property_function_state(
        &self,
        substance: &str,
        kind: PropertyKind,
    ) -> DerivedValueState<'_, dyn Fn(f64) -> f64 + Send + Sync> {
        if kind.is_transport() {
            self.transport_function_state(substance, kind.function_type())
        } else {
            self.therm_function_state(substance, kind.function_type())
        }
    }

    /// Returns the readiness state for a canonical property symbolic expression.
    pub fn property_symbolic_state(
        &self,
        substance: &str,
        kind: PropertyKind,
    ) -> DerivedValueState<'_, Expr> {
        if kind.is_transport() {
            self.transport_symbolic_state(substance, kind.symbolic_type())
        } else {
            self.therm_symbolic_state(substance, kind.symbolic_type())
        }
    }

    /// Returns the canonical transport value cache.
    pub fn transport_values(&self) -> &HashMap<String, HashMap<DataType, Option<f64>>> {
        &self.transport_map_of_properties_values
    }

    /// Returns the readiness state of a thermodynamic numeric cache entry.
    pub fn therm_value_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, f64> {
        self.derived_value_state(&self.therm_map_of_properties_values, substance, data_type)
    }

    /// Returns the readiness state of a transport numeric cache entry.
    pub fn transport_value_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, f64> {
        self.derived_value_state(
            &self.transport_map_of_properties_values,
            substance,
            data_type,
        )
    }

    /// Returns the readiness state of a thermodynamic symbolic cache entry.
    pub fn therm_symbolic_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, Expr> {
        self.derived_boxed_value_state(&self.therm_map_of_sym, substance, data_type)
    }

    /// Returns the readiness state of a transport symbolic cache entry.
    pub fn transport_symbolic_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, Expr> {
        self.derived_boxed_value_state(&self.transport_map_of_sym, substance, data_type)
    }

    /// Returns the canonical transport closure cache.
    pub fn transport_functions(
        &self,
    ) -> &HashMap<String, HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>> {
        &self.transport_map_of_fun
    }

    /// Returns the readiness state of a transport closure cache entry.
    pub fn transport_function_state(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'_, dyn Fn(f64) -> f64 + Send + Sync> {
        self.derived_function_state(&self.transport_map_of_fun, substance, data_type)
    }

    /// Returns the canonical transport symbolic cache.
    pub fn transport_symbolic(&self) -> &HashMap<String, HashMap<DataType, Option<Box<Expr>>>> {
        &self.transport_map_of_sym
    }

    /// Returns the molar mass map.
    pub fn molar_masses(&self) -> &HashMap<String, f64> {
        &self.molar_mass_by_substance
    }

    /// Returns the molar mass map under a canonical accessor name.
    pub fn molar_mass_map(&self) -> &HashMap<String, f64> {
        &self.molar_mass_by_substance
    }

    /// Returns a single molar mass by substance name.
    pub fn molar_mass_of(&self, substance: &str) -> Option<f64> {
        self.molar_mass_by_substance.get(substance).copied()
    }

    /// Returns the unique elements list.
    pub fn unique_elements(&self) -> &[String] {
        &self.unique_elements
    }

    /// Returns the element composition matrix, if it has been built.
    pub fn element_composition_matrix(&self) -> Option<&DMatrix<f64>> {
        self.elem_composition_matrix.as_ref()
    }

    /// Returns the diffusion model, if it has been initialized.
    pub fn diffusion_data(&self) -> Option<&MultiSubstanceDiffusion> {
        self.diffusion_data.as_ref()
    }

    /// Returns the immutable thermodynamic and transport repository handle.
    pub fn thermo_data(&self) -> &ThermoData {
        &self.thermo_data
    }

    /// Returns the diagnostics sink.
    pub fn logger(&self) -> &SimpleExceptionLogger {
        &self.logger
    }

    /// Maps a nested cache entry to an explicit readiness state.
    fn derived_value_state<'a>(
        &'a self,
        cache: &'a HashMap<String, HashMap<DataType, Option<f64>>>,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'a, f64> {
        cache
            .get(substance)
            .and_then(|property_map| property_map.get(&data_type))
            .and_then(|value| value.as_ref())
            .map_or(DerivedValueState::NotCalculated, DerivedValueState::Ready)
    }

    /// Maps a nested boxed-expression cache entry to an explicit readiness state.
    fn derived_boxed_value_state<'a, T>(
        &'a self,
        cache: &'a HashMap<String, HashMap<DataType, Option<Box<T>>>>,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'a, T> {
        cache
            .get(substance)
            .and_then(|property_map| property_map.get(&data_type))
            .and_then(|value| value.as_ref())
            .map(|value| value.as_ref())
            .map_or(DerivedValueState::NotCalculated, DerivedValueState::Ready)
    }

    /// Maps a nested boxed-closure cache entry to an explicit readiness state.
    fn derived_function_state<'a>(
        &'a self,
        cache: &'a HashMap<
            String,
            HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>,
        >,
        substance: &str,
        data_type: DataType,
    ) -> DerivedValueState<'a, dyn Fn(f64) -> f64 + Send + Sync> {
        cache
            .get(substance)
            .and_then(|property_map| property_map.get(&data_type))
            .and_then(|value| value.as_deref())
            .map_or(DerivedValueState::NotCalculated, DerivedValueState::Ready)
    }

    pub(crate) fn sync_search_state_view(&mut self, substance: &str) {
        if let Some(state) = self.search_states.get(substance) {
            self.search_results
                .insert(substance.to_string(), state.to_compat_map());
        }
    }

    pub(crate) fn store_search_result(
        &mut self,
        substance: &str,
        found_kind: WhatIsFound,
        library: String,
        priority_type: LibraryPriority,
        data: Value,
        calculator: CalculatorType,
        overwrite: bool,
    ) {
        let search_result = SearchResult {
            library,
            record_key: substance.to_string(),
            priority_type,
            data,
            calculator: Some(calculator),
        };
        let state = self
            .search_states
            .entry(substance.to_string())
            .or_insert_with(SubstanceSearchState::new);
        state.set_found(found_kind, search_result, overwrite);
        self.sync_search_state_view(substance);
    }

    /// Stores a result selected through physical-state-aware record lookup.
    fn store_resolved_search_result(
        &mut self,
        substance: &str,
        found_kind: WhatIsFound,
        library: String,
        record_key: String,
        priority_type: LibraryPriority,
        data: Value,
        calculator: CalculatorType,
        overwrite: bool,
    ) {
        let search_result = SearchResult {
            library,
            record_key,
            priority_type,
            data,
            calculator: Some(calculator),
        };
        let state = self
            .search_states
            .entry(substance.to_string())
            .or_insert_with(SubstanceSearchState::new);
        state.set_found(found_kind, search_result, overwrite);
        self.sync_search_state_view(substance);
    }

    /// Sets an explicit physical-state requirement for one requested substance.
    ///
    /// Changing the requirement changes record selection, therefore all
    /// lookup-derived values are invalidated transactionally before the next
    /// search. A state is a selection constraint, not a request to mutate the
    /// JSON record name by hand.
    pub fn set_substance_physical_state(&mut self, substance: String, state: Phases) {
        self.physical_state_requirements.insert(substance, state);
        self.invalidate_lookup_dependent_state();
    }

    /// Replaces all physical-state requirements as one lookup transaction.
    ///
    /// This is the phase-system entry point: it prevents a caller from
    /// accidentally publishing search results after only a subset of one
    /// phase's component constraints has been updated.
    pub fn set_physical_state_requirements(&mut self, states: HashMap<String, Phases>) {
        self.physical_state_requirements = states;
        self.invalidate_lookup_dependent_state();
    }

    /// Clears an explicit physical-state requirement for one substance.
    pub fn clear_substance_physical_state(&mut self, substance: &str) {
        if self.physical_state_requirements.remove(substance).is_some() {
            self.invalidate_lookup_dependent_state();
        }
    }

    /// Returns the current physical-state requirement, if one was set.
    pub fn substance_physical_state(&self, substance: &str) -> Option<Phases> {
        self.physical_state_requirements.get(substance).copied()
    }

    fn resolve_library_record(
        &self,
        library: &str,
        substance: &str,
    ) -> SubsDataResult<Option<ResolvedThermoRecord>> {
        let query = crate::Thermodynamics::physical_state::ThermoRecordQuery::new(substance)
            .with_physical_state_opt(self.substance_physical_state(substance));
        let resolved = self
            .thermo_data
            .resolve_thermo_record(library, &query)
            .map_err(|error| SubsDataError::CalculationFailed {
                substance: substance.to_string(),
                operation: "physical-state-aware library lookup".to_string(),
                reason: error.to_string(),
                source: None,
            })?;
        Ok(resolved)
    }

    pub(crate) fn insert_not_found_if_absent(&mut self, substance: &str) {
        let state = self
            .search_states
            .entry(substance.to_string())
            .or_insert_with(SubstanceSearchState::new);
        state.mark_missing();
        self.sync_search_state_view(substance);
    }

    pub fn new() -> Self {
        Self {
            P: None,
            P_unit: None,
            T: None,
            Molar_mass_unit: None,
            substances: Vec::new(),
            library_priorities: HashMap::new(),
            explicit_search_instructions: HashMap::new(),
            search_states: HashMap::new(),
            search_results: HashMap::new(),
            ro_map: None,
            ro_map_sym: None,
            map_of_phases: HashMap::new(),
            physical_state_requirements: HashMap::new(),
            thermo_data: ThermoData::new(),
            therm_map_of_properties_values: HashMap::new(),
            therm_map_of_fun: HashMap::new(),
            therm_map_of_sym: HashMap::new(),
            transport_map_of_properties_values: HashMap::new(),
            transport_map_of_fun: HashMap::new(),
            transport_map_of_sym: HashMap::new(),
            elem_composition_matrix: None,
            molar_mass_by_substance: HashMap::new(),
            diffusion_data: None,
            unique_elements: Vec::new(),
            logger: SimpleExceptionLogger::new(),
        }
    }

    pub fn empty() -> Self {
        Self {
            P: None,
            P_unit: None,
            T: None,
            Molar_mass_unit: None,
            substances: Vec::new(),
            library_priorities: HashMap::new(),
            explicit_search_instructions: HashMap::new(),
            search_states: HashMap::new(),
            search_results: HashMap::new(),
            ro_map: None,
            ro_map_sym: None,
            map_of_phases: HashMap::new(),
            physical_state_requirements: HashMap::new(),
            thermo_data: ThermoData::empty(),
            therm_map_of_properties_values: HashMap::new(),
            therm_map_of_fun: HashMap::new(),
            therm_map_of_sym: HashMap::new(),
            transport_map_of_properties_values: HashMap::new(),
            transport_map_of_fun: HashMap::new(),
            transport_map_of_sym: HashMap::new(),
            elem_composition_matrix: None,
            molar_mass_by_substance: HashMap::new(),
            diffusion_data: None,
            unique_elements: Vec::new(),
            logger: SimpleExceptionLogger::new(),
        }
    }
    ////////////////////////////SETTERS////////////////////////////////////
    pub fn set_substances(&mut self, substances: Vec<String>) {
        self.invalidate_substance_dependent_state();
        self.substances = substances;
    }

    /// Dependency graph for cache invalidation:
    /// - substances -> lookup results, thermo caches, transport caches,
    ///   phase maps, composition, molar mass, diffusion, derived views
    /// - lookup instructions -> search results and every derived cache
    /// - temperature -> thermo and transport property caches
    /// - pressure / supplied molar masses -> transport caches and diffusion
    /// This keeps the invalidation policy explicit and centralized.
    /// Clears all caches that are tied to the current substance list.
    ///
    /// This keeps stale search results and derived properties from surviving
    /// after the user changes the set of substances under analysis.
    fn invalidate_substance_dependent_state(&mut self) {
        self.search_states.clear();
        self.search_results.clear();
        self.therm_map_of_properties_values.clear();
        self.therm_map_of_fun.clear();
        self.therm_map_of_sym.clear();
        self.transport_map_of_properties_values.clear();
        self.transport_map_of_fun.clear();
        self.transport_map_of_sym.clear();
        self.elem_composition_matrix = None;
        self.molar_mass_by_substance.clear();
        self.diffusion_data = None;
        self.unique_elements.clear();
        self.ro_map = None;
        self.ro_map_sym = None;
        self.map_of_phases.clear();
        self.physical_state_requirements.clear();
    }

    /// Clears lookup-related caches after changing library routing rules.
    ///
    /// Library priorities and explicit lookup instructions affect which
    /// records are visible, so search results and every derived cache must be
    /// rebuilt after those inputs change.
    fn invalidate_lookup_dependent_state(&mut self) {
        self.search_states.clear();
        self.search_results.clear();
        self.therm_map_of_properties_values.clear();
        self.therm_map_of_fun.clear();
        self.therm_map_of_sym.clear();
        self.transport_map_of_properties_values.clear();
        self.transport_map_of_fun.clear();
        self.transport_map_of_sym.clear();
        self.elem_composition_matrix = None;
        self.molar_mass_by_substance.clear();
        self.diffusion_data = None;
        self.unique_elements.clear();
        self.ro_map = None;
        self.ro_map_sym = None;
    }

    /// Clears caches that depend on temperature-dependent calculations.
    fn invalidate_temperature_dependent_state(&mut self) {
        self.therm_map_of_properties_values.clear();
        self.transport_map_of_properties_values.clear();
        self.transport_map_of_fun.clear();
        self.transport_map_of_sym.clear();
        self.ro_map = None;
        self.ro_map_sym = None;
        self.diffusion_data = None;
    }

    /// Clears caches that depend on pressure or manually supplied molar masses.
    fn invalidate_pressure_and_molar_mass_dependent_state(&mut self) {
        self.transport_map_of_properties_values.clear();
        self.transport_map_of_fun.clear();
        self.transport_map_of_sym.clear();
        self.ro_map = None;
        self.ro_map_sym = None;
        self.diffusion_data = None;
    }

    ////////////////////////SEARCH AND PRIORITY HANDLING///////////////////////////////

    /// Sets the priority level for a single library
    ///
    /// Libraries with `Priority` level are searched first, followed by `Permitted` libraries.
    /// `Explicit` priority is used for direct substance-to-library mappings.
    ///
    /// # Arguments
    ///
    /// * `library` - Name of the library (e.g., "NASA_gas", "CEA", "Aramco_transport")
    /// * `priority` - Priority level for this library
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_library_priority("NASA_gas".to_string(), LibraryPriority::Priority);
    /// subs_data.set_library_priority("NIST".to_string(), LibraryPriority::Permitted);
    /// ```
    pub fn set_library_priority(&mut self, library: String, priority: LibraryPriority) {
        let canonical_library = ThermoData::canonical_library_name(&library);
        self.library_priorities.insert(canonical_library, priority);
        self.invalidate_lookup_dependent_state();
    }

    /// Sets the priority level for a typed library identifier.
    pub fn set_library_priority_id(&mut self, library: LibraryId, priority: LibraryPriority) {
        self.set_library_priority(library.canonical_name().to_string(), priority);
    }

    /// Sets the same priority level for multiple libraries simultaneously
    ///
    /// Convenient method for bulk priority assignment, commonly used to set
    /// all preferred libraries to `Priority` level at once.
    ///
    /// # Arguments
    ///
    /// * `libraries` - Vector of library names
    /// * `priority` - Priority level to assign to all libraries
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_multiple_library_priorities(
    ///     vec!["NASA_gas".to_string(), "NASA_cond".to_string()],
    ///     LibraryPriority::Priority
    /// );
    /// ```
    pub fn set_multiple_library_priorities(
        &mut self,
        libraries: Vec<String>,
        priority: LibraryPriority,
    ) {
        for library in libraries {
            let canonical_library = ThermoData::canonical_library_name(&library);
            self.library_priorities
                .insert(canonical_library, priority.clone());
        }
        self.invalidate_lookup_dependent_state();
    }

    /// Sets the same priority for multiple typed library identifiers.
    pub fn set_multiple_library_priorities_ids(
        &mut self,
        libraries: Vec<LibraryId>,
        priority: LibraryPriority,
    ) {
        for library in libraries {
            self.set_library_priority_id(library, priority);
        }
    }

    /// Sets explicit search instructions for direct substance-to-library mapping
    ///
    /// Bypasses the normal priority-based search for specified substances,
    /// forcing them to be searched only in their designated libraries.
    ///
    /// # Arguments
    ///
    /// * `direct_search` - HashMap mapping substance names to specific library names
    ///
    /// # Example
    ///
    /// ```rust
    /// let mut explicit_map = HashMap::new();
    /// explicit_map.insert("H2O".to_string(), "NIST".to_string());
    /// explicit_map.insert("CO2".to_string(), "NASA_gas".to_string());
    /// subs_data.set_explicis_searh_instructions(explicit_map);
    /// ```
    pub fn set_explicit_search_instructions(&mut self, direct_search: HashMap<String, String>) {
        self.explicit_search_instructions = direct_search
            .into_iter()
            .map(|(substance, library)| (substance, ThermoData::canonical_library_name(&library)))
            .collect();
        self.invalidate_lookup_dependent_state();
    }

    /// Legacy spelling retained for compatibility with older callers.
    #[deprecated(note = "Use set_explicit_search_instructions")]
    pub fn set_explicis_searh_instructions(&mut self, direct_search: HashMap<String, String>) {
        self.set_explicit_search_instructions(direct_search);
    }

    /// Sets a single explicit search instruction using typed library identifiers.
    pub fn set_explicit_search_instruction(&mut self, substance: String, library: LibraryId) {
        self.explicit_search_instructions
            .insert(substance, library.canonical_name().to_string());
        self.invalidate_lookup_dependent_state();
    }

    /// Creates the appropriate calculator instance for a given library
    ///
    /// Determines whether a library contains thermodynamic or transport data
    /// and initializes the corresponding calculator type. This is used internally
    /// during the search process to prepare calculators for property calculations.
    ///
    /// # Arguments
    ///
    /// * `library` - Name of the library to create calculator for
    ///
    /// # Returns
    ///
    /// * `Ok(CalculatorType)` - Appropriate calculator for the library
    /// * `Err(SubsDataError::UnsupportedLibrary)` if the library type is not recognized
    ///
    /// # Supported Libraries
    ///
    /// **Thermodynamic**: NASA, NIST, NASA7, NUIG_thermo, Cantera databases
    /// **Transport**: CEA, Aramco_transport
    pub fn create_calculator(&self, library: &str) -> SubsDataResult<CalculatorType> {
        let canonical_library = ThermoData::canonical_library_name(library);
        match ThermoData::library_capability(&canonical_library) {
            Some(LibraryCapability::Thermo) => Ok(CalculatorType::Thermo(create_thermal_by_name(
                &canonical_library,
            ))),
            Some(LibraryCapability::Transport) => Ok(CalculatorType::Transport(
                create_transport_calculator_by_name(&canonical_library)?,
            )),
            None => Err(SubsDataError::UnsupportedLibrary(library.to_string())),
        }
    }

    /// Performs comprehensive search for all substances across configured libraries
    ///
    /// Implements a hierarchical search strategy:
    /// 1. **Priority Libraries**: Searches high-priority databases first
    /// 2. **Permitted Libraries**: Falls back to secondary databases if not found
    /// 3. **Explicit Instructions**: Processes direct substance-to-library mappings
    /// 4. **Result Storage**: Stores all findings with provenance information
    ///
    /// For each found substance, initializes the appropriate calculator and stores
    /// the complete search result including library name, priority type, raw data,
    /// and calculator instance.
    ///
    /// # Returns
    ///
    /// * `Ok(())` - All searches completed successfully
    /// * `Err(SubsDataError)` - Error during search or calculator initialization
    ///
    /// # Example
    ///
    /// ```rust
    /// let mut subs_data = SubsData::new();
    /// subs_data.set_substances(vec!["CO".to_string(), "H2O".to_string()]);
    /// subs_data.set_multiple_library_priorities(
    ///     vec!["NASA_gas".to_string()],
    ///     LibraryPriority::Priority
    /// );
    /// subs_data.search_substances()?;
    /// ```
    pub fn search_substances(&mut self) -> SubsDataResult<()> {
        let previous_search_states = self.search_states.clone();
        let result = (|| -> SubsDataResult<()> {
            // Get all priority libraries
            let mut priority_libs: Vec<String> = self
                .library_priorities
                .iter()
                .filter(move |&(_, &ref p)| *p == LibraryPriority::Priority)
                .map(|(lib, _)| ThermoData::canonical_library_name(lib))
                .collect();
            priority_libs.sort();
            priority_libs.dedup();
            // Get all permitted libraries
            let mut permitted_libs: Vec<String> = self
                .library_priorities
                .iter()
                .filter(move |&(_, &ref p)| *p == LibraryPriority::Permitted)
                .map(|(lib, _)| ThermoData::canonical_library_name(lib))
                .collect();
            permitted_libs.sort();
            permitted_libs.dedup();
            let subs_to_search_explicit_instruction: HashSet<String> =
                self.explicit_search_instructions.keys().cloned().collect();
            let mut subs_to_search = self.substances.clone();
            // Use a set for explicit lookup filtering so repeated membership
            // checks stay linear instead of quadratic as the input grows.
            subs_to_search.retain(|subs_i| !subs_to_search_explicit_instruction.contains(subs_i));
            // Search for each substance
            for substance in subs_to_search {
                let mut found_thermo = false;
                let mut found_transport = false;
                // Try priority libraries first
                for lib in &priority_libs {
                    if let Some(resolved_record) = self.resolve_library_record(lib, &substance)? {
                        let substance_data = &resolved_record.data;
                        // Create the calculator for this property
                        let mut calculator = self.create_calculator(lib)?;

                        match &mut calculator {
                            CalculatorType::Thermo(thermo) => {
                                if found_thermo {
                                    continue;
                                }
                                thermo.from_serde(substance_data.clone())?;
                                self.store_resolved_search_result(
                                    &substance,
                                    WhatIsFound::Thermo,
                                    lib.clone(),
                                    resolved_record.record_key.clone(),
                                    LibraryPriority::Priority,
                                    substance_data.clone(),
                                    calculator.clone(),
                                    false,
                                );
                                found_thermo = true;
                            }
                            CalculatorType::Transport(transport) => {
                                if found_transport {
                                    continue;
                                }
                                transport.from_serde(substance_data.clone())?;
                                self.store_resolved_search_result(
                                    &substance,
                                    WhatIsFound::Transport,
                                    lib.clone(),
                                    resolved_record.record_key.clone(),
                                    LibraryPriority::Priority,
                                    substance_data.clone(),
                                    calculator.clone(),
                                    false,
                                );
                                found_transport = true;
                            }
                        }
                    }
                }

                // If one of the properties is still missing, try permitted libraries.
                if !found_thermo || !found_transport {
                    for lib in &permitted_libs {
                        if let Some(resolved_record) =
                            self.resolve_library_record(lib, &substance)?
                        {
                            let substance_data = &resolved_record.data;
                            let mut calculator = self.create_calculator(lib)?;
                            match &mut calculator {
                                CalculatorType::Thermo(thermo) => {
                                    if found_thermo {
                                        continue;
                                    }
                                    thermo.newinstance()?;
                                    thermo.from_serde(substance_data.clone())?;
                                    self.store_resolved_search_result(
                                        &substance,
                                        WhatIsFound::Thermo,
                                        lib.clone(),
                                        resolved_record.record_key.clone(),
                                        LibraryPriority::Permitted,
                                        substance_data.clone(),
                                        calculator.clone(),
                                        false,
                                    );
                                    found_thermo = true;
                                }
                                CalculatorType::Transport(transport) => {
                                    if found_transport {
                                        continue;
                                    }
                                    //   let _ = transport.newinstance();
                                    transport.from_serde(substance_data.clone())?;
                                    self.store_resolved_search_result(
                                        &substance,
                                        WhatIsFound::Transport,
                                        lib.clone(),
                                        resolved_record.record_key.clone(),
                                        LibraryPriority::Permitted,
                                        substance_data.clone(),
                                        calculator.clone(),
                                        false,
                                    );
                                    found_transport = true;
                                }
                            }
                        }
                    }
                }

                // If neither property was found, mark the substance as missing.
                if !found_thermo && !found_transport {
                    self.insert_not_found_if_absent(&substance);
                }
            } // end of for loop
            // now let's check if we have any explicit search instructions
            let explicit_search_instructions: Vec<(String, String)> = self
                .explicit_search_instructions
                .iter()
                .map(|(substance, library)| {
                    (
                        substance.clone(),
                        ThermoData::canonical_library_name(library),
                    )
                })
                .collect();
            for (substance, library) in explicit_search_instructions {
                let mut found_thermo = false;
                let mut found_transport = false;
                if let Some(resolved_record) = self.resolve_library_record(&library, &substance)? {
                    let substance_data = &resolved_record.data;
                    // Create the calculator for this property
                    let mut calculator = self.create_calculator(&library)?;
                    match &mut calculator {
                        CalculatorType::Thermo(thermo) => {
                            thermo.from_serde(substance_data.clone())?;
                            self.store_resolved_search_result(
                                &substance,
                                WhatIsFound::Thermo,
                                library.clone(),
                                resolved_record.record_key.clone(),
                                LibraryPriority::Explicit,
                                substance_data.clone(),
                                calculator.clone(),
                                true,
                            );
                            found_thermo = true;
                        }
                        CalculatorType::Transport(transport) => {
                            transport.from_serde(substance_data.clone())?;
                            self.store_resolved_search_result(
                                &substance,
                                WhatIsFound::Transport,
                                library.clone(),
                                resolved_record.record_key.clone(),
                                LibraryPriority::Explicit,
                                substance_data.clone(),
                                calculator.clone(),
                                true,
                            );
                            found_transport = true;
                        }
                    }
                }
                // If explicit lookup missed both properties, keep any previous hits intact.
                if !found_thermo && !found_transport && !self.search_states.contains_key(&substance)
                {
                    self.insert_not_found_if_absent(&substance);
                }
            } // for loop
            Ok(())
        })();
        if result.is_err() {
            self.search_states = previous_search_states;
            self.search_results = self.get_all_results();
        }
        result
    }
    /////////////////////////////////////GETTERS////////////////////////////////////
    /// Retrieves a compatibility snapshot of the complete search results for a
    /// specific substance.
    ///
    /// The snapshot is reconstructed from the canonical typed search state.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance to query
    ///
    /// # Returns
    ///
    /// * `Some(HashMap)` - Map of data types to search results
    /// * `None` - Substance was not searched or found
    ///
    /// # Example
    ///
    /// ```rust
    /// if let Some(results) = subs_data.get_substance_result("CO2") {
    ///     if let Some(Some(thermo_result)) = results.get(&WhatIsFound::Thermo) {
    ///         println!(
    ///             "Found CO2 thermodynamic data in: {}",
    ///             thermo_result.library()
    ///         );
    ///     }
    /// }
    /// ```
    pub fn get_substance_result(
        &self,
        substance: &str,
    ) -> Option<HashMap<WhatIsFound, Option<SearchResult>>> {
        self.search_states
            .get(substance)
            .map(SubstanceSearchState::to_compat_map)
    }

    /// Retrieve the canonical typed search state for one substance.
    pub fn get_substance_search_state(&self, substance: &str) -> Option<&SubstanceSearchState> {
        self.search_states.get(substance)
    }

    /// Retrieve the canonical found record for one substance and property.
    ///
    /// This is the direct read-only accessor for canonical lookup data.
    pub fn get_search_result(
        &self,
        substance: &str,
        found_kind: WhatIsFound,
    ) -> Option<&SearchResult> {
        self.get_canonical_search_result(substance, found_kind)
    }

    /// Retrieve the canonical found record for one substance and property.
    ///
    /// This keeps read paths anchored to the typed search state instead of the
    /// legacy compatibility snapshot.
    pub(crate) fn get_canonical_search_result(
        &self,
        substance: &str,
        found_kind: WhatIsFound,
    ) -> Option<&SearchResult> {
        self.search_states.get(substance)?.result(found_kind)
    }

    /// Retrieve the source library name for one canonical found record.
    pub(crate) fn get_canonical_search_library(
        &self,
        substance: &str,
        found_kind: WhatIsFound,
    ) -> SubsDataResult<String> {
        self.get_canonical_search_result(substance, found_kind)
            .map(|result| result.library().to_string())
            .ok_or_else(|| SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: match found_kind {
                    WhatIsFound::Thermo => "Thermo".to_string(),
                    WhatIsFound::Transport => "Transport".to_string(),
                    WhatIsFound::NotFound => "NotFound".to_string(),
                },
            })
    }

    /// Returns a compatibility snapshot of all search results for all substances.
    ///
    /// The returned map is reconstructed from the canonical typed state.
    ///
    /// # Returns
    ///
    /// Reference to the complete search results HashMap
    pub fn get_all_results(&self) -> HashMap<String, HashMap<WhatIsFound, Option<SearchResult>>> {
        self.search_states
            .iter()
            .map(|(substance, state)| (substance.clone(), state.to_compat_map()))
            .collect()
    }

    /// Retrieve all canonical typed search states.
    pub fn get_all_search_states(&self) -> &HashMap<String, SubstanceSearchState> {
        &self.search_states
    }

    /// Returns a mutable canonical search result for the requested property.
    ///
    /// This helper intentionally bypasses the legacy compatibility map.
    pub(crate) fn get_canonical_search_result_mut(
        &mut self,
        substance: &str,
        found_kind: WhatIsFound,
    ) -> SubsDataResult<&mut SearchResult> {
        let state = self
            .search_states
            .get_mut(substance)
            .ok_or_else(|| SubsDataError::SubstanceNotFound(substance.to_string()))?;

        // If the canonical search state says the substance is missing entirely,
        // keep that contract visible to all mutable lookup callers.
        if state.is_missing() {
            return Err(SubsDataError::SubstanceNotFound(substance.to_string()));
        }

        let slot = match found_kind {
            WhatIsFound::Thermo => &mut state.thermo,
            WhatIsFound::Transport => &mut state.transport,
            WhatIsFound::NotFound => {
                return Err(SubsDataError::CalculatorNotAvailable {
                    substance: substance.to_string(),
                    calc_type: "NotFound".to_string(),
                });
            }
        };

        match slot {
            PropertySearchState::Found(result) => Ok(result),
            PropertySearchState::NotSearched
            | PropertySearchState::Missing
            | PropertySearchState::Failed(_) => Err(SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: match found_kind {
                    WhatIsFound::Thermo => "Thermo".to_string(),
                    WhatIsFound::Transport => "Transport".to_string(),
                    WhatIsFound::NotFound => "NotFound".to_string(),
                },
            }),
        }
    }

    /// Returns a list of substances that were not found in any library
    ///
    /// Useful for identifying missing data and determining which substances
    /// need alternative data sources or manual input.
    ///
    /// # Returns
    ///
    /// Vector of substance names that were not found
    ///
    /// # Example
    ///
    /// ```rust
    /// let not_found = subs_data.get_not_found_substances();
    /// if !not_found.is_empty() {
    ///     println!("Missing data for: {:?}", not_found);
    /// }
    /// ```
    pub fn get_not_found_substances(&self) -> Vec<String> {
        self.search_states
            .iter()
            .filter(|(_, state)| state.is_missing())
            .map(|(substance, _)| substance.clone())
            .collect()
    }

    /// Returns substances found in high-priority libraries
    ///
    /// Lists substances that were successfully found in libraries marked with
    /// `LibraryPriority::Priority`, indicating the most trusted data sources.
    ///
    /// # Returns
    ///
    /// Vector of substance names found in priority libraries
    ///
    /// # Example
    ///
    /// ```rust
    /// let priority_found = subs_data.get_priority_found_substances();
    /// println!("High-quality data available for: {:?}", priority_found);
    /// ```
    pub fn get_priority_found_substances(&self) -> Vec<String> {
        self.search_states
            .iter()
            .filter(|(_, state)| state.has_priority_found())
            .map(|(substance, _)| substance.clone())
            .collect()
    }

    /// Search substances by elements and populate search results
    pub fn search_by_elements(&mut self, elements: Vec<String>) -> SubsDataResult<Vec<String>> {
        let found_substances = self.thermo_data.search_by_elements(elements);
        self.substances = found_substances.clone();
        self.populate_element_search_results(found_substances)
    }

    /// Search substances containing only specified elements and populate search results
    pub fn search_by_elements_only(
        &mut self,
        elements: Vec<String>,
    ) -> SubsDataResult<Vec<String>> {
        let found_substances = self.thermo_data.search_by_elements_only(elements);
        self.substances = found_substances.clone();
        self.populate_element_search_results(found_substances)
    }

    /// Common function to populate search results from element-based searches
    fn populate_element_search_results(
        &mut self,
        found_substances: Vec<String>,
    ) -> SubsDataResult<Vec<String>> {
        for substance in &found_substances {
            let library_data_pairs: Vec<(String, Value)> = self
                .thermo_data
                .hashmap_of_thermo_data
                .get(substance)
                .map(|thermo_data| {
                    let mut libraries: Vec<String> = thermo_data.keys().cloned().collect();
                    libraries.sort();
                    libraries
                        .into_iter()
                        .filter_map(|library| {
                            thermo_data
                                .get(&library)
                                .cloned()
                                .map(|data| (library, data))
                        })
                        .collect()
                })
                .unwrap_or_default();

            for (library, data) in library_data_pairs {
                let mut calculator = self.create_calculator(&library)?;

                if let CalculatorType::Thermo(thermo) = &mut calculator {
                    thermo.from_serde(data.clone())?;
                    let search_result = SearchResult {
                        library: library.clone(),
                        record_key: substance.clone(),
                        priority_type: LibraryPriority::Priority,
                        data: data.clone(),
                        calculator: Some(calculator.clone()),
                    };
                    let state = self
                        .search_states
                        .entry(substance.clone())
                        .or_insert_with(SubstanceSearchState::new);
                    state.set_found(WhatIsFound::Thermo, search_result, true);
                    self.sync_search_state_view(substance);
                }
            }
        }
        Ok(found_substances)
    }
    ///
    /// Vector of substance names found in permitted libraries
    ///
    /// # Example
    ///
    /// ```rust
    /// let permitted_found = subs_data.get_permitted_found_substances();
    /// if !permitted_found.is_empty() {
    ///     println!("Secondary data sources used for: {:?}", permitted_found);
    /// }
    /// ```
    pub fn get_permitted_found_substances(&self) -> Vec<String> {
        self.search_states
            .iter()
            .filter(|(_, state)| state.has_permitted_found())
            .map(|(substance, _)| substance.clone())
            .collect()
    }

    /// Retrieves a function closure for calculating thermodynamic properties
    ///
    /// Returns a temperature-dependent function for the specified property type.
    /// These functions are created from polynomial fits or other mathematical models.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance
    /// * `data_type` - Type of property function (e.g., `DataType::Cp_fun`)
    ///
    /// # Returns
    ///
    /// * `Some(function)` - Closure that takes temperature and returns property value
    /// * `None` - Function not available for this substance/property combination
    ///
    /// # Example
    ///
    /// ```rust
    /// if let Some(cp_func) = subs_data.get_thermo_function("CO2", DataType::Cp_fun) {
    ///     let cp_at_500k = cp_func(500.0); // Calculate Cp at 500 K
    /// }
    /// ```
    pub fn get_thermo_function(
        &self,
        substance: &str,
        data_type: DataType,
    ) -> Option<&(dyn Fn(f64) -> f64 + Send + Sync)> {
        self.therm_map_of_fun
            .get(substance)
            .and_then(|map| map.get(&data_type))
            .and_then(|opt| opt.as_deref())
    }

    /// Retrieves a symbolic expression for thermodynamic property calculations
    ///
    /// Returns a mathematical expression that can be manipulated symbolically,
    /// differentiated, integrated, or converted to other forms for analytical work.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance
    /// * `data_type` - Type of property expression (e.g., `DataType::Cp_sym`)
    ///
    /// # Returns
    ///
    /// * `Some(expression)` - Symbolic mathematical expression
    /// * `None` - Expression not available for this substance/property combination
    ///
    /// # Example
    ///
    /// ```rust
    /// if let Some(cp_expr) = subs_data.get_thermo_symbolic("H2O", DataType::Cp_sym) {
    ///     let cp_func = cp_expr.lambdify1D(); // Convert to evaluable function
    ///     let cp_derivative = cp_expr.diff("T"); // Symbolic differentiation
    /// }
    /// ```
    pub fn get_thermo_symbolic(&self, substance: &str, data_type: DataType) -> Option<&Expr> {
        self.therm_map_of_sym
            .get(substance)
            .and_then(|map| map.get(&data_type))
            .and_then(|opt| opt.as_deref())
    }
    ///////////////////////////////////////THERMAL PROPERTIES////////////////////////////////////

    /// Generic calculator access pattern - reduces boilerplate for calculator operations
    pub(crate) fn with_thermo_calculator<F, R>(
        &mut self,
        substance: &str,
        f: F,
    ) -> SubsDataResult<R>
    where
        F: FnOnce(&mut dyn ThermoCalculator) -> Result<R, ThermoError>,
    {
        let search_result = self.get_canonical_search_result_mut(substance, WhatIsFound::Thermo)?;
        let calculator = search_result.calculator_mut().ok_or_else(|| {
            SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: "Thermo".to_string(),
            }
        })?;

        match calculator {
            CalculatorType::Thermo(thermo) => {
                let result: SubsDataResult<R> = f(thermo).map_err(Into::into);
                if let Err(ref e) = result {
                    crate::log_error!(self.logger, e, substance, "with_thermo_calculator");
                }
                result
            }
            CalculatorType::Transport(_) => Err(SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: "Thermo (found Transport instead)".to_string(),
            }),
        }
    }

    /// Extracts polynomial coefficients for thermodynamic calculations at a specific temperature
    ///
    /// Determines which temperature range applies and extracts the appropriate polynomial
    /// coefficients from the substance's thermodynamic data. This is typically required
    /// before performing property calculations.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance
    /// * `temperature` - Temperature [K] to determine coefficient range
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Coefficients extracted successfully
    /// * `Err(SubsDataError)` - Invalid temperature or substance not found
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.extract_thermal_coeffs("CO2", 500.0)?;
    /// ```
    pub fn extract_thermal_coeffs(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<()> {
        if temperature <= 0.0 {
            let error = SubsDataError::InvalidTemperature(temperature);
            return self.log_and_propagate(Err(error), substance, "extract_thermal_coeffs");
        }
        let result = self.with_thermo_calculator(substance, |thermo| {
            thermo.extract_model_coefficients(temperature)
        });
        self.log_and_propagate(result, substance, "extract_thermal_coeffs")
    }

    /// Extracts polynomial coefficients for all substances at a specific temperature
    ///
    /// Batch operation that extracts thermodynamic polynomial coefficients for all
    /// substances in the substance list. Useful for preparing multiple substances
    /// for property calculations at the same temperature.
    ///
    /// # Arguments
    ///
    /// * `temperature` - Temperature [K] to determine coefficient ranges
    ///
    /// # Returns
    ///
    /// * `Ok(())` - All coefficients extracted successfully
    /// * `Err(SubsDataError)` - Invalid temperature or extraction failure
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.extract_all_thermal_coeffs(298.15)?; // Standard conditions
    /// ```
    pub fn extract_all_thermal_coeffs(&mut self, temperature: f64) -> SubsDataResult<()> {
        if temperature <= 0.0 {
            let error = SubsDataError::InvalidTemperature(temperature);
            crate::log_error!(
                self.logger,
                &error,
                "ALL_SUBSTANCES",
                "extract_all_thermal_coeffs"
            );
            return Err(error);
        }
        for substance in self.substances.clone() {
            let result = self.extract_thermal_coeffs(&substance, temperature);
            self.log_and_propagate(result, &substance, "extract_all_thermal_coeffs")?;
        }
        Ok(())
    }
    /// we have right to calculate thermal properties, construct functions and symbolc expressions
    /// only if the thermal polynomial coefficients corresponds to the given temperature
    /// so this method if the current coefficients for this substance are suitable for the given temperature
    pub fn is_coeffs_valid_for_T(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<bool> {
        if temperature <= 0.0 {
            let error = SubsDataError::InvalidTemperature(temperature);
            return self.log_and_propagate(Err(error), substance, "is_coeffs_valid_for_T");
        }
        let result = self.with_thermo_calculator(substance, |thermo| {
            thermo.is_coeffs_valid_for_T(temperature)
        });
        self.log_and_propagate(result, substance, "is_coeffs_valid_for_T")
    }
    /// we have right to calculate thermal properties, construct functions and symbolc expressions
    /// only if the thermal polynomial coefficients corresponds to the given temperature
    /// so this method checks if the current coefficients for this substance are suitable for the given temperature
    /// and if no the coeddicients are renewed
    pub fn extract_coeffs_if_current_coeffs_not_valid(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<bool> {
        let check_flag = self.is_coeffs_valid_for_T(substance, temperature);
        match check_flag {
            Ok(flag) => match flag {
                true => return Ok(true),
                false => {
                    self.extract_thermal_coeffs(substance, temperature)?;
                    return self.log_and_propagate(
                        check_flag,
                        substance,
                        "extract_coeffs_if_current_coeffs_not_valid",
                    );
                }
            },
            Err(error) => {
                crate::log_error!(
                    self.logger,
                    &error,
                    "ALL_SUBSTANCES",
                    "extract_coeffs_if_current_coeffs_not_valid"
                );
                return Err(error);
            }
        }
    }
    /// we have right to calculate thermal properties, construct functions and symbolc expressions
    /// only if the thermal polynomial coefficients corresponds to the given temperature
    /// so this method checks if the current coefficients for this substance are suitable for the given temperature
    /// and if no the coeddicients are renewed for all substances
    pub fn extract_coeffs_if_current_coeffs_not_valid_for_all_subs(
        &mut self,
        temperature: f64,
    ) -> SubsDataResult<Vec<String>> {
        let mut subs_with_changed_coeffs = Vec::new();
        for substance in self.substances.clone() {
            let flag = self.extract_coeffs_if_current_coeffs_not_valid(&substance, temperature)?;
            if flag {
                subs_with_changed_coeffs.push(substance.clone());
            }
        }
        Ok(subs_with_changed_coeffs)
    }

    /// Calculates )
    /// Sets the temperature range for thermodynamic calculations for a specific substance
    ///
    /// Defines the valid temperature interval for property calculations. This affects
    /// which polynomial coefficients are used and can enable temperature range validation.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance
    /// * `T_min` - Minimum temperature [K]
    /// * `T_max` - Maximum temperature [K]
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Temperature range set successfully
    /// * `Err(SubsDataError)` - Invalid range or substance not found
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T_range_for_thermo("H2O", 273.15, 2000.0)?;
    /// ```
    pub fn set_T_range_for_thermo(
        &mut self,
        substance: &str,
        T_min: f64,
        T_max: f64,
    ) -> SubsDataResult<()> {
        let result =
            self.with_thermo_calculator(substance, |thermo| thermo.set_T_interval(T_min, T_max));
        self.log_and_propagate(result, substance, "set_T_range_for_thermo")
    }
    /// Sets the same temperature range for all substances' thermodynamic calculations
    ///
    /// Batch operation that applies the same temperature limits to all substances.
    /// Useful when all substances will be used within the same temperature range.
    ///
    /// # Arguments
    ///
    /// * `T_min` - Minimum temperature [K] for all substances
    /// * `T_max` - Maximum temperature [K] for all substances
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Temperature ranges set for all substances
    /// * `Err(SubsDataError)` - Invalid range or setting failure
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T_range_for_all_thermo(300.0, 1500.0)?; // Combustion range
    /// ```
    pub fn set_T_range_for_all_thermo(&mut self, T_min: f64, T_max: f64) -> SubsDataResult<()> {
        for substance in self.substances.clone() {
            let result = self.set_T_range_for_thermo(&substance, T_min, T_max);
            self.log_and_propagate(result, &substance, "set_T_range_for_all_thermo")?;
        }
        Ok(())
    }
    /// Parses and validates thermodynamic polynomial coefficients for a substance
    ///
    /// Processes the raw coefficient data from the library, validates the polynomial
    /// structure, and prepares the coefficients for property calculations.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance to parse coefficients for
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Coefficients parsed successfully
    /// * `Err(SubsDataError)` - Parsing failure or substance not found
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.parse_thermal_coeffs("CH4")?;
    /// ```
    pub fn parse_thermal_coeffs(&mut self, substance: &str) -> SubsDataResult<()> {
        let result = self.with_thermo_calculator(substance, |thermo| thermo.parse_coefficients());
        self.log_and_propagate(result, substance, "parse_thermal_coeffs")
    }

    /// Parses thermodynamic coefficients for all substances in the substance list
    ///
    /// Batch operation that processes polynomial coefficients for all substances,
    /// preparing them for property calculations. This is typically done once after
    /// the search process is complete.
    ///
    /// # Returns
    ///
    /// * `Ok(())` - All coefficients parsed successfully
    /// * `Err(SubsDataError)` - Parsing failure for any substance
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.search_substances()?;
    /// subs_data.parse_all_thermal_coeffs()?;
    /// ```
    pub fn parse_all_thermal_coeffs(&mut self) -> SubsDataResult<()> {
        for substance in self.substances.clone() {
            let result = self.parse_thermal_coeffs(&substance);
            self.log_and_propagate(result, &substance, "parse_all_thermal_coeffs")?;
        }
        Ok(())
    }

    /// Fits polynomial coefficients to the specified temperature interval for a substance
    ///
    /// Adjusts or refits the thermodynamic polynomial coefficients to optimize accuracy
    /// within the specified temperature range. This can improve calculation precision
    /// for specific temperature intervals.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance to fit coefficients for
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Coefficients fitted successfully
    /// * `Err(SubsDataError)` - Fitting failure or substance not found
    ///
    /// # Note
    ///
    /// Temperature interval must be set using `set_T_range_for_thermo` before calling this method.
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T_range_for_thermo("O2", 500.0, 1000.0)?;
    /// subs_data.fitting_thermal_coeffs_for_T_interval("O2")?;
    /// ```
    pub fn fitting_thermal_coeffs_for_T_interval(&mut self, substance: &str) -> SubsDataResult<()> {
        let result =
            self.with_thermo_calculator(substance, |thermo| thermo.fitting_coeffs_for_T_interval());
        self.log_and_propagate(result, substance, "fitting_thermal_coeffs_for_T_interval")
    }

    /// Fits polynomial coefficients for all substances to their specified temperature intervals
    ///
    /// Batch operation that optimizes thermodynamic polynomial coefficients for all substances
    /// within their respective temperature ranges. Improves calculation accuracy across
    /// the entire substance set.
    ///
    /// # Returns
    ///
    /// * `Ok(())` - All coefficients fitted successfully
    /// * `Err(SubsDataError)` - Fitting failure for any substance
    ///
    /// # Note
    ///
    /// Temperature intervals must be set for all substances before calling this method.
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T_range_for_all_thermo(400.0, 1200.0)?;
    /// subs_data.fitting_all_thermal_coeffs_for_T_interval()?;
    /// ```
    pub fn fitting_all_thermal_coeffs_for_T_interval(&mut self) -> SubsDataResult<()> {
        for substance in self.substances.clone() {
            let result = self.fitting_thermal_coeffs_for_T_interval(&substance);
            self.log_and_propagate(
                result,
                &substance,
                "fitting_all_thermal_coeffs_for_T_interval",
            )?;
        }
        Ok(())
    }

    /// Calculates integral mean values for thermodynamic properties over temperature range
    ///
    /// Computes temperature-averaged values of thermodynamic properties (Cp, dH, dS)
    /// over the specified temperature interval. Useful for mean property calculations
    /// in process design and analysis.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance to calculate integral means for
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Integral means calculated successfully
    /// * `Err(SubsDataError)` - Calculation failure or substance not found
    ///
    /// # Note
    ///
    /// Temperature interval must be set before calling this method.
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T_range_for_thermo("N2", 300.0, 800.0)?;
    /// subs_data.integr_mean("N2")?;
    /// ```
    pub fn integr_mean(&mut self, substance: &str) -> SubsDataResult<()> {
        let result = self.with_thermo_calculator(substance, |thermo| thermo.integr_mean());
        self.log_and_propagate(result, substance, "integr_mean")
    }

    /// Calculates thermodynamic properties (Cp, dH, dS) for a substance at given temperature
    ///
    /// Computes heat capacity, enthalpy change, and entropy change for the specified
    /// substance at the given temperature using the appropriate thermodynamic model.
    ///
    /// # Arguments
    ///
    /// * `substance` - Name of the substance
    /// * `temperature` - Temperature [K] for property calculation
    ///
    /// # Returns
    ///
    /// * `Ok((cp, dh, ds))` - Tuple of (heat capacity [J/mol·K], enthalpy change [J/mol], entropy change [J/mol·K])
    /// * `Err(SubsDataError)` - Calculation failure, invalid temperature, or substance not found
    ///
    /// # Example
    ///
    /// ```rust
    /// let (cp, dh, ds) = subs_data.calculate_thermo_properties("CO2", 500.0)?;
    /// println!("Cp: {:.2} J/mol·K, dH: {:.2} J/mol, dS: {:.2} J/mol·K", cp, dh, ds);
    /// ```
    pub fn calculate_thermo_properties(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<(f64, f64, f64)> {
        let result = self._calculate_thermo_properties_internal(substance, temperature);
        self.log_and_propagate(result, substance, "calculate_thermo_properties")
    }

    fn _calculate_thermo_properties_internal(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<(f64, f64, f64)> {
        if temperature <= 0.0 {
            return Err(SubsDataError::InvalidTemperature(temperature));
        }

        let state = self
            .get_substance_search_state(substance)
            .ok_or_else(|| SubsDataError::SubstanceNotFound(substance.to_string()))?;

        // Check if the canonical typed state explicitly marks this substance as missing.
        if state.is_missing() {
            return Err(SubsDataError::SubstanceNotFound(substance.to_string()));
        }

        let search_result = self
            .get_canonical_search_result(substance, WhatIsFound::Thermo)
            .ok_or_else(|| SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: "Thermo".to_string(),
            })?;

        match search_result {
            SearchResult {
                calculator: Some(CalculatorType::Thermo(thermo)),
                ..
            } => {
                let mut thermo = thermo.clone();

                thermo.calculate_Cp_dH_dS(temperature)?;

                let cp = thermo.get_Cp()?;
                let dh = thermo.get_dh()?;
                let ds = thermo.get_ds()?;

                Self::validate_finite("Cp", cp)?;
                Self::validate_finite("dH", dh)?;
                Self::validate_finite("dS", ds)?;

                Ok((cp, dh, ds))
            }
            SearchResult {
                calculator: Some(CalculatorType::Transport(_)),
                ..
            } => Err(SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: "Thermo (found Transport instead)".to_string(),
            }),
            SearchResult {
                calculator: None, ..
            } => Err(SubsDataError::CalculatorNotAvailable {
                substance: substance.to_string(),
                calc_type: "Thermo".to_string(),
            }),
        }
    }

    /// Sets molar masses for substances with optional unit specification
    ///
    /// Assigns molar mass values to substances, which are used in density calculations,
    /// unit conversions, and other mass-based property calculations.
    ///
    /// # Arguments
    ///
    /// * `M_map` - HashMap mapping substance names to their molar masses
    /// * `M_unit` - Optional unit specification (e.g., "g/mol", "kg/mol")
    ///
    /// # Example
    ///
    /// ```rust
    /// let mut molar_masses = HashMap::new();
    /// molar_masses.insert("CO2".to_string(), 44.01);
    /// molar_masses.insert("H2O".to_string(), 18.015);
    /// subs_data.set_M(molar_masses, Some("g/mol".to_string()));
    /// ```
    fn validate_finite(field: &str, value: f64) -> SubsDataResult<()> {
        if value.is_finite() {
            Ok(())
        } else {
            Err(SubsDataError::InvalidPhysicalValue {
                field: field.to_string(),
                value,
            })
        }
    }

    fn validate_finite_positive(field: &str, value: f64) -> SubsDataResult<()> {
        if value.is_finite() && value > 0.0 {
            Ok(())
        } else {
            Err(SubsDataError::InvalidPhysicalValue {
                field: field.to_string(),
                value,
            })
        }
    }

    pub fn set_M(
        &mut self,
        M_map: HashMap<String, f64>,
        M_unit: Option<String>,
    ) -> SubsDataResult<()> {
        for (substance, molar_mass) in &M_map {
            Self::validate_finite_positive(&format!("molar_mass:{}", substance), *molar_mass)?;
        }
        self.molar_mass_by_substance = M_map;
        self.Molar_mass_unit = M_unit;
        self.invalidate_pressure_and_molar_mass_dependent_state();
        Ok(())
    }
    /// Sets system pressure with optional unit specification
    ///
    /// Defines the pressure for property calculations. This affects density calculations,
    /// phase equilibrium, and pressure-dependent transport properties.
    ///
    /// # Arguments
    ///
    /// * `P` - Pressure value
    /// * `P_unit` - Optional unit specification (e.g., "Pa", "atm", "bar")
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_P(101325.0, Some("Pa".to_string())); // Standard atmospheric pressure
    /// subs_data.set_P(1.0, Some("atm".to_string()));      // Alternative specification
    /// ```
    pub fn set_P(&mut self, P: f64, P_unit: Option<String>) -> SubsDataResult<()> {
        Self::validate_finite_positive("pressure", P)?;
        self.P = Some(P);
        self.P_unit = P_unit;
        self.invalidate_pressure_and_molar_mass_dependent_state();
        Ok(())
    }

    /// Sets system temperature for property calculations
    ///
    /// Defines the default temperature for property calculations. This is used
    /// when temperature is not explicitly specified in calculation methods.
    ///
    /// # Arguments
    ///
    /// * `T` - Temperature [K]
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.set_T(298.15); // Standard temperature (25°C)
    /// ```
    pub fn set_T(&mut self, T: f64) -> SubsDataResult<()> {
        Self::validate_finite_positive("temperature", T)?;

        // Avoid invalidating derived caches when the user re-applies the same
        // temperature value. This keeps repeated UI updates and roundtrips cheap.
        if self.T == Some(T) {
            return Ok(());
        }

        self.T = Some(T);
        self.invalidate_temperature_dependent_state();
        Ok(())
    }

    /// Generic utility for building property maps with error handling
    ///
    /// Internal utility function that reduces boilerplate code for property calculations
    /// across multiple substances. Handles errors gracefully by inserting empty values
    /// for failed calculations while continuing with successful ones.
    ///
    /// # Type Parameters
    ///
    /// * `T` - Type of property values (f64, function closures, expressions)
    /// * `F` - Calculator function type
    ///
    /// # Arguments
    ///
    /// * `calculator_fn` - Function that calculates properties for a single substance
    /// * `target_map` - Map to store calculated properties
    /// * `property_types` - Types of properties to calculate
    /// * `empty_value_fn` - Function to generate empty values for failed calculations
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Property map built successfully
    /// * `Err(SubsDataError)` - The first substance-level failure encountered
    pub fn build_property_map<T, F>(
        &mut self,
        calculator_fn: F,
        target_map: &mut HashMap<String, HashMap<DataType, Option<T>>>,
        _property_types: &[DataType],
        _empty_value_fn: fn() -> Option<T>,
    ) -> SubsDataResult<()>
    where
        F: Fn(&mut Self, &str) -> SubsDataResult<HashMap<DataType, Option<T>>>,
    {
        let mut rebuilt_map = HashMap::new();
        for substance in self.substances.clone() {
            let properties = calculator_fn(self, &substance)?;
            rebuilt_map.insert(substance, properties);
        }
        *target_map = rebuilt_map;
        Ok(())
    }

    /// Helper function to calculate thermo properties for a single substance
    fn calculate_single_thermo_properties(
        &mut self,
        substance: &str,
        temperature: f64,
    ) -> SubsDataResult<HashMap<DataType, Option<f64>>> {
        let (cp, dh, ds) = self._calculate_thermo_properties_internal(substance, temperature)?;
        let mut property_map = HashMap::new();
        property_map.insert(DataType::Cp, Some(cp));
        property_map.insert(DataType::dH, Some(dh));
        property_map.insert(DataType::dS, Some(ds));
        Ok(property_map)
    }

    /// Calculates and stores thermodynamic property values for all substances at specified temperature
    ///
    /// Batch calculation that computes Cp, dH, and dS for all substances at the given temperature
    /// and stores the results in `therm_map_of_properties_values`. This is the primary method
    /// for obtaining numerical property values.
    ///
    /// # Arguments
    ///
    /// * `temperature` - Temperature [K] for property calculations
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Properties calculated and stored successfully
    /// * `Err(SubsDataError)` - Invalid temperature or calculation failure
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.search_substances()?;
    /// subs_data.calculate_therm_map_of_properties(500.0)?;
    ///
    /// // Access calculated values
    /// if let Some(co2_props) = subs_data.therm_map_of_properties_values.get("CO2") {
    ///     if let Some(Some(cp)) = co2_props.get(&DataType::Cp) {
    ///         println!("CO2 heat capacity at 500K: {:.2} J/mol·K", cp);
    ///     }
    /// }
    /// ```
    pub fn calculate_therm_map_of_properties(&mut self, temperature: f64) -> SubsDataResult<()> {
        let mut temp_map = HashMap::new();
        for substance in self.substances.clone() {
            let result = self.calculate_single_thermo_properties(&substance, temperature);
            let properties =
                self.log_and_propagate(result, &substance, "calculate_therm_map_of_properties")?;
            temp_map.insert(substance, properties);
        }
        self.therm_map_of_properties_values = temp_map;
        Ok(())
    }

    /// Helper function to calculate thermo functions for a single substance
    fn calculate_single_thermo_functions(
        &mut self,
        substance: &str,
    ) -> SubsDataResult<HashMap<DataType, Option<Box<dyn Fn(f64) -> f64 + Send + Sync>>>> {
        self.with_thermo_calculator(substance, |thermo| {
            thermo.create_closures_Cp_dH_dS()?;
            let mut function_map = HashMap::new();

            function_map.insert(DataType::Cp_fun, thermo.get_C_fun().ok());
            function_map.insert(DataType::dH_fun, thermo.get_dh_fun().ok());
            function_map.insert(DataType::dS_fun, thermo.get_ds_fun().ok());

            Ok(function_map)
        })
    }

    /// Creates and stores function closures for thermodynamic property calculations
    ///
    /// Generates temperature-dependent function closures (Cp(T), dH(T), dS(T)) for all substances
    /// and stores them in `therm_map_of_fun`. These functions can be used for efficient
    /// property calculations at multiple temperatures.
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Function closures created and stored successfully
    /// * `Err(SubsDataError)` - Function creation failure
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.search_substances()?;
    /// subs_data.calculate_therm_map_of_fun()?;
    ///
    /// // Use generated functions
    /// if let Some(cp_func) = subs_data.get_thermo_function("H2O", DataType::Cp_fun) {
    ///     let cp_300k = cp_func(300.0);
    ///     let cp_400k = cp_func(400.0);
    ///     println!("H2O Cp: {:.2} at 300K, {:.2} at 400K", cp_300k, cp_400k);
    /// }
    /// ```
    pub fn calculate_therm_map_of_fun(&mut self) -> SubsDataResult<()> {
        let mut temp_map = HashMap::new();
        for substance in self.substances.clone() {
            let result = self.calculate_single_thermo_functions(&substance);
            let properties =
                self.log_and_propagate(result, &substance, "calculate_therm_map_of_fun")?;
            temp_map.insert(substance, properties);
        }
        self.therm_map_of_fun = temp_map;
        Ok(())
    }

    /// Helper function to calculate thermo symbolic expressions for a single substance
    fn calculate_single_thermo_symbolic(
        &mut self,
        substance: &str,
    ) -> SubsDataResult<HashMap<DataType, Option<Box<Expr>>>> {
        self.with_thermo_calculator(substance, |thermo| {
            thermo.create_sym_Cp_dH_dS()?;
            let mut sym_map = HashMap::new();

            sym_map.insert(DataType::Cp_sym, thermo.get_Cp_sym().ok().map(Box::new));
            sym_map.insert(DataType::dH_sym, thermo.get_dh_sym().ok().map(Box::new));
            sym_map.insert(DataType::dS_sym, thermo.get_ds_sym().ok().map(Box::new));

            Ok(sym_map)
        })
    }

    /// Creates and stores symbolic expressions for thermodynamic property calculations
    ///
    /// Generates symbolic mathematical expressions for thermodynamic properties (Cp, dH, dS)
    /// and stores them in `therm_map_of_sym`. These expressions can be manipulated symbolically,
    /// differentiated, integrated, or used for analytical work.
    ///
    /// # Returns
    ///
    /// * `Ok(())` - Symbolic expressions created and stored successfully
    /// * `Err(SubsDataError)` - Expression creation failure
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.search_substances()?;
    /// subs_data.calculate_therm_map_of_sym()?;
    ///
    /// // Use symbolic expressions
    /// if let Some(cp_expr) = subs_data.get_thermo_symbolic("CH4", DataType::Cp_sym) {
    ///     let cp_func = cp_expr.lambdify1D();           // Convert to function
    ///     let cp_derivative = cp_expr.diff("T");        // Symbolic differentiation
    ///     let cp_integral = cp_expr.integrate("T");     // Symbolic integration
    /// }
    /// ```
    pub fn calculate_therm_map_of_sym(&mut self) -> SubsDataResult<()> {
        let mut temp_map = HashMap::new();
        for substance in self.substances.clone() {
            let result = self.calculate_single_thermo_symbolic(&substance);
            let properties =
                self.log_and_propagate(result, &substance, "calculate_therm_map_of_sym")?;
            temp_map.insert(substance, properties);
        }
        self.therm_map_of_sym = temp_map;
        Ok(())
    }

    /// Internal helper method for error logging and propagation
    ///
    /// Logs errors to the internal logger while preserving the original error for propagation.
    /// This enables comprehensive error tracking while maintaining normal error handling flow.
    ///
    /// # Type Parameters
    ///
    /// * `T` - Type of the result value
    ///
    /// # Arguments
    ///
    /// * `result` - Result to log and propagate
    /// * `substance` - Name of substance associated with the operation
    /// * `function` - Name of function where error occurred
    ///
    /// # Returns
    ///
    /// The original result, unchanged, after logging any errors
    pub fn log_and_propagate<T>(
        &mut self,
        result: SubsDataResult<T>,
        substance: &str,
        function: &str,
    ) -> SubsDataResult<T> {
        if let Err(ref e) = result {
            crate::log_error!(self.logger, e, substance, function);
        }
        result
    }

    /// Prints all logged errors to the console
    ///
    /// Displays a comprehensive list of all errors that occurred during substance
    /// processing, including substance names, function names, and error details.
    /// Useful for debugging and error analysis.
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.search_substances().ok(); // May generate errors
    /// subs_data.calculate_therm_map_of_properties(500.0).ok();
    /// subs_data.print_error_logs(); // Display any errors that occurred
    /// ```
    pub fn print_error_logs(&self) {
        use crate::Thermodynamics::User_substances_error::ExceptionLogger;
        self.logger.print_logs();
    }

    /// Clears all logged errors from the internal logger
    ///
    /// Removes all previously logged errors, providing a clean slate for subsequent
    /// operations. Useful when reusing the same `SubsData` instance for multiple
    /// calculation sessions.
    ///
    /// # Example
    ///
    /// ```rust
    /// subs_data.print_error_logs();  // Show current errors
    /// subs_data.clear_error_logs();  // Clear error history
    /// // Perform new calculations with clean error log
    /// ```
    pub fn clear_error_logs(&mut self) {
        use crate::Thermodynamics::User_substances_error::ExceptionLogger;
        self.logger.clear_logs();
    }
}