darkmatter 0.0.2

Dark matter simulation engine — gravitational fields, particle dynamics, halo stability, and cosmological constants
Documentation
pub fn virial_ratio(kinetic_energy: f64, potential_energy: f64) -> f64 {
    2.0 * kinetic_energy / potential_energy.abs()
}

pub fn is_virialized(kinetic_energy: f64, potential_energy: f64) -> bool {
    let ratio = virial_ratio(kinetic_energy, potential_energy);
    (ratio - 1.0).abs() < 0.1
}

pub fn virial_temperature(total_mass: f64, radius: f64) -> f64 {
    let g = crate::constants::physical::G;
    let k_b = crate::constants::physical::K_B;
    let mu = 0.6;
    let m_p = 1.672_621_923_7e-27;
    mu * m_p * g * total_mass / (2.0 * k_b * radius)
}

pub fn crossing_time(radius: f64, velocity_dispersion: f64) -> f64 {
    2.0 * radius / velocity_dispersion
}

pub fn relaxation_time(n_particles: f64, crossing: f64) -> f64 {
    n_particles / (8.0 * n_particles.ln()) * crossing
}

pub fn dynamical_friction_timescale(
    m_satellite: f64,
    m_host: f64,
    radius: f64,
    v_circ: f64,
    coulomb_log: f64,
) -> f64 {
    let g = crate::constants::physical::G;
    1.17 * radius * radius * v_circ / (g * m_satellite * coulomb_log) * (m_host / m_satellite)
}

pub fn tidal_stripping_criterion(rho_sub: f64, rho_host_at_r: f64) -> bool {
    rho_sub > rho_host_at_r
}

pub fn concentration_from_mass(mass_solar: f64) -> f64 {
    let log_m = (mass_solar / 1e12).log10();
    10.0_f64.powf(0.905 - 0.101 * log_m)
}

pub fn halo_spin_parameter(angular_momentum: f64, energy: f64, mass: f64) -> f64 {
    let g = crate::constants::physical::G;
    angular_momentum * energy.abs().sqrt() / (g * mass.powf(2.5))
}

pub fn nfw_scale_radius(virial_radius: f64, concentration: f64) -> f64 {
    virial_radius / concentration
}

pub fn virial_radius(mass: f64, redshift: f64) -> f64 {
    let g = crate::constants::physical::G;
    let h_si = crate::constants::cosmic::HUBBLE_CONSTANT_SI;
    let hz = h_si
        * (crate::constants::cosmic::OMEGA_MATTER * (1.0 + redshift).powi(3)
            + crate::constants::cosmic::OMEGA_LAMBDA)
            .sqrt();
    let rho_crit = 3.0 * hz * hz / (8.0 * std::f64::consts::PI * g);
    let delta_vir = 200.0;
    (3.0 * mass / (4.0 * std::f64::consts::PI * delta_vir * rho_crit)).powf(1.0 / 3.0)
}