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

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

pub fn thermal_velocity(t_gev: f64, m_gev: f64) -> f64 {
    (3.0 * t_gev / m_gev).sqrt() * C
}

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

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

pub fn dm_heating_rate(rho_dm: f64, sigma_v: f64, m_dm: f64, f_heat: f64) -> f64 {
    let n = rho_dm / m_dm;
    0.5 * n * n * sigma_v * m_dm * C * C * f_heat
}

pub fn dm_baryon_scattering_heating(
    sigma_dm_b: f64,
    rho_dm: f64,
    rho_b: f64,
    m_dm: f64,
    m_b: f64,
    t_dm: f64,
    t_b: f64,
) -> f64 {
    let n_dm = rho_dm / m_dm;
    let n_b = rho_b / m_b;
    let mu = m_dm * m_b / (m_dm + m_b);
    let v_th = (K_B * (t_dm / m_dm + t_b / m_b)).sqrt();
    n_dm * n_b * sigma_dm_b * v_th * mu * K_B * (t_dm - t_b) / (m_dm + m_b)
}

pub fn self_interaction_transfer_rate(sigma_over_m: f64, rho: f64, v_disp: f64) -> f64 {
    sigma_over_m * rho * v_disp
}

pub fn sidm_core_formation_time(sigma_over_m: f64, rho_core: f64, v_disp: f64) -> f64 {
    1.0 / (sigma_over_m * rho_core * v_disp)
}

pub fn sidm_isothermal_core_density(
    sigma_over_m: f64,
    rho_s: f64,
    r_s: f64,
    v_max: f64,
    age: f64,
) -> f64 {
    let rate = sigma_over_m * rho_s * v_max;
    if rate * age > 1.0 {
        rho_s * r_s / (rate * age * r_s)
    } else {
        rho_s
    }
}

pub fn heat_capacity_dm_gas(n: f64) -> f64 {
    1.5 * n * K_B
}

pub const ADIABATIC_INDEX_DM: f64 = 5.0 / 3.0;

pub fn sound_speed_dm(t_kelvin: f64, m_dm_kg: f64) -> f64 {
    let gamma = ADIABATIC_INDEX_DM;
    (gamma * K_B * t_kelvin / m_dm_kg).sqrt()
}

pub fn phase_space_density(rho: f64, sigma: f64) -> f64 {
    rho / (sigma * sigma * sigma)
}

pub fn tremaine_gunn_bound(m_dm_ev: f64) -> f64 {
    let m_kg = m_dm_ev * 1.783e-36;
    let hbar_cubed = HBAR * HBAR * HBAR;
    m_kg.powi(4) / (2.0 * PI * PI * hbar_cubed)
}