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

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

pub fn void_density_profile_hsw(delta_c: f64, r: f64, r_v: f64) -> f64 {
    delta_c * (1.0 - (r / r_v).powi(2))
}

pub fn void_effective_radius(volume: f64) -> f64 {
    (3.0 * volume / (4.0 * PI)).powf(1.0 / 3.0)
}

pub const VOID_UNDERDENSITY_MEAN: f64 = -0.8;

pub fn void_expansion_velocity(r: f64, r_v: f64, hubble: f64, delta_v: f64) -> f64 {
    hubble * r * (1.0 - delta_v / 3.0) * (r / r_v)
}

pub fn void_number_density(r_min_mpc: f64) -> f64 {
    let r0 = 15.0;
    1e-6 * (r_min_mpc / r0).powf(-3.5)
}

pub fn void_size_function_svdw(
    delta_v: f64,
    delta_c: f64,
    sigma_v: f64,
    d_sigma_dr: f64,
    r_v: f64,
    rho_mean: f64,
) -> f64 {
    let nu_v = delta_v.abs() / sigma_v;
    let d = delta_c / delta_v.abs();
    let volume = 4.0 / 3.0 * PI * r_v.powi(3);
    let f_ln_sigma =
        (2.0 / PI).sqrt() * nu_v * (-nu_v * nu_v / 2.0).exp() * (1.0 + d.powi(2)).powf(-1.0);
    rho_mean / (volume * volume) * f_ln_sigma * (d_sigma_dr / sigma_v).abs() * r_v
}

pub fn alcock_paczynski_void(r_parallel: f64, r_perp: f64) -> f64 {
    r_parallel / r_perp
}

pub fn void_galaxy_ccf_amplitude(delta_mean: f64) -> f64 {
    -delta_mean
}

pub fn isw_signal_from_void(r_v_mpc: f64, delta_v: f64, hubble: f64) -> f64 {
    let mpc = 3.0857e22;
    let r_v = r_v_mpc * mpc;
    -2.0 * delta_v * hubble * r_v / 3.0
}

pub fn lensing_convergence_void(delta_v: f64, r_v: f64, sigma_cr: f64, rho_mean: f64) -> f64 {
    rho_mean * delta_v * r_v / sigma_cr
}

pub fn void_dm_density(rho_mean: f64, delta_v: f64) -> f64 {
    rho_mean * (1.0 + delta_v)
}

pub fn void_volume_fraction(delta_threshold: f64) -> f64 {
    if delta_threshold < -0.5 { 0.6 } else { 0.4 }
}

pub fn fifth_force_enhancement_in_void(coupling: f64) -> f64 {
    1.0 + 2.0 * coupling * coupling
}

pub fn void_mass_enclosed(rho_mean: f64, delta_v: f64, r_v: f64) -> f64 {
    4.0 / 3.0 * PI * r_v.powi(3) * rho_mean * (1.0 + delta_v)
}

pub fn void_gravitational_potential(rho_mean: f64, delta_v: f64, r_v: f64) -> f64 {
    let m = void_mass_enclosed(rho_mean, delta_v, r_v);
    -G * m / r_v
}

pub fn void_tidal_field(rho_mean: f64, delta_v: f64, r: f64, r_v: f64) -> f64 {
    let m = void_mass_enclosed(rho_mean, delta_v, r_v);
    2.0 * G * m / (r * r * r)
}

pub fn void_profile_ridged(rho_mean: f64, delta_c: f64, r: f64, r_v: f64, r_wall: f64) -> f64 {
    if r < r_v {
        rho_mean * (1.0 + delta_c * (1.0 - (r / r_v).powi(2)))
    } else if r < r_wall {
        let overdensity = -delta_c * r_v.powi(3) / (r_wall.powi(3) - r_v.powi(3));
        rho_mean * (1.0 + overdensity)
    } else {
        rho_mean
    }
}