darkmatter 0.0.1

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

use crate::constants::physical::G;
use std::f64::consts::PI;

pub fn blumenthal_contraction_radius(
    r_initial: f64,
    m_dm_initial: f64,
    m_baryon_final: f64,
    m_dm_conserved: f64,
) -> f64 {
    let m_total_initial = m_dm_initial;
    let m_total_final = m_dm_conserved + m_baryon_final;
    if m_total_final < 1e-100 {
        return r_initial;
    }
    r_initial * m_total_initial / m_total_final
}

pub fn adiabatic_contraction_nfw(
    r_f: f64,
    m_nfw_fn: &dyn Fn(f64) -> f64,
    m_baryon_fn: &dyn Fn(f64) -> f64,
    tol: f64,
    max_iter: usize,
) -> f64 {
    let mut r_i = r_f;
    for _ in 0..max_iter {
        let m_i = m_nfw_fn(r_i);
        let m_b = m_baryon_fn(r_f);
        let r_new = r_f * (m_i + m_b) / m_i;
        if (r_new - r_i).abs() / r_i.max(1e-30) < tol {
            return r_new;
        }
        r_i = r_new;
    }
    r_i
}

pub fn gnedin_contraction_radius(
    r_initial: f64,
    r_bar: f64,
    m_dm_initial: f64,
    m_baryon_final: f64,
    a_gnedin: f64,
) -> f64 {
    let r_eff = a_gnedin * r_initial + (1.0 - a_gnedin) * r_bar;
    let m_total = m_dm_initial + m_baryon_final;
    if m_total < 1e-100 {
        return r_initial;
    }
    r_eff * m_dm_initial / m_total
}

pub fn contracted_density_profile(
    r_f: f64,
    r_f_dr: f64,
    rho_initial_fn: &dyn Fn(f64) -> f64,
    m_nfw_fn: &dyn Fn(f64) -> f64,
    m_baryon_fn: &dyn Fn(f64) -> f64,
    tol: f64,
    max_iter: usize,
) -> f64 {
    let r_i_1 = adiabatic_contraction_nfw(r_f, m_nfw_fn, m_baryon_fn, tol, max_iter);
    let r_i_2 = adiabatic_contraction_nfw(r_f_dr, m_nfw_fn, m_baryon_fn, tol, max_iter);
    let dr_f = r_f_dr - r_f;
    let dr_i = r_i_2 - r_i_1;
    if dr_f.abs() < 1e-50 {
        return rho_initial_fn(r_i_1);
    }
    rho_initial_fn(r_i_1) * (r_i_1 / r_f).powi(2) * (dr_i / dr_f).abs()
}

pub fn baryon_feedback_expansion(r: f64, r_core_kpc: f64, expansion_factor: f64) -> f64 {
    let r_core = r_core_kpc * 3.0857e19;
    r * (1.0 + expansion_factor * (-r * r / (2.0 * r_core * r_core)).exp())
}

pub fn cusp_to_core_transformation(rho_nfw: f64, r: f64, r_core_kpc: f64) -> f64 {
    let r_core = r_core_kpc * 3.0857e19;
    let suppression = r * r / (r * r + r_core * r_core);
    rho_nfw * suppression
}

pub fn dekel_zhao_contracted_profile(
    r: f64,
    rho_s: f64,
    r_s: f64,
    alpha: f64,
    beta: f64,
    gamma: f64,
) -> f64 {
    let x = r / r_s;
    rho_s / (x.powf(gamma) * (1.0 + x.powf(alpha)).powf((beta - gamma) / alpha))
}

pub fn pseudo_adiabatic_contraction(
    r_final: f64,
    m_dm_initial_fn: &dyn Fn(f64) -> f64,
    m_baryon_fn: &dyn Fn(f64) -> f64,
    nu: f64,
) -> f64 {
    let m_b = m_baryon_fn(r_final);
    let m_dm = m_dm_initial_fn(r_final);
    let r_initial = r_final * (m_dm + m_b) / m_dm;
    r_final + nu * (r_initial - r_final)
}

pub fn response_to_baryon_potential(
    r: f64,
    rho_dm_fn: &dyn Fn(f64) -> f64,
    phi_baryon_fn: &dyn Fn(f64) -> f64,
    sigma_v: f64,
) -> f64 {
    let phi_b = phi_baryon_fn(r);
    let rho_0 = rho_dm_fn(r);
    rho_0 * (-phi_b / (sigma_v * sigma_v)).exp()
}

pub fn concentration_mass_relation_contracted(m_200_solar: f64, f_baryon: f64) -> f64 {
    let log_m = (m_200_solar / 1e12).log10();
    let c_base = 10.0_f64.powf(0.905 - 0.101 * log_m);
    c_base * (1.0 + 2.0 * f_baryon)
}

pub fn circular_velocity_contracted(m_enclosed_solar: f64, r_kpc: f64) -> f64 {
    let m = m_enclosed_solar * 1.989e30;
    let r = r_kpc * 3.0857e19;
    (G * m / r).sqrt() / 1e3
}

pub fn nfw_enclosed_mass_contracted(rho_s: f64, r_s: f64, r: f64) -> f64 {
    let x = r / r_s;
    4.0 * PI * rho_s * r_s.powi(3) * ((1.0 + x).ln() - x / (1.0 + x))
}