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

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

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

impl PrimordialBlackHole {
    pub fn from_solar_masses(m_solar: f64) -> Self {
        Self {
            mass_kg: m_solar * SOLAR_MASS,
        }
    }

    pub fn schwarzschild_radius(&self) -> f64 {
        2.0 * G * self.mass_kg / (C * C)
    }

    pub fn hawking_temperature_gev(&self) -> f64 {
        let hbar_si = 1.0546e-34;
        let kb = 1.381e-23;
        let t_kelvin = hbar_si * C * C * C / (8.0 * PI * G * self.mass_kg * kb);
        t_kelvin * 8.617e-14 * 1e-9
    }

    pub fn hawking_luminosity(&self) -> f64 {
        let hbar_si = 1.0546e-34;
        hbar_si * C * C * C * C * C * C / (15360.0 * PI * G * G * self.mass_kg * self.mass_kg)
    }

    pub fn evaporation_time_seconds(&self) -> f64 {
        let hbar_si = 1.0546e-34;
        5120.0 * PI * G * G * self.mass_kg.powi(3) / (hbar_si * C * C * C * C)
    }

    pub fn is_stable_now(&self) -> bool {
        self.evaporation_time_seconds() > 4.35e17
    }

    pub fn mass_solar(&self) -> f64 {
        self.mass_kg / SOLAR_MASS
    }

    pub fn formation_mass_from_horizon(t_seconds: f64) -> f64 {
        C * C * C * t_seconds / (2.0 * G)
    }

    pub fn microlensing_einstein_radius(&self, d_l: f64, d_s: f64) -> f64 {
        let d_ls = d_s - d_l;
        if d_ls <= 0.0 || d_s <= 0.0 {
            return 0.0;
        }
        (4.0 * G * self.mass_kg * d_ls / (C * C * d_l * d_s)).sqrt()
    }

    pub fn microlensing_timescale(&self, d_l: f64, d_s: f64, v_transverse: f64) -> f64 {
        let r_e = self.microlensing_einstein_radius(d_l, d_s);
        if v_transverse <= 0.0 {
            return f64::INFINITY;
        }
        r_e / v_transverse
    }

    pub fn accretion_luminosity(&self, accretion_rate_kg_s: f64, efficiency: f64) -> f64 {
        efficiency * accretion_rate_kg_s * C * C
    }

    pub fn bondi_accretion_rate(&self, rho_gas: f64, cs_gas: f64) -> f64 {
        4.0 * PI * G * G * self.mass_kg * self.mass_kg * rho_gas / (cs_gas * cs_gas * cs_gas)
    }
}

pub fn pbh_dm_fraction_constraint_evaporation(mass_kg: f64) -> f64 {
    let pbh = PrimordialBlackHole { mass_kg };
    let t_evap = pbh.evaporation_time_seconds();
    let age_universe = 4.35e17;
    if t_evap < age_universe { 0.0 } else { 1.0 }
}

pub fn pbh_mass_range_dm_solar() -> (f64, f64) {
    (1e-16, 1e3)
}

pub fn pbh_abundance_beta(mass_kg: f64, f_dm: f64) -> f64 {
    let m_eq = 2.8e17 * SOLAR_MASS;
    f_dm * (mass_kg / m_eq).sqrt()
}

pub fn extended_mass_function_lognormal(mass: f64, mass_c: f64, sigma: f64) -> f64 {
    let log_ratio = (mass / mass_c).ln();
    (-log_ratio * log_ratio / (2.0 * sigma * sigma)).exp() / (mass * sigma * (2.0 * PI).sqrt())
}

pub fn gravitational_wave_from_merger(m1_kg: f64, m2_kg: f64, distance: f64) -> f64 {
    let m_total = m1_kg + m2_kg;
    let mu = m1_kg * m2_kg / m_total;
    let mc = mu.powf(3.0 / 5.0) * m_total.powf(2.0 / 5.0);
    4.0 * G * mc / (C * C * distance) * (G * mc / (C * C * C)).powf(5.0 / 3.0)
}

pub fn pbh_merger_rate_per_gpc3_yr(f_dm: f64, mass_solar: f64) -> f64 {
    let base_rate = 1.3e6;
    base_rate * f_dm.powf(53.0 / 37.0) * mass_solar.powf(-32.0 / 37.0)
}