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

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

pub const RHO_LOCAL_GEV_CM3: f64 = 0.3;
pub const RHO_LOCAL_KG_M3: f64 = 0.3 * 1.783e-27 * 1e6;

pub const V0_KM_S: f64 = 220.0;
pub const V_SUN_KM_S: f64 = 232.0;
pub const V_ESC_KM_S: f64 = 544.0;
pub const V_EARTH_ORBITAL_KM_S: f64 = 29.8;

pub const V0_M_S: f64 = V0_KM_S * 1e3;
pub const V_ESC_M_S: f64 = V_ESC_KM_S * 1e3;

pub fn maxwell_boltzmann_speed(v: f64, v0: f64) -> f64 {
    let norm = 4.0 * PI / (PI.sqrt() * v0 * v0 * v0);
    norm * v * v * (-v * v / (v0 * v0)).exp()
}

pub fn mean_speed(v0: f64) -> f64 {
    2.0 * v0 / PI.sqrt()
}

pub fn rms_speed(v0: f64) -> f64 {
    v0 * (3.0 / 2.0_f64).sqrt()
}

pub fn velocity_distribution_galactic_frame(v: f64, v0: f64, v_esc: f64) -> f64 {
    if v > v_esc {
        return 0.0;
    }
    let norm_factor = crate::utils::math::erf_approx(v_esc / v0)
        - 2.0 * v_esc / (PI.sqrt() * v0) * (-v_esc * v_esc / (v0 * v0)).exp();
    if norm_factor <= 0.0 {
        return 0.0;
    }
    let f = (-v * v / (v0 * v0)).exp() / (PI.sqrt() * v0 * v0 * v0 * norm_factor);
    4.0 * PI * v * v * f
}

pub fn velocity_in_lab_frame(v_galactic: f64, v_sun: f64, v_earth: f64, phase_days: f64) -> f64 {
    let omega = 2.0 * PI / 365.25;
    let cos_angle = (omega * phase_days).cos();
    let v_lab_sq = v_galactic * v_galactic
        + v_sun * v_sun
        + v_earth * v_earth
        + 2.0 * v_galactic * v_sun
        + 2.0 * v_earth * (v_galactic + v_sun) * 0.49 * cos_angle;
    v_lab_sq.abs().sqrt()
}

pub const MODULATION_PEAK_DAY: f64 = 152.0;

pub fn annual_modulation_fraction(v_earth: f64, v_sun: f64, v0: f64) -> f64 {
    v_earth * v_sun / (v0 * v0)
}

pub fn dm_flux_per_cm2_s(rho_gev_cm3: f64, m_dm_gev: f64, v_km_s: f64) -> f64 {
    let n = rho_gev_cm3 / m_dm_gev;
    n * v_km_s * 1e5
}

pub fn dm_wind_direction_galactic() -> (f64, f64, f64) {
    (0.0, -V_SUN_KM_S * 1e3, 0.0)
}

pub fn local_dm_column_density(path_length_kpc: f64) -> f64 {
    RHO_LOCAL_KG_M3 * path_length_kpc * KPC
}

pub fn phase_space_density_local() -> f64 {
    let rho = RHO_LOCAL_GEV_CM3;
    let sigma = V0_KM_S * 1e5;
    rho / (sigma * sigma * sigma)
}

pub fn number_density_local(m_dm_gev: f64) -> f64 {
    RHO_LOCAL_GEV_CM3 / m_dm_gev * 1e6
}