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

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

pub fn tidal_tensor_diagonal(m_host: f64, r: f64) -> f64 {
    -G * m_host / (r * r * r)
}

pub fn tidal_tensor_radial(m_host: f64, r: f64) -> f64 {
    2.0 * G * m_host / (r * r * r)
}

pub fn tidal_radius(m_sub: f64, m_host: f64, d: f64) -> f64 {
    d * (m_sub / (3.0 * m_host)).powf(1.0 / 3.0)
}

pub fn roche_lobe_radius(m_sub: f64, m_host: f64, d: f64) -> f64 {
    let q = m_sub / m_host;
    let q_13 = q.powf(1.0 / 3.0);
    d * 0.462 * q_13 / (0.6 * q_13 * q_13 + (1.0 + q_13).ln())
}

pub fn tidal_heating_energy(
    m_sub: f64,
    r_half: f64,
    m_host: f64,
    r_peri: f64,
    v_peri_km_s: f64,
) -> f64 {
    let v_peri = v_peri_km_s * 1e3;
    let omega_tidal = (G * m_host / (r_peri * r_peri * r_peri)).sqrt();
    let dt = r_peri / v_peri;
    let i_sub = 0.4 * m_sub * r_half * r_half;
    i_sub * omega_tidal * omega_tidal * dt * dt
}

pub fn mass_loss_rate_king(m_sub: f64, r_tidal: f64, r_half: f64, t_orbit_gyr: f64) -> f64 {
    let t_orbit = t_orbit_gyr * 3.156e16;
    if r_tidal > 10.0 * r_half {
        return 0.0;
    }
    let fraction_lost = ((r_half / r_tidal).powi(3)).min(1.0);
    m_sub * fraction_lost / t_orbit
}

pub fn tidal_stripping_mass_remaining(
    m_initial: f64,
    rho_host_at_r: f64,
    rho_sub_at_half: f64,
    n_orbits: f64,
) -> f64 {
    if rho_host_at_r < rho_sub_at_half {
        return m_initial;
    }
    let eta = (rho_host_at_r / rho_sub_at_half).sqrt();
    m_initial * (-eta * n_orbits * 0.1).exp()
}

pub fn tidal_shock_energy(
    m_sub: f64,
    r_half: f64,
    m_perturber: f64,
    impact_param: f64,
    v_rel_km_s: f64,
) -> f64 {
    let v_rel = v_rel_km_s * 1e3;
    let b = impact_param;
    2.0 * G * G * m_perturber * m_perturber * m_sub * r_half * r_half
        / (3.0 * v_rel * v_rel * b.powi(4))
}

pub fn disk_shocking_heating(m_sub: f64, r_half: f64, sigma_disk_kg_m2: f64, v_z_km_s: f64) -> f64 {
    let v_z = v_z_km_s * 1e3;
    let g_z = 2.0 * PI * G * sigma_disk_kg_m2;
    let delta_v = g_z * r_half / v_z;
    0.5 * m_sub * delta_v * delta_v / (m_sub + 1e-100) * m_sub
}

pub fn survival_criterion(
    m_sub: f64,
    r_half: f64,
    m_host: f64,
    r_peri: f64,
    v_peri_km_s: f64,
) -> bool {
    let e_tidal = tidal_heating_energy(m_sub, r_half, m_host, r_peri, v_peri_km_s);
    let e_bind = 0.4 * G * m_sub * m_sub / r_half;
    e_tidal < 0.3 * e_bind
}

pub fn tidal_stream_width_kpc(m_sub: f64, m_host: f64, r_orbit_kpc: f64) -> f64 {
    let r_t = tidal_radius(m_sub, m_host, r_orbit_kpc * 3.0857e19);
    r_t / 3.0857e19
}

pub fn tidal_stream_velocity_dispersion_km_s(m_sub: f64, r_tidal_kpc: f64) -> f64 {
    let r_t = r_tidal_kpc * 3.0857e19;
    (G * m_sub / r_t).sqrt() / 1e3
}

pub fn orbital_decay_from_tidal(
    m_sub: f64,
    m_host: f64,
    r_orbit: f64,
    e_tidal_per_orbit: f64,
) -> f64 {
    let e_orbit = -G * m_host * m_sub / (2.0 * r_orbit);
    if e_orbit.abs() < 1e-100 {
        return 0.0;
    }
    -r_orbit * e_tidal_per_orbit / (2.0 * e_orbit.abs())
}