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

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

pub fn annihilation_luminosity(rho: f64, sigma_v: f64, m_dm: f64, volume: f64) -> f64 {
    let n = rho / m_dm;
    0.5 * n * n * sigma_v * m_dm * C * C * volume
}

pub fn j_factor_nfw(rho_s: f64, r_s: f64, distance: f64, angle_rad: f64) -> f64 {
    let r_max = distance * angle_rad;
    let steps = 200;
    let dr = r_max / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let r = (i as f64 + 0.5) * dr;
        let rho = crate::halos::profiles::nfw(r, rho_s, r_s);
        sum += rho * rho * dr;
    }
    2.0 * PI * angle_rad * angle_rad * sum
}

pub fn d_factor_nfw(rho_s: f64, r_s: f64, distance: f64, angle_rad: f64) -> f64 {
    let r_max = distance * angle_rad;
    let steps = 200;
    let dr = r_max / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let r = (i as f64 + 0.5) * dr;
        let rho = crate::halos::profiles::nfw(r, rho_s, r_s);
        sum += rho * dr;
    }
    2.0 * PI * angle_rad * angle_rad * sum
}

pub fn gamma_ray_flux(j_factor: f64, sigma_v: f64, m_dm: f64, n_gamma: f64, distance: f64) -> f64 {
    sigma_v * n_gamma * j_factor / (8.0 * PI * m_dm * m_dm * distance * distance)
}

pub fn decay_flux(d_factor: f64, decay_rate: f64, m_dm: f64, n_gamma: f64, distance: f64) -> f64 {
    decay_rate * n_gamma * d_factor / (4.0 * PI * m_dm * distance * distance)
}

pub fn positron_injection_rate(rho: f64, sigma_v: f64, m_dm: f64, branching: f64) -> f64 {
    let n = rho / m_dm;
    0.5 * n * n * sigma_v * branching
}

pub fn neutrino_flux_from_sun(
    capture_rate: f64,
    annihilation_fraction: f64,
    n_nu_per_ann: f64,
    distance_sun: f64,
) -> f64 {
    capture_rate * annihilation_fraction * n_nu_per_ann / (4.0 * PI * distance_sun * distance_sun)
}

pub fn sommerfeld_enhancement(alpha_coupling: f64, v: f64) -> f64 {
    let epsilon = alpha_coupling / v * C;
    if epsilon < 0.01 {
        return 1.0;
    }
    let pi_eps = PI * epsilon;
    pi_eps / (1.0 - (-pi_eps).exp())
}

pub const THERMAL_RELIC_SIGMA_V: f64 = 3e-26;

pub fn boost_factor_substructure(concentration: f64) -> f64 {
    let c = concentration;
    let fc = (1.0 + c).ln() - c / (1.0 + c);
    c * c * c / (9.0 * fc * fc) * (1.0 - 1.0 / (1.0 + c).powi(3))
}

pub fn inverse_compton_power(n_dm: f64, sigma_v: f64, e_electron: f64, u_photon: f64) -> f64 {
    let sigma_t = 6.6524e-29;
    let rate = 0.5 * n_dm * n_dm * sigma_v;
    rate * (4.0 / 3.0) * sigma_t * C * (e_electron / (0.511e6 * 1.602e-19)).powi(2) * u_photon
}