andam 0.1.2

A comprehensive Rust library for cosmological calculations, visualizations and simulations
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

//! Dark energy models beyond cosmological constant

use crate::dynamics::Universe;

/// Dark energy equation of state
#[derive(Debug, Clone)]
pub enum DarkEnergyModel {
    /// Cosmological constant: w = -1
    Lambda,
    /// Constant w
    ConstantW(f64),
    /// CPL parametrization: w(a) = w_0 + w_a(1-a)
    CPL { w_0: f64, w_a: f64 },
    /// Early dark energy
    EarlyDE { w_0: f64, omega_e: f64 },
}

impl DarkEnergyModel {
    /// Equation of state at scale factor a
    pub fn w(&self, a: f64) -> f64 {
        match self {
            DarkEnergyModel::Lambda => -1.0,
            DarkEnergyModel::ConstantW(w) => *w,
            DarkEnergyModel::CPL { w_0, w_a } => {
                w_0 + w_a * (1.0 - a)
            },
            DarkEnergyModel::EarlyDE { w_0, .. } => *w_0,
        }
    }

    /// Dark energy density evolution
    pub fn rho_de(&self, a: f64, omega_de_0: f64) -> f64 {
        match self {
            DarkEnergyModel::Lambda => omega_de_0,
            DarkEnergyModel::ConstantW(w) => {
                omega_de_0 * a.powf(-3.0 * (1.0 + w))
            },
            DarkEnergyModel::CPL { w_0, w_a } => {
                // ρ_DE ∝ exp(-3∫(1+w(a'))da'/a')
                let integral = (1.0 + w_0) * a.ln() + w_a * (a - 1.0);
                omega_de_0 * (-3.0 * integral).exp()
            },
            DarkEnergyModel::EarlyDE { w_0, .. } => {
                // Simplified early DE
                omega_de_0 * a.powf(-3.0 * (1.0 + w_0))
            },
        }
    }

    /// Modified Hubble parameter
    pub fn hubble_modified(&self, a: f64, _universe: &Universe) -> f64 {
        let omega_m = 0.3; // Get from universe
        let omega_de_0 = 0.7;

        let rho_m = omega_m * a.powf(-3.0);
        let rho_de = self.rho_de(a, omega_de_0);

        (rho_m + rho_de).sqrt()
    }
}

/// Compare models
pub fn compare_dark_energy_models(
    models: Vec<(&str, DarkEnergyModel)>,
    universe: &Universe,
) -> Vec<(String, Vec<(f64, f64)>)> {
    let mut results = Vec::new();

    for (name, model) in models {
        let mut data = Vec::new();

        for i in 0..100 {
            let a = 0.1 + (i as f64) * 0.009;
            let h = model.hubble_modified(a, universe);
            data.push((a, h));
        }

        results.push((name.to_string(), data));
    }

    results
}

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

    #[test]
    fn test_lambda_cdm() {
        let model = DarkEnergyModel::Lambda;
        assert_eq!(model.w(0.5), -1.0);
        assert_eq!(model.rho_de(0.5, 0.7), 0.7);
    }

    #[test]
    fn test_constant_w() {
        let model = DarkEnergyModel::ConstantW(-0.9);
        assert_eq!(model.w(0.5), -0.9);

        // Check density evolution
        let rho = model.rho_de(0.5, 1.0);
        assert!(rho > 0.0);
    }

    #[test]
    fn test_cpl() {
        let model = DarkEnergyModel::CPL { w_0: -0.9, w_a: -0.1 };
        let w_at_a = model.w(0.5);
        assert!((w_at_a - (-0.9 - 0.1 * 0.5)).abs() < 1e-10);
    }
}