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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

use crate::constants::cosmic::{CMB_TEMPERATURE_K, HUBBLE_CONSTANT_SI, OMEGA_BARYON, OMEGA_DM};
use crate::constants::physical::{C, HBAR, K_B};
use std::f64::consts::{E, PI};

pub fn isocurvature_power_spectrum_cdm(k_mpc: f64, amplitude: f64, n_iso: f64) -> f64 {
    let k_pivot = 0.05;
    amplitude * (k_mpc / k_pivot).powf(n_iso - 1.0)
}

pub fn isocurvature_power_spectrum_baryon(k_mpc: f64, amplitude: f64, n_iso: f64) -> f64 {
    isocurvature_power_spectrum_cdm(k_mpc, amplitude, n_iso)
}

pub fn isocurvature_fraction(alpha_iso: f64, a_s: f64) -> f64 {
    alpha_iso / (1.0 + alpha_iso) * a_s
}

pub fn silk_damping_scale_mpc(omega_b_h2: f64, omega_m_h2: f64) -> f64 {
    8.130 * (omega_m_h2 / 0.15).powf(-0.5) * (omega_b_h2 / 0.023).powf(-0.5)
}

pub fn silk_damping_mass_solar(omega_b_h2: f64, omega_m_h2: f64) -> f64 {
    let rd_mpc = silk_damping_scale_mpc(omega_b_h2, omega_m_h2);
    let rd_m = rd_mpc * 3.0857e22;
    let rho_crit =
        3.0 * HUBBLE_CONSTANT_SI * HUBBLE_CONSTANT_SI / (8.0 * PI * crate::constants::physical::G);
    let rho_m = rho_crit * (omega_m_h2 / (HUBBLE_CONSTANT_SI * 3.0857e22 / (1e5)).powi(2));
    4.0 / 3.0 * PI * (rd_m / 2.0).powi(3) * rho_m / 1.989e30
}

pub fn silk_damping_envelope(k_mpc: f64, k_silk_mpc: f64) -> f64 {
    (-(k_mpc / k_silk_mpc).powi(2)).exp()
}

pub fn bao_correlation_function(r_mpc: f64, r_bao_mpc: f64, sigma_bao_mpc: f64) -> f64 {
    let dr = r_mpc - r_bao_mpc;
    0.05 * (-dr * dr / (2.0 * sigma_bao_mpc * sigma_bao_mpc)).exp()
        / (2.0 * PI * sigma_bao_mpc * sigma_bao_mpc).sqrt()
        + xi_smooth(r_mpc)
}

fn xi_smooth(r_mpc: f64) -> f64 {
    if r_mpc < 1.0 {
        return 1.0;
    }
    (1.0 / r_mpc).powf(1.8)
}

pub fn bao_peak_position(omega_m_h2: f64, omega_b_h2: f64) -> f64 {
    let omega_r_h2 = 4.15e-5 * (CMB_TEMPERATURE_K / 2.7255).powi(4);
    let z_eq = omega_m_h2 / omega_r_h2 - 1.0;
    let z_d = 1059.68;
    let rb = 3.0 * omega_b_h2 / (4.0 * omega_r_h2 * z_d);
    let c_s_avg = C / (3.0 * (1.0 + rb)).sqrt();
    let radiation_correction = 1.0 / (1.0 + (1.0 + z_d) / z_eq).sqrt();
    c_s_avg * radiation_correction / HUBBLE_CONSTANT_SI / (1.0 + z_d).sqrt()
        * (2.0 / (3.0 * omega_m_h2 * 100.0 * 1e3 / 3.0857e22).sqrt())
        / 3.0857e22
}

pub fn transfer_function_eisenstein_hu(
    k_mpc: f64,
    omega_m_h2: f64,
    omega_b_h2: f64,
    h: f64,
) -> f64 {
    let theta = CMB_TEMPERATURE_K / 2.7;
    let z_eq = 2.5e4 * omega_m_h2 * theta.powf(-4.0);
    let k_eq = 0.0746 * omega_m_h2 * theta.powf(-2.0);

    let b1 = 0.313 * omega_m_h2.powf(-0.419) * (1.0 + 0.607 * omega_m_h2.powf(0.674));
    let b2 = 0.238 * omega_m_h2.powf(0.223);
    let z_d = 1291.0 * omega_m_h2.powf(0.251) / (1.0 + 0.659 * omega_m_h2.powf(0.828))
        * (1.0 + b1 * omega_b_h2.powf(b2));

    let r_d = 31.5e3 * omega_b_h2 * theta.powf(-4.0) / z_d;
    let r_eq = 31.5e3 * omega_b_h2 * theta.powf(-4.0) / z_eq;
    let s =
        2.0 / (3.0 * k_eq) * (6.0 / r_eq).sqrt() * ((1.0 + r_d).sqrt() + (r_d + r_eq).sqrt()).ln()
            / (1.0 + (r_eq).sqrt()).ln();

    let q = k_mpc * theta * theta / (omega_m_h2 * h);
    let f_b = omega_b_h2 / omega_m_h2;
    let f_c = 1.0 - f_b;

    let alpha_c = f_c.powf(-0.567) * (1.0 + (f_b / f_c).powf(0.5));
    let l = (E + 1.8 * alpha_c * q).ln();
    let c_val = 14.2 / alpha_c + 386.0 / (1.0 + 69.9 * q.powf(1.08));
    let t0 = l / (l + c_val * q * q);

    let ks = k_mpc * s;
    let j0 = ks.sin() / ks;
    f_c * t0 + f_b * j0 * (-0.5 * (k_mpc / 5.0).powi(2)).exp()
}

pub fn matter_power_spectrum(
    k_mpc: f64,
    a_s: f64,
    n_s: f64,
    omega_m_h2: f64,
    omega_b_h2: f64,
    h: f64,
) -> f64 {
    let k_pivot = 0.05;
    let p_prim = a_s * (k_mpc / k_pivot).powf(n_s - 1.0);
    let t = transfer_function_eisenstein_hu(k_mpc, omega_m_h2, omega_b_h2, h);
    p_prim * t * t * k_mpc
}

pub fn growth_factor_approx(z: f64, omega_m: f64, omega_lambda: f64) -> f64 {
    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 * 5.0 * omega_m_z
        / (2.0
            * (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 sigma_8_from_power(
    a_s: f64,
    n_s: f64,
    omega_m_h2: f64,
    omega_b_h2: f64,
    h: f64,
    n_steps: usize,
) -> f64 {
    let r = 8.0;
    let k_min = 1e-4_f64;
    let k_max = 10.0_f64;
    let dlnk = (k_max / k_min).ln() / n_steps as f64;
    let mut integral = 0.0;
    for i in 0..n_steps {
        let lnk = k_min.ln() + (i as f64 + 0.5) * dlnk;
        let k = lnk.exp();
        let pk = matter_power_spectrum(k, a_s, n_s, omega_m_h2, omega_b_h2, h);
        let x = k * r;
        let w = 3.0 * (x.sin() - x * x.cos()) / (x * x * x);
        integral += pk * w * w * k * k * k * dlnk;
    }
    (integral / (2.0 * PI * PI)).sqrt()
}

pub fn baryon_to_dm_ratio() -> f64 {
    OMEGA_BARYON / OMEGA_DM
}

pub fn de_broglie_wavelength_dm_mpc(m_dm_gev: f64, v_km_s: f64) -> f64 {
    let m_kg = m_dm_gev * 1.782_661_92e-27;
    let v = v_km_s * 1e3;
    HBAR / (m_kg * v) / 3.0857e22
}

pub fn thermal_velocity_dm_km_s(m_dm_gev: f64, t_kev: f64) -> f64 {
    let m_kg = m_dm_gev * 1.782_661_92e-27;
    let t_k = t_kev * 1e3 * 1.602e-19 / K_B;
    (3.0 * K_B * t_k / m_kg).sqrt() / 1e3
}