darkmatter 0.0.1

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
139
140
141
142
143
144

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

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct FuzzyDarkMatter {
    pub mass_ev: f64,
}

impl FuzzyDarkMatter {
    pub fn mass_kg(&self) -> f64 {
        self.mass_ev * 1.782_661_92e-36
    }

    pub fn m22(&self) -> f64 {
        self.mass_ev / 1e-22
    }

    pub fn compton_wavelength_m(&self) -> f64 {
        HBAR / (self.mass_kg() * C)
    }

    pub fn de_broglie_wavelength_kpc(&self, v_km_s: f64) -> f64 {
        let v = v_km_s * 1e3;
        HBAR / (self.mass_kg() * v) / 3.0857e19
    }

    pub fn jeans_length_kpc(&self, rho_kg_m3: f64) -> f64 {
        let m = self.mass_kg();
        let lambda_j = (PI * HBAR * HBAR / (m * m * G * rho_kg_m3)).powf(0.25);
        lambda_j / 3.0857e19
    }

    pub fn jeans_mass_solar(&self, rho_kg_m3: f64) -> f64 {
        let l_j = self.jeans_length_kpc(rho_kg_m3) * 3.0857e19;
        4.0 / 3.0 * PI * rho_kg_m3 * (l_j / 2.0).powi(3) / 1.989e30
    }

    pub fn half_mode_mass_solar(&self) -> f64 {
        let m22 = self.m22();
        4.4e7 * m22.powf(-4.0 / 3.0)
    }

    pub fn soliton_core_radius_kpc(&self, m_halo_solar: f64) -> f64 {
        let m22 = self.m22();
        1.6 * m22.powf(-1.0) * (m_halo_solar / 1e9).powf(-1.0 / 3.0)
    }
}

pub fn soliton_density_profile(r_kpc: f64, rho_c: f64, r_c_kpc: f64) -> f64 {
    let x = r_kpc / r_c_kpc;
    rho_c / (1.0 + 0.091 * x * x).powf(8.0)
}

pub fn soliton_central_density(mass_ev: f64, r_c_kpc: f64) -> f64 {
    let m22 = mass_ev / 1e-22;
    let r_c_pc = r_c_kpc * 1e3;
    1.9e7 * m22.powf(-2.0) * (r_c_pc / 100.0).powf(-4.0) * 1.989e30 / 3.0857e19_f64.powi(3)
}

pub fn soliton_mass_solar(rho_c: f64, r_c_kpc: f64) -> f64 {
    let r_c = r_c_kpc * 3.0857e19;
    let volume_factor = 4.0 * PI * r_c.powi(3) * 2.42;
    rho_c * volume_factor / 1.989e30
}

pub fn core_halo_mass_relation_solar(mass_ev: f64, m_halo_solar: f64) -> f64 {
    let m22 = mass_ev / 1e-22;
    1.4e9 * m22.powf(-1.0) * (m_halo_solar / 1e12).powf(1.0 / 3.0)
}

pub fn quantum_pressure(mass_ev: f64, rho_kg_m3: f64, r_m: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    HBAR * HBAR / (2.0 * m * m) * rho_kg_m3 / (r_m * r_m)
}

pub fn gravitational_pressure(rho_kg_m3: f64, r_m: f64) -> f64 {
    G * rho_kg_m3 * rho_kg_m3 * r_m * r_m
}

pub fn quantum_jeans_wavenumber(mass_ev: f64, rho_kg_m3: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    (16.0 * PI * G * rho_kg_m3 * m * m / (HBAR * HBAR)).powf(0.25)
}

pub fn suppression_transfer_function(k_mpc: f64, mass_ev: f64) -> f64 {
    let m22 = mass_ev / 1e-22;
    let k_j = 9.0 * m22.powf(0.5);
    let x = k_mpc / k_j;
    (1.0 + x.powi(2)).powf(-4.0).cos().powi(2) + (-(x / 3.0).powi(2)).exp() * 0.5
}

pub fn interference_fringe_scale_kpc(mass_ev: f64, v_km_s: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    let v = v_km_s * 1e3;
    let lambda_db = HBAR / (m * v);
    lambda_db / 3.0857e19
}

pub fn granule_mass_solar(mass_ev: f64, rho_kg_m3: f64, v_km_s: f64) -> f64 {
    let l = interference_fringe_scale_kpc(mass_ev, v_km_s) * 3.0857e19;
    4.0 / 3.0 * PI * rho_kg_m3 * (l / 2.0).powi(3) / 1.989e30
}

pub fn dynamical_heating_rate(mass_ev: f64, rho_dm_kg_m3: f64, v_km_s: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    let v = v_km_s * 1e3;
    let t_relax = m * v * v / (G * G * rho_dm_kg_m3 * rho_dm_kg_m3) * v / HBAR;
    if t_relax < 1e-50 {
        return f64::INFINITY;
    }
    v * v / t_relax
}

pub fn superradiance_rate(mass_ev: f64, m_bh_solar: f64, spin_bh: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    let m_bh = m_bh_solar * 1.989e30;
    let r_g = G * m_bh / (C * C);
    let alpha = G * m_bh * m / (HBAR * C);
    if alpha > 0.5 {
        return 0.0;
    }
    let omega_h = spin_bh * C / (2.0 * r_g * (1.0 + (1.0 - spin_bh * spin_bh).sqrt()));
    let omega_r = alpha * alpha * m * C * C / (2.0 * HBAR);
    if omega_r >= omega_h {
        return 0.0;
    }
    alpha.powf(7.0) * (omega_h - omega_r)
}

pub fn occupation_number(rho_kg_m3: f64, mass_ev: f64, v_km_s: f64) -> f64 {
    let m = mass_ev * 1.782_661_92e-36;
    let v = v_km_s * 1e3;
    let n_density = rho_kg_m3 / m;
    let p = m * v;
    let phase_vol = (2.0 * PI * HBAR).powi(3);
    n_density * phase_vol / (4.0 * PI * p * p * p / 3.0)
}

pub fn condensation_fraction(t_k: f64, t_c_k: f64) -> f64 {
    if t_k >= t_c_k {
        return 0.0;
    }
    1.0 - (t_k / t_c_k).powf(1.5)
}