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

use crate::constants::physical::{C, G, K_B};
use std::f64::consts::PI;

pub fn beta_model_density(rho_0: f64, r: f64, r_c: f64, beta: f64) -> f64 {
    rho_0 * (1.0 + (r / r_c).powi(2)).powf(-3.0 * beta / 2.0)
}

pub fn beta_model_surface_brightness(s_0: f64, theta: f64, theta_c: f64, beta: f64) -> f64 {
    s_0 * (1.0 + (theta / theta_c).powi(2)).powf(-3.0 * beta + 0.5)
}

pub fn beta_model_integrated_mass(rho_0: f64, r_c: f64, beta: f64, r: f64) -> f64 {
    let x = r / r_c;
    4.0 * PI * rho_0 * r_c.powi(3) * hypergeometric_integral(x, beta)
}

fn hypergeometric_integral(x: f64, beta: f64) -> f64 {
    let steps = 500;
    let dx = x / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let xi = (i as f64 + 0.5) * dx;
        sum += xi * xi * (1.0 + xi * xi).powf(-3.0 * beta / 2.0) * dx;
    }
    sum
}

pub fn xray_luminosity_bremsstrahlung(n_e: f64, t_kelvin: f64, volume: f64) -> f64 {
    1.4e-27 * t_kelvin.sqrt() * n_e * n_e * volume
}

pub fn cooling_time(n_e: f64, t_kelvin: f64) -> f64 {
    let lambda = 1.4e-27 * t_kelvin.sqrt();
    1.5 * n_e * K_B * t_kelvin / (n_e * n_e * lambda)
}

pub fn cooling_radius(r_200: f64, t_kelvin: f64, n_e_0: f64, age: f64) -> f64 {
    let t_cool_center = cooling_time(n_e_0, t_kelvin);
    if t_cool_center > age {
        return 0.0;
    }
    r_200 * (t_cool_center / age).sqrt()
}

pub fn sz_y_parameter(n_e: f64, t_kelvin: f64, path_length: f64) -> f64 {
    let sigma_t = 6.6524e-29;
    let me_c2 = 0.511e6 * 1.602e-19;
    sigma_t / me_c2 * n_e * K_B * t_kelvin * path_length
}

pub fn sz_temperature_decrement(y: f64, t_cmb: f64) -> f64 {
    -2.0 * y * t_cmb
}

pub fn sz_kinetic_decrement(tau: f64, v_peculiar: f64, t_cmb: f64) -> f64 {
    -tau * v_peculiar / C * t_cmb
}

pub fn sz_relativistic_correction(y: f64, t_e_kev: f64, t_cmb: f64) -> f64 {
    let x_e = t_e_kev / 511.0;
    -2.0 * y * t_cmb * (1.0 + 3.75 * x_e)
}

pub fn thermal_pressure(n_e: f64, t_kelvin: f64) -> f64 {
    n_e * K_B * t_kelvin
}

pub fn entropy_profile(t_kev: f64, n_e_cm3: f64) -> f64 {
    t_kev / n_e_cm3.powf(2.0 / 3.0)
}

pub const GAS_MASS_FRACTION_UNIVERSAL: f64 = 0.13;

pub const ICM_METALLICITY_TYPICAL: f64 = 0.3;

pub fn iron_mass_in_icm(m_gas: f64, metallicity_solar: f64) -> f64 {
    let iron_solar_mass_fraction = 1.17e-3;
    m_gas * metallicity_solar * iron_solar_mass_fraction
}

pub fn hydrostatic_equilibrium_temperature(m_total: f64, r: f64) -> f64 {
    let mu = 0.6;
    let mp = 1.673e-27;
    G * m_total * mu * mp / (2.0 * K_B * r)
}

pub fn comptonization_parameter_integrated(
    n_e_profile: impl Fn(f64) -> f64,
    t_profile: impl Fn(f64) -> f64,
    r_max: f64,
    steps: usize,
) -> f64 {
    let sigma_t = 6.6524e-29;
    let me_c2 = 0.511e6 * 1.602e-19;
    let dr = r_max / steps as f64;
    let mut y = 0.0;
    for i in 0..steps {
        let r = (i as f64 + 0.5) * dr;
        y += sigma_t * K_B * n_e_profile(r) * t_profile(r) / me_c2 * dr;
    }
    y
}

pub fn xray_surface_brightness_profile(
    n_e_profile: impl Fn(f64) -> f64,
    t_profile: impl Fn(f64) -> f64,
    r_projected: f64,
    r_max: f64,
    steps: usize,
) -> f64 {
    let dl = r_max / steps as f64;
    let mut sb = 0.0;
    for i in 0..steps {
        let l = (i as f64 + 0.5) * dl;
        let r = (r_projected * r_projected + l * l).sqrt();
        if r < r_max {
            let n_e = n_e_profile(r);
            let t = t_profile(r);
            sb += 1.4e-27 * t.sqrt() * n_e * n_e * dl;
        }
    }
    2.0 * sb
}

pub fn cooling_flow_mass_rate(l_x: f64, t_kelvin: f64) -> f64 {
    let mu = 0.6;
    let mp = 1.673e-27;
    2.0 * mu * mp * l_x / (5.0 * K_B * t_kelvin)
}

pub fn bondi_accretion_onto_bcg(m_bcg: f64, rho_icm: f64, cs: f64) -> f64 {
    4.0 * PI * G * G * m_bcg * m_bcg * rho_icm / (cs * cs * cs)
}