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

use crate::constants::cosmic::{HUBBLE_CONSTANT_SI, OMEGA_BARYON, OMEGA_MATTER};
use crate::constants::physical::{C, G, K_B, SIGMA_SB};
use std::f64::consts::PI;

pub fn ionization_fraction_equilibrium(
    n_h_cm3: f64,
    gamma_ion_s: f64,
    alpha_rec_cm3_s: f64,
    t_gas_k: f64,
) -> f64 {
    let alpha = recombination_coefficient_case_b(t_gas_k);
    let alpha_eff = if alpha_rec_cm3_s > 0.0 {
        alpha_rec_cm3_s
    } else {
        alpha
    };
    let discriminant = (gamma_ion_s / (alpha_eff * n_h_cm3)).min(1e20);
    discriminant / (1.0 + discriminant)
}

pub fn recombination_coefficient_case_b(t_gas_k: f64) -> f64 {
    2.6e-13 * (t_gas_k / 1e4).powf(-0.76)
}

pub fn photoionization_rate_dm_annihilation(
    sigma_v_cm3_s: f64,
    rho_dm_gev_cm3: f64,
    m_dm_gev: f64,
    f_ion: f64,
    e_ion_ev: f64,
) -> f64 {
    let n_dm = rho_dm_gev_cm3 / m_dm_gev;
    let energy_rate = 0.5 * n_dm * n_dm * sigma_v_cm3_s * 2.0 * m_dm_gev * 1e9;
    f_ion * energy_rate / (e_ion_ev * 1e6)
}

pub fn photoionization_rate_dm_decay(
    gamma_decay_s: f64,
    rho_dm_gev_cm3: f64,
    m_dm_gev: f64,
    f_ion: f64,
    e_ion_ev: f64,
) -> f64 {
    let n_dm = rho_dm_gev_cm3 / m_dm_gev;
    let energy_rate = n_dm * gamma_decay_s * m_dm_gev * 1e9;
    f_ion * energy_rate / (e_ion_ev * 1e6)
}

pub fn optical_depth_reionization(x_e_fn: &dyn Fn(f64) -> f64, z_max: f64, n_steps: usize) -> f64 {
    let sigma_t = 6.652e-29;
    let n_h0 = 1.9e-7;
    let dz = z_max / n_steps as f64;
    let mut tau = 0.0;
    for i in 0..n_steps {
        let z = (i as f64 + 0.5) * dz;
        let x_e = x_e_fn(z);
        let h_z =
            HUBBLE_CONSTANT_SI * (OMEGA_MATTER * (1.0 + z).powi(3) + (1.0 - OMEGA_MATTER)).sqrt();
        let n_e = x_e * n_h0 * (1.0 + z).powi(3);
        tau += sigma_t * n_e * C / h_z * dz / (1.0 + z);
    }
    tau
}

pub fn clumping_factor(z: f64) -> f64 {
    let c = 26.2917;
    let d = 6.14;
    c * ((1.0 + z) / d).powf(-1.1)
}

pub fn ionizing_photon_rate_dm(
    sigma_v_cm3_s: f64,
    rho_dm_gev_cm3: f64,
    m_dm_gev: f64,
    f_gamma: f64,
) -> f64 {
    let n_dm = rho_dm_gev_cm3 / m_dm_gev;
    0.5 * n_dm * n_dm * sigma_v_cm3_s * f_gamma
}

pub fn stromgren_radius_dm_kpc(q_ion_s: f64, n_h_cm3: f64, alpha_rec: f64) -> f64 {
    let r_cm = (3.0 * q_ion_s / (4.0 * PI * n_h_cm3 * n_h_cm3 * alpha_rec)).powf(1.0 / 3.0);
    r_cm / 3.0857e21
}

pub fn reionization_redshift_approx(n_ion_per_baryon: f64, f_esc: f64, c_hii: f64) -> f64 {
    let f_eff = n_ion_per_baryon * f_esc;
    if f_eff < 1.0 {
        return 6.0;
    }
    6.0 + 4.0 * ((f_eff / c_hii).ln()).max(0.0)
}

pub fn dm_heating_igm_temperature(
    t_igm_k: f64,
    sigma_v_cm3_s: f64,
    rho_dm_gev_cm3: f64,
    m_dm_gev: f64,
    f_heat: f64,
    n_b_cm3: f64,
    dt_s: f64,
) -> f64 {
    let n_dm = rho_dm_gev_cm3 / m_dm_gev;
    let heat_rate_ev_cm3_s = 0.5 * n_dm * n_dm * sigma_v_cm3_s * 2.0 * m_dm_gev * 1e9 * f_heat;
    let heat_per_baryon_ev = heat_rate_ev_cm3_s / n_b_cm3 * dt_s;
    let delta_t = heat_per_baryon_ev * 1.602e-19 / (1.5 * K_B);
    t_igm_k + delta_t
}

pub fn neutral_fraction_post_reionization(z: f64, z_reion: f64) -> f64 {
    if z > z_reion {
        return 1.0;
    }
    1e-4 * ((1.0 + z) / 4.0).powi(4)
}

pub fn lyman_alpha_optical_depth(x_hi: f64, z: f64) -> f64 {
    let tau_gp = 7.16e5 * ((1.0 + z) / 10.0).powf(1.5);
    tau_gp * x_hi
}

pub fn bubble_size_distribution_peak_mpc(z: f64, xi_ion: f64) -> f64 {
    let r_0 = 5.0;
    r_0 * (xi_ion / 40.0).powf(1.0 / 3.0) * ((1.0 + z) / 8.0).powf(-1.0)
}

pub fn mean_baryon_density_kg_m3(z: f64) -> f64 {
    let rho_crit = 3.0 * HUBBLE_CONSTANT_SI * HUBBLE_CONSTANT_SI / (8.0 * PI * G);
    rho_crit * OMEGA_BARYON * (1.0 + z).powi(3)
}

pub fn dm_annihilation_effective_temperature(
    sigma_v_cm3_s: f64,
    rho_dm_gev_cm3: f64,
    m_dm_gev: f64,
    r_sphere_m: f64,
) -> f64 {
    let n_dm = rho_dm_gev_cm3 / m_dm_gev;
    let power_density_gev = 0.5 * n_dm * n_dm * sigma_v_cm3_s * 2.0 * m_dm_gev;
    let luminosity_w = power_density_gev * 1.602e-10 * 1e6 * 4.0 / 3.0 * PI * r_sphere_m.powi(3);
    let area = 4.0 * PI * r_sphere_m * r_sphere_m;
    (luminosity_w / (SIGMA_SB * area)).powf(0.25)
}