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
131
132
133
134
135
136
137
138

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

pub fn two_point_correlation(r_mpc: f64, r_0: f64, gamma: f64) -> f64 {
    (r_mpc / r_0).powf(-gamma)
}

pub const CORRELATION_LENGTH_DM: f64 = 5.0;

pub fn galaxy_bias(xi_galaxy: f64, xi_dm: f64) -> f64 {
    if xi_dm > 0.0 {
        (xi_galaxy / xi_dm).sqrt()
    } else {
        1.0
    }
}

pub fn bias_from_peak_height(nu: f64) -> f64 {
    let delta_c = 1.686;
    1.0 + (nu * nu - 1.0) / delta_c
}

pub fn power_spectrum_to_correlation(
    p_k: impl Fn(f64) -> f64,
    r: f64,
    k_min: f64,
    k_max: f64,
    steps: usize,
) -> f64 {
    let dk = (k_max / k_min).ln() / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let ln_k = k_min.ln() + (i as f64 + 0.5) * dk;
        let k = ln_k.exp();
        let kr = k * r;
        let sinc = if kr.abs() < 1e-10 { 1.0 } else { kr.sin() / kr };
        sum += k * k * k * p_k(k) * sinc * dk / (2.0 * PI * PI);
    }
    sum
}

pub fn sigma_8_from_power_spectrum(
    p_k: impl Fn(f64) -> f64,
    k_min: f64,
    k_max: f64,
    steps: usize,
) -> f64 {
    let r = 8.0;
    let dk = (k_max / k_min).ln() / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let ln_k = k_min.ln() + (i as f64 + 0.5) * dk;
        let k = ln_k.exp();
        let x = k * r;
        let w = if x < 0.01 {
            1.0 - x * x / 10.0
        } else {
            3.0 * (x.sin() - x * x.cos()) / (x * x * x)
        };
        sum += k * k * k * p_k(k) * w * w * dk / (2.0 * PI * PI);
    }
    sum.sqrt()
}

pub const BAO_SCALE_MPC: f64 = 147.0;

pub fn angular_diameter_distance(z: f64, h0_km_s_mpc: f64, omega_m: f64) -> f64 {
    let c_km_s = C * 1e-3;
    let mpc = 3.0857e22;
    let steps = 1000;
    let dz = z / steps as f64;
    let mut integral = 0.0;
    for i in 0..steps {
        let zi = (i as f64 + 0.5) * dz;
        let ez = (omega_m * (1.0 + zi).powi(3) + (1.0 - omega_m)).sqrt();
        integral += dz / ez;
    }
    c_km_s / h0_km_s_mpc * integral / (1.0 + z) * mpc * 1e3
}

pub fn luminosity_distance(z: f64, h0_km_s_mpc: f64, omega_m: f64) -> f64 {
    angular_diameter_distance(z, h0_km_s_mpc, omega_m) * (1.0 + z) * (1.0 + z)
}

pub fn comoving_distance(z: f64, h0_km_s_mpc: f64, omega_m: f64) -> f64 {
    angular_diameter_distance(z, h0_km_s_mpc, omega_m) * (1.0 + z)
}

pub fn three_point_correlation_hierarchical(xi_12: f64, xi_23: f64, xi_31: f64, q: f64) -> f64 {
    q * (xi_12 * xi_23 + xi_23 * xi_31 + xi_31 * xi_12)
}

pub fn halo_model_one_halo_power(
    k: f64,
    m_min: f64,
    m_max: f64,
    dn_dm: impl Fn(f64) -> f64,
    u_k_m: impl Fn(f64, f64) -> f64,
    rho_mean: f64,
    steps: usize,
) -> f64 {
    let dlnm = ((m_max / m_min).ln()) / steps as f64;
    let mut sum = 0.0;
    for i in 0..steps {
        let lnm = m_min.ln() + (i as f64 + 0.5) * dlnm;
        let m = lnm.exp();
        let uk = u_k_m(k, m);
        sum += dn_dm(m) * m * m * uk * uk * dlnm / (rho_mean * rho_mean);
    }
    sum
}

pub fn halo_model_two_halo_power(k: f64, p_lin: f64, bias1: f64, bias2: f64) -> f64 {
    let scale_dependent_correction = 1.0 / (1.0 + (k * 0.1).powi(2));
    bias1 * bias2 * p_lin * scale_dependent_correction
}

pub fn nfw_fourier_profile(k: f64, m: f64, rho_s: f64, r_s: f64, c: f64) -> f64 {
    let r_vir = c * r_s;
    let m_nfw = 4.0 * PI * rho_s * r_s.powi(3) * ((1.0 + c).ln() - c / (1.0 + c));
    let norm = m / m_nfw;
    let steps = 200;
    let dr = r_vir / steps as f64;
    let mut integral = 0.0;
    for i in 0..steps {
        let r = (i as f64 + 0.5) * dr;
        let x = r / r_s;
        let rho = rho_s / (x * (1.0 + x) * (1.0 + x));
        let kr = k * r;
        let sinc = if kr < 1e-10 { 1.0 } else { kr.sin() / kr };
        integral += 4.0 * PI * r * r * rho * sinc * dr;
    }
    norm * integral / m
}

pub fn gravitational_redshift_cluster(m: f64, r: f64) -> f64 {
    G * m / (r * C * C)
}