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::G;
use std::f64::consts::PI;

pub fn nfw(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    rho_s / (x * (1.0 + x) * (1.0 + x))
}

pub fn nfw_enclosed_mass(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    4.0 * PI * rho_s * r_s * r_s * r_s * ((1.0 + x).ln() - x / (1.0 + x))
}

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

pub fn nfw_concentration_duffy(mass_solar: f64, redshift: f64) -> f64 {
    let pivot = 2e12;
    5.71 * (mass_solar / pivot).powf(-0.084) * (1.0 + redshift).powf(-0.47)
}

pub fn hernquist(r: f64, total_mass: f64, a: f64) -> f64 {
    total_mass * a / (2.0 * PI * r * (r + a).powi(3))
}

pub fn hernquist_enclosed_mass(r: f64, total_mass: f64, a: f64) -> f64 {
    total_mass * r * r / ((r + a) * (r + a))
}

pub fn hernquist_potential(r: f64, total_mass: f64, a: f64) -> f64 {
    -G * total_mass / (r + a)
}

pub fn burkert(r: f64, rho_0: f64, r_0: f64) -> f64 {
    rho_0 / ((1.0 + r / r_0) * (1.0 + (r / r_0) * (r / r_0)))
}

pub fn burkert_enclosed_mass(r: f64, rho_0: f64, r_0: f64) -> f64 {
    let x = r / r_0;
    PI * rho_0 * r_0 * r_0 * r_0 * ((1.0 + x * x).ln() + 2.0 * (1.0 + x).ln() - 2.0 * x.atan())
}

pub fn isothermal(sigma_v: f64, r: f64) -> f64 {
    sigma_v * sigma_v / (2.0 * PI * G * r * r)
}

pub fn isothermal_enclosed_mass(sigma_v: f64, r: f64) -> f64 {
    2.0 * sigma_v * sigma_v * r / G
}

pub fn einasto(r: f64, rho_e: f64, r_e: f64, alpha: f64) -> f64 {
    let d_n = 3.0 * alpha - 1.0 / 3.0 + 0.0079 / alpha;
    rho_e * (-d_n * ((r / r_e).powf(1.0 / alpha) - 1.0)).exp()
}

pub fn moore(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    rho_s / (x.powf(1.5) * (1.0 + x).powf(1.5))
}

pub fn pseudo_isothermal(r: f64, rho_0: f64, r_c: f64) -> f64 {
    rho_0 / (1.0 + (r / r_c) * (r / r_c))
}

pub fn pseudo_isothermal_enclosed_mass(r: f64, rho_0: f64, r_c: f64) -> f64 {
    4.0 * PI * rho_0 * r_c * r_c * (r - r_c * (r / r_c).atan())
}

pub fn plummer(r: f64, total_mass: f64, b: f64) -> f64 {
    3.0 * total_mass / (4.0 * PI * b * b * b) * (1.0 + (r / b) * (r / b)).powf(-2.5)
}

pub fn dekel_zhao(r: f64, rho_s: f64, r_s: f64, alpha: f64, beta: f64, gamma: f64) -> f64 {
    let x = r / r_s;
    rho_s / (x.powf(gamma) * (1.0 + x.powf(alpha)).powf((beta - gamma) / alpha))
}

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct HaloProfile {
    pub virial_mass: f64,
    pub virial_radius: f64,
    pub concentration: f64,
    pub scale_radius: f64,
    pub rho_s: f64,
}

impl HaloProfile {
    pub fn from_nfw(virial_mass: f64, concentration: f64, redshift: f64) -> Self {
        let h_si = crate::constants::cosmic::HUBBLE_CONSTANT_SI;
        let omega_m = crate::constants::cosmic::OMEGA_MATTER;
        let omega_l = crate::constants::cosmic::OMEGA_LAMBDA;
        let hz = h_si * (omega_m * (1.0 + redshift).powi(3) + omega_l).sqrt();
        let rho_crit = 3.0 * hz * hz / (8.0 * PI * G);
        let delta_vir = 200.0;
        let r_vir = (3.0 * virial_mass / (4.0 * PI * delta_vir * rho_crit)).powf(1.0 / 3.0);
        let r_s = r_vir / concentration;
        let c = concentration;
        let rho_s = delta_vir * rho_crit * c * c * c / (3.0 * ((1.0 + c).ln() - c / (1.0 + c)));
        Self {
            virial_mass,
            virial_radius: r_vir,
            concentration,
            scale_radius: r_s,
            rho_s,
        }
    }

    pub fn density_at(&self, r: f64) -> f64 {
        nfw(r, self.rho_s, self.scale_radius)
    }

    pub fn enclosed_mass_at(&self, r: f64) -> f64 {
        nfw_enclosed_mass(r, self.rho_s, self.scale_radius)
    }

    pub fn circular_velocity_at(&self, r: f64) -> f64 {
        (G * self.enclosed_mass_at(r) / r).sqrt()
    }

    pub fn potential_at(&self, r: f64) -> f64 {
        nfw_potential(r, self.rho_s, self.scale_radius)
    }

    pub fn v_max(&self) -> f64 {
        let r_max = 2.163 * self.scale_radius;
        self.circular_velocity_at(r_max)
    }
}