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

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

pub fn dynamical_friction_timescale(
    m_satellite: f64,
    m_cluster: f64,
    r_orbit: f64,
    v_orbit: f64,
    ln_lambda: f64,
) -> f64 {
    let rho = 3.0 * m_cluster / (4.0 * PI * r_orbit.powi(3));
    let coulomb = 4.0 * PI * G * G * m_satellite * rho * ln_lambda / (v_orbit * v_orbit * v_orbit);
    m_satellite * v_orbit / coulomb
}

pub fn tidal_stripping_radius(m_satellite: f64, m_cluster: f64, d_cluster_center: f64) -> f64 {
    d_cluster_center * (m_satellite / (3.0 * m_cluster)).powf(1.0 / 3.0)
}

pub fn ram_pressure_stripping_condition(
    rho_icm: f64,
    v_galaxy: f64,
    sigma_star: f64,
    sigma_gas: f64,
) -> f64 {
    let p_ram = rho_icm * v_galaxy * v_galaxy;
    let p_restoring = 2.0 * PI * G * sigma_star * sigma_gas;
    p_ram / p_restoring
}

pub fn galaxy_harassment_energy(
    m_perturber: f64,
    b_impact: f64,
    v_encounter: f64,
    r_galaxy: f64,
) -> f64 {
    let delta_v = 2.0 * G * m_perturber / (b_impact * v_encounter);
    0.5 * delta_v * delta_v * r_galaxy
}

pub fn merger_timescale_lacey_cole(
    m1: f64,
    m2: f64,
    r_vir: f64,
    v_vir: f64,
    ln_lambda: f64,
) -> f64 {
    let mass_ratio = m2 / m1;
    let f_orbit = 0.5;
    let orbital_energy = r_vir * v_vir * v_vir / (G * (m1 + m2));
    f_orbit * 1.17 / mass_ratio / ln_lambda * orbital_energy * r_vir / v_vir
}

pub fn bullet_cluster_velocity_estimate(m1: f64, m2: f64, r_sep: f64) -> f64 {
    (2.0 * G * (m1 + m2) / r_sep).sqrt()
}

pub fn bullet_cluster_constraints_sidm(v_collision: f64, separation_dm_gas: f64) -> f64 {
    separation_dm_gas / (v_collision * v_collision) * C * C
}

pub fn dm_self_interaction_offset(sigma_over_m: f64, sigma_dm: f64, v_collision: f64) -> f64 {
    let tau = sigma_over_m * sigma_dm;
    tau * v_collision
}

pub const SUBHALO_MASS_FUNCTION_SLOPE: f64 = -0.9;

pub fn subhalo_abundance_in_cluster(m_cluster_solar: f64, m_sub_min_solar: f64) -> f64 {
    let alpha = SUBHALO_MASS_FUNCTION_SLOPE;
    let m_ratio = m_sub_min_solar / m_cluster_solar;
    0.01 * m_ratio.powf(alpha + 1.0) * m_cluster_solar / m_sub_min_solar
}

pub fn velocity_dispersion_from_mass(m_200: f64, r_200: f64) -> f64 {
    (G * m_200 / (3.0 * r_200)).sqrt()
}

pub fn crossing_time(r_vir: f64, sigma_v: f64) -> f64 {
    2.0 * r_vir / sigma_v
}

pub fn relaxation_criterion(offset: f64, r_200: f64) -> bool {
    offset / r_200 < 0.07
}

pub fn bautz_morgan_type(bcg_dominance: f64) -> &'static str {
    if bcg_dominance > 2.0 {
        "I"
    } else if bcg_dominance > 1.5 {
        "I-II"
    } else if bcg_dominance > 1.0 {
        "II"
    } else {
        "III"
    }
}

pub fn infall_velocity_at_turnaround(m_cluster: f64, r_turnaround: f64) -> f64 {
    (2.0 * G * m_cluster / r_turnaround).sqrt()
}

pub fn virial_shock_temperature(v_infall: f64) -> f64 {
    let mu = 0.6;
    let mp = 1.673e-27;
    let k_b = 1.381e-23;
    3.0 * mu * mp * v_infall * v_infall / (16.0 * k_b)
}

pub fn cluster_pair_binding_energy(m1: f64, m2: f64, separation: f64) -> f64 {
    G * m1 * m2 / separation
}

pub fn supercluster_mass_estimate(n_clusters: f64, avg_mass: f64) -> f64 {
    n_clusters * avg_mass * SOLAR_MASS
}

pub fn filament_bridge_dm_density(m1: f64, m2: f64, separation: f64, d_from_center: f64) -> f64 {
    let m_eff = (m1 * m2).sqrt();
    let rho = G * m_eff / (C * C * separation * separation);
    rho * (-d_from_center * d_from_center / (0.1 * separation * separation)).exp()
}