darkmatter 0.0.2

Dark matter simulation engine — gravitational fields, particle dynamics, halo stability, and cosmological constants
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

use crate::constants::cosmic::{HUBBLE_CONSTANT, OMEGA_BARYON, OMEGA_DM, OMEGA_LAMBDA};
use crate::constants::physical::G;
use std::f64::consts::PI;

pub fn hubble_parameter(z: f64) -> f64 {
    let omega_m = OMEGA_DM + OMEGA_BARYON;
    let omega_r = 9.1e-5;
    let h_sq = HUBBLE_CONSTANT
        * HUBBLE_CONSTANT
        * (omega_r * (1.0 + z).powi(4) + omega_m * (1.0 + z).powi(3) + OMEGA_LAMBDA);
    h_sq.sqrt()
}

pub fn linear_growth_factor(z: f64) -> f64 {
    let omega_m = OMEGA_DM + OMEGA_BARYON;
    let a = 1.0 / (1.0 + z);
    let omega_m_z = omega_m * (1.0 + z).powi(3) / (omega_m * (1.0 + z).powi(3) + OMEGA_LAMBDA);
    let omega_l_z = OMEGA_LAMBDA / (omega_m * (1.0 + z).powi(3) + OMEGA_LAMBDA);
    a * 2.5 * omega_m_z
        / (omega_m_z.powf(4.0 / 7.0) - omega_l_z
            + (1.0 + omega_m_z / 2.0) * (1.0 + omega_l_z / 70.0))
}

pub fn growth_rate_f(z: f64) -> f64 {
    let omega_m = OMEGA_DM + OMEGA_BARYON;
    let omega_m_z = omega_m * (1.0 + z).powi(3) / (omega_m * (1.0 + z).powi(3) + OMEGA_LAMBDA);
    omega_m_z.powf(0.55)
}

pub fn matter_power_spectrum_primordial(k: f64, a_s: f64, n_s: f64, k_pivot: f64) -> f64 {
    a_s * (k / k_pivot).powf(n_s - 1.0)
}

pub fn transfer_function_bbks(k_mpc: f64, omega_m_h2: f64) -> f64 {
    let q = k_mpc / (omega_m_h2 * (2.7255_f64 / 2.7).powi(2));
    let t = (2.34 * q).ln() / (2.34 * q)
        * (1.0 + 3.89 * q + (16.1 * q).powi(2) + (5.46 * q).powi(3) + (6.71 * q).powi(4))
            .powf(-0.25);
    if q < 1e-10 { 1.0 } else { t }
}

pub fn sigma_r_squared(
    r: f64,
    p_k: impl Fn(f64) -> f64,
    k_min: f64,
    k_max: f64,
    steps: usize,
) -> f64 {
    let dk = (k_max / k_min).ln() / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let ln_k = (k_min).ln() + (i as f64 + 0.5) * dk;
        let k = ln_k.exp();
        let x = k * r;
        let w = if x < 0.01 {
            1.0 - x * x / 10.0
        } else {
            3.0 * (x.sin() - x * x.cos()) / (x * x * x)
        };
        sum += k * k * k * p_k(k) * w * w * dk / (2.0 * PI * PI);
    }
    sum
}

pub fn press_schechter_mass_function(
    m: f64,
    rho_m: f64,
    sigma: f64,
    d_sigma_d_m: f64,
    delta_c: f64,
) -> f64 {
    let nu = delta_c / sigma;
    let f_nu = (2.0 / PI).sqrt() * nu * (-nu * nu / 2.0).exp();
    rho_m / (m * m) * f_nu * (d_sigma_d_m / sigma).abs() * m
}

pub fn sheth_tormen_mass_function(
    m: f64,
    rho_m: f64,
    sigma: f64,
    d_sigma_d_m: f64,
    delta_c: f64,
) -> f64 {
    let a_st: f64 = 0.707;
    let p_st = 0.3;
    let aa = 0.3222;
    let nu = delta_c / sigma;
    let nu_prime = a_st.sqrt() * nu;
    let f_nu = aa
        * (1.0 + 1.0 / nu_prime.powf(2.0 * p_st))
        * (2.0 / PI).sqrt()
        * nu_prime
        * (-nu_prime * nu_prime / 2.0).exp();
    rho_m / (m * m) * f_nu * (d_sigma_d_m / sigma).abs() * m
}

pub fn collapse_redshift(delta: f64, z_initial: f64) -> f64 {
    let d_i = linear_growth_factor(z_initial);
    let d_c = 1.686;
    let d_coll = d_i * delta / d_c;
    let mut z = 0.0;
    let mut step = 50.0;
    while step > 0.001 {
        if linear_growth_factor(z) > d_coll {
            z += step;
        }
        step *= 0.5;
        if linear_growth_factor(z) < d_coll {
            z -= step;
        }
    }
    z.max(0.0)
}

pub fn jeans_length(cs: f64, rho: f64) -> f64 {
    cs * (PI / (G * rho)).sqrt()
}

pub fn jeans_mass(cs: f64, rho: f64) -> f64 {
    let lj = jeans_length(cs, rho);
    4.0 / 3.0 * PI * rho * (lj / 2.0).powi(3)
}

pub fn free_streaming_length(m_dm_ev: f64, _t_decoupling_gev: f64) -> f64 {
    let mpc = 3.0857e22;
    0.2 * (30.0 / (m_dm_ev * 1e-9)).sqrt() * mpc
}

pub fn warm_dm_cutoff_mass(m_dm_kev: f64) -> f64 {
    1e10 * (1.0 / m_dm_kev).powf(3.33) * 1.989e30
}