lcdm-core 0.1.0

// lcdm-core/src/spec.rs

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

/// High-level specification for the neutrino model.
///
/// This enum is designed to prevent double-counting between massless (radiation-like)
/// and massive (non-cold dark matter) neutrino components.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub enum NeutrinoSpec {
    /// Effective model (recommended):
    ///
    /// - `masses_ev` enumerates the explicit species being modeled (may include 0.0).
    /// - `n_eff` rescales the thermal bath so the total relativistic energy density matches
    ///   the desired effective number of species.
    ///
    /// In this model, the neutrino background is represented by a single effective temperature
    /// derived from `n_eff` and the number of explicit species.
    Effective {
        /// Effective number of relativistic species (e.g., 3.046).
        n_eff: f64,
        /// Neutrino masses in eV (length N; may include 0.0 for massless species).
        masses_ev: Vec<f64>,
    },

    /// Split model (CLASS-style):
    ///
    /// Massless (ultra-relativistic) and massive (ncdm) components are strictly separated:
    /// - `n_ur` represents the massless-only contribution.
    /// - `masses_ev` represents massive species only (should typically be strictly m > 0).
    /// - `temp_factors` allows per-species temperature scaling relative to a base temperature.
    Split {
        /// Effective number of massless (UR) species (e.g., 3.046 if no massive species).
        n_ur: f64,
        /// Masses (eV) for the massive (ncdm) species (length N_ncdm).
        masses_ev: Vec<f64>,
        /// Temperature ratios T_i / T_base for each massive species (typically all 1.0).
        temp_factors: Vec<f64>,
    },
}

/// Accuracy presets for numerical routines (table sizes, integration steps, etc.).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub enum AccuracyPreset {
    /// Prioritize speed (smaller tables / fewer integration steps).
    Fast,
    /// Balanced default settings.
    #[default]
    Default,
    /// Higher precision (larger tables / more integration steps).
    Precise,
    /// Research/paper-grade settings (very large tables / many integration steps).
    Paper,
}

impl AccuracyPreset {
    /// Number of tabulation points for the neutrino background table.
    pub fn neutrino_n_points(&self) -> usize {
        match self {
            Self::Fast => 256,
            Self::Default => 512,
            Self::Precise => 1024,
            Self::Paper => 2048,
        }
    }

    /// Composite Simpson integration step count for neutrino integrals.
    pub fn neutrino_integ_steps(&self) -> usize {
        match self {
            Self::Fast => 600,
            Self::Default => 1200,
            Self::Precise => 2000,
            Self::Paper => 4000,
        }
    }

    /// Upper cutoff for the dimensionless momentum integral (y_max).
    pub fn neutrino_y_max(&self) -> f64 {
        match self {
            Self::Fast => 35.0_f64,
            Self::Default => 45.0_f64,
            Self::Precise => 50.0_f64,
            Self::Paper => 60.0_f64,
        }
    }
}

/// Single source of truth for cosmological composition at z = 0.
///
/// This struct contains only **physical densities** (ω ≡ Ω h²) and derived
/// temperatures. Any hot-path computations (E(z), distances, growth, etc.)
/// should use values from here to avoid divergence across components.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub struct DerivedParams {
    // --- Basic echo ---
    /// Dimensionless Hubble parameter h.
    pub h: f64,
    /// h² (cached).
    pub h2: f64,
    /// CMB temperature today in Kelvin.
    pub t_cmb0_k: f64,

    // --- Neutrino temperatures ---
    /// Base neutrino temperature:
    ///   T_ν0,base = (4/11)^(1/3) * T_cmb0
    ///
    /// This is the background temperature used for the Split model.
    pub t_nu0_base_k: f64,
    /// Effective neutrino temperature used for the Effective model.
    ///
    /// This scales the bath by (N_eff / N_species)^(1/4) so the total relativistic
    /// energy density matches `n_eff`.
    pub t_nu0_k: f64,

    // --- Physical densities (z = 0) ---
    /// Physical baryon density ω_b0.
    pub omega_b0: f64,
    /// Physical cold dark matter density ω_cdm0.
    pub omega_cdm0: f64,
    /// Total matter physical density ω_m0 = ω_b0 + ω_cdm0.
    pub omega_m0: f64,

    /// Physical curvature density ω_k0.
    pub omega_k0: f64,
    /// Physical dark energy density ω_de0 (Λ in ΛCDM).
    pub omega_de0: f64,

    /// Physical photon density ω_γ0.
    pub omega_g0: f64,

    // --- Neutrino components (diagnostics / split info) ---
    /// Massless-only (UR) contribution used in the Split model; 0.0 in Effective mode.
    pub omega_ur0: f64,
    /// Total neutrino physical density ω_ν0 (massive + any species included via the model).
    pub omega_nu0: f64,

    // --- Sums ---
    /// Sum of present-day physical densities for photons + UR + neutrinos.
    ///
    /// Note: this does not imply strict (1+z)^4 scaling for the massive part at all z.
    pub omega_gnu0: f64,

    /// Sum of all non-dark-energy components today:
    ///   ω_non_de0 = ω_γ0 + ω_ur0 + ω_ν0 + ω_m0 + ω_k0
    pub omega_non_de0: f64,
}

impl DerivedParams {
    /// Construct `DerivedParams` and precompute cached sums.
    pub fn new(
        h: f64,
        t_cmb0: f64,
        omega_b: f64,
        omega_cdm: f64,
        omega_k: f64,
        omega_de: f64,
        omega_g: f64,
        omega_ur: f64,
        omega_nu: f64,
        t_nu0_base: f64,
        t_nu0: f64,
    ) -> Self {
        let h2 = h * h;
        let omega_m = omega_b + omega_cdm;
        let omega_non_de = omega_g + omega_ur + omega_nu + omega_m + omega_k;

        Self {
            h,
            h2,
            t_cmb0_k: t_cmb0,
            t_nu0_base_k: t_nu0_base,
            t_nu0_k: t_nu0,
            omega_b0: omega_b,
            omega_cdm0: omega_cdm,
            omega_m0: omega_m,
            omega_k0: omega_k,
            omega_de0: omega_de,
            omega_g0: omega_g,
            omega_ur0: omega_ur,
            omega_nu0: omega_nu,
            omega_gnu0: omega_g + omega_ur + omega_nu,
            omega_non_de0: omega_non_de,
        }
    }
}