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

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

pub fn stellar_to_halo_mass_behroozi(m_halo_solar: f64) -> f64 {
    let m1 = 1e12;
    let epsilon = 0.02;
    let alpha = -1.412;
    let delta_val = 0.608;
    let gamma_val = 1.124;
    let x = (m_halo_solar / m1).log10();
    let f_x = -x.abs().powf(alpha) + delta_val * (x.exp() / (1.0 + x.exp())).powf(gamma_val);
    m1 * epsilon * 10.0_f64.powf(f_x)
}

pub fn abundance_matching_mstar(m_halo_solar: f64) -> f64 {
    let m_pivot = 1e12;
    let m_star_pivot = 5e10;
    let alpha_low = 2.0;
    let alpha_high = 0.5;
    if m_halo_solar < m_pivot {
        m_star_pivot * (m_halo_solar / m_pivot).powf(alpha_low)
    } else {
        m_star_pivot * (m_halo_solar / m_pivot).powf(alpha_high)
    }
}

pub fn halo_spin_parameter(j: f64, e: f64, m: f64) -> f64 {
    j * e.abs().sqrt() / (G * m.powf(2.5))
}

pub fn disk_scale_radius(lambda: f64, r_vir: f64, f_r: f64) -> f64 {
    lambda / (2.0_f64).sqrt() * r_vir * f_r
}

pub fn disk_scale_height(r_d: f64, ratio: f64) -> f64 {
    r_d * ratio
}

pub fn tully_fisher(v_flat_km_s: f64) -> f64 {
    let log_l = 4.0 * v_flat_km_s.log10() - 1.5;
    10.0_f64.powf(log_l)
}

pub fn faber_jackson(sigma_km_s: f64) -> f64 {
    let log_l = 4.0 * sigma_km_s.log10() - 0.5;
    10.0_f64.powf(log_l)
}

pub fn adiabatic_contraction_blumenthal(r_i: f64, m_dm_i: f64, m_baryon: f64, r_f: f64) -> f64 {
    r_i * m_dm_i / (m_dm_i + m_baryon * (1.0 - r_f / r_i).max(0.0))
}

pub fn nfw_to_burkert_from_sidm(
    rho_s: f64,
    r_s: f64,
    sigma_over_m: f64,
    age_gyr: f64,
) -> (f64, f64) {
    let age = age_gyr * 3.156e16;
    let v_max = (4.625 * G * rho_s * r_s * r_s).sqrt();
    let rate = sigma_over_m * rho_s * v_max;
    let core_ratio = (rate * age).min(3.0);
    let r_core = r_s * core_ratio;
    let rho_core = rho_s * r_s / (r_core + r_s);
    (rho_core, r_core)
}

pub fn dwarf_spheroidal_mass_within_half_light(sigma_v: f64, r_half: f64) -> f64 {
    4.0 * sigma_v * sigma_v * r_half / G
}

pub fn wolf_mass_estimator(sigma_los: f64, r_half: f64) -> f64 {
    930.0 * sigma_los * sigma_los * r_half
}

pub fn dm_fraction_within_r(m_total: f64, m_stars: f64, m_gas: f64) -> f64 {
    (m_total - m_stars - m_gas).max(0.0) / m_total
}

pub fn too_big_to_fail_threshold_vmax(v_max_km_s: f64) -> bool {
    v_max_km_s > 30.0
}

pub fn missing_satellites_predicted(m_host_solar: f64, m_sub_min_solar: f64) -> f64 {
    0.008 * (m_host_solar / m_sub_min_solar).powf(0.9)
}

pub fn virial_mass_from_circular_velocity(v_c: f64, r_vir: f64) -> f64 {
    v_c * v_c * r_vir / G
}

pub fn angular_momentum_from_spin(lambda: f64, m_vir: f64, r_vir: f64) -> f64 {
    let v_vir = (G * m_vir / r_vir).sqrt();
    lambda * m_vir * v_vir * r_vir
}

pub fn specific_angular_momentum(j_total: f64, m: f64) -> f64 {
    j_total / m
}

pub fn maximum_disk_hypothesis(v_max: f64, r_d: f64) -> f64 {
    0.85 * v_max * v_max * r_d / G
}

pub fn schwarzschild_radius_galactic(m: f64) -> f64 {
    2.0 * G * m / (C * C)
}

pub fn galactic_escape_velocity(m_vir: f64, r: f64) -> f64 {
    (2.0 * G * m_vir / r).sqrt()
}

pub fn dm_annihilation_luminosity_galaxy(
    rho_dm_profile: impl Fn(f64) -> f64,
    sigma_v: f64,
    m_dm: f64,
    r_max: f64,
    steps: usize,
) -> f64 {
    let dr = r_max / steps as f64;
    let mut luminosity = 0.0;
    for i in 0..steps {
        let r = (i as f64 + 0.5) * dr;
        let rho = rho_dm_profile(r);
        let n = rho / m_dm;
        let shell_volume = 4.0 * PI * r * r * dr;
        luminosity += 0.5 * n * n * sigma_v * m_dm * C * C * shell_volume;
    }
    luminosity
}