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

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

pub fn circular_velocity_nfw(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    let m_enc = 4.0 * PI * rho_s * r_s.powi(3) * ((1.0 + x).ln() - x / (1.0 + x));
    (G * m_enc / r).sqrt()
}

pub fn circular_velocity_disk(r: f64, m_disk: f64, r_d: f64) -> f64 {
    let y = r / (2.0 * r_d);
    let v_sq = 4.0 * PI * G * m_disk / r_d
        * y
        * y
        * (crate::utils::math::bessel_i0(y) * crate::utils::math::bessel_k0(y)
            - crate::utils::math::bessel_i1(y) * crate::utils::math::bessel_k1(y));
    v_sq.abs().sqrt()
}

pub fn circular_velocity_bulge(r: f64, m_bulge: f64, a: f64) -> f64 {
    (G * m_bulge * r / ((r + a) * (r + a))).sqrt()
}

pub fn circular_velocity_gas(r: f64, m_gas: f64, r_gas: f64) -> f64 {
    let y = r / (2.0 * r_gas);
    let v_sq = 4.0 * PI * G * m_gas / r_gas
        * y
        * y
        * (crate::utils::math::bessel_i0(y) * crate::utils::math::bessel_k0(y)
            - crate::utils::math::bessel_i1(y) * crate::utils::math::bessel_k1(y));
    v_sq.abs().sqrt()
}

pub fn total_rotation_curve(
    r: f64,
    rho_s: f64,
    r_s: f64,
    m_disk: f64,
    r_d: f64,
    m_bulge: f64,
    a_bulge: f64,
) -> f64 {
    let v_dm = circular_velocity_nfw(r, rho_s, r_s);
    let v_disk = circular_velocity_disk(r, m_disk, r_d);
    let v_bulge = circular_velocity_bulge(r, m_bulge, a_bulge);
    (v_dm * v_dm + v_disk * v_disk + v_bulge * v_bulge).sqrt()
}

pub fn flat_rotation_velocity(v_outer: f64) -> f64 {
    v_outer
}

pub fn dm_fraction_from_rotation(v_total: f64, v_baryon: f64) -> f64 {
    let v_dm_sq = v_total * v_total - v_baryon * v_baryon;
    if v_dm_sq > 0.0 {
        v_dm_sq / (v_total * v_total)
    } else {
        0.0
    }
}

pub fn mond_interpolation_simple(a: f64, a_0: f64) -> f64 {
    let x = a / a_0;
    x / (1.0 + x * x).sqrt()
}

pub const MOND_A0: f64 = 1.2e-10;

pub fn radial_acceleration_relation(g_obs: f64, g_bar: f64) -> f64 {
    let a0 = MOND_A0;
    let predicted = g_bar / (1.0 - (-((g_bar / a0).sqrt())).exp());
    let residual = g_obs - predicted;
    predicted + residual * 0.0
}

pub fn mond_effective_mass(g_obs: f64, r: f64) -> f64 {
    g_obs * r * r / G
}

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

pub fn dm_halo_mass_from_flat_v(v_flat: f64, r_vir: f64) -> f64 {
    v_flat * v_flat * r_vir / G
}

pub fn enclosed_mass_from_rotation(v_c: f64, r: f64) -> f64 {
    v_c * v_c * r / G
}

pub fn asymptotic_velocity_nfw(rho_s: f64, r_s: f64) -> f64 {
    (4.0 * PI * G * rho_s * r_s * r_s).sqrt()
}

pub fn v_max_nfw(rho_s: f64, r_s: f64) -> f64 {
    let factor = (2.163_f64).ln() - 2.163 / 3.163;
    (4.0 * PI * G * rho_s * r_s * r_s * factor / 2.163).sqrt()
}

pub fn relativistic_rotation_correction(v: f64) -> f64 {
    let beta = v / C;
    v / (1.0 - beta * beta).sqrt()
}