use crate::constants::physical::{C, G, SOLAR_MASS};
use std::f64::consts::PI;
pub fn capture_rate_sun(rho_local_kg_m3: f64, v_dm: f64, sigma_si_cm2: f64, m_dm_kg: f64) -> f64 {
let r_sun = 6.957e8;
let m_sun = SOLAR_MASS;
let m_h = 1.673e-27;
let n_h = m_sun / m_h;
let sigma_si = sigma_si_cm2 * 1e-4;
let n_dm = rho_local_kg_m3 / m_dm_kg;
let v_esc_sun = (2.0 * G * m_sun / r_sun).sqrt();
let geometric = PI * r_sun * r_sun;
let grav_focus = 1.0 + (v_esc_sun / v_dm).powi(2);
let prob = (n_h * sigma_si / (PI * r_sun * r_sun)).min(1.0);
n_dm * v_dm * geometric * grav_focus * prob
}
pub fn capture_rate_planet(
rho_local_kg_m3: f64,
v_dm: f64,
sigma_si_cm2: f64,
m_dm_kg: f64,
planet_mass: f64,
planet_radius: f64,
) -> f64 {
let m_nucleon = 1.673e-27;
let n_nucleons = planet_mass / m_nucleon;
let sigma_si = sigma_si_cm2 * 1e-4;
let n_dm = rho_local_kg_m3 / m_dm_kg;
let v_esc = (2.0 * G * planet_mass / planet_radius).sqrt();
let geometric = PI * planet_radius * planet_radius;
let grav_focus = 1.0 + (v_esc / v_dm).powi(2);
let prob = (n_nucleons * sigma_si / (PI * planet_radius * planet_radius)).min(1.0);
n_dm * v_dm * geometric * grav_focus * prob
}
pub fn equilibrium_annihilation_rate(capture_rate: f64, sigma_v: f64, v_eff: f64) -> f64 {
let ca = sigma_v / v_eff.powi(3);
let tau = 1.0 / (capture_rate * ca).sqrt();
0.5 * capture_rate * (capture_rate * ca * tau * tau / 1.0).tanh()
}
pub fn dm_accumulated_number(
capture_rate: f64,
age_seconds: f64,
sigma_v: f64,
volume: f64,
) -> f64 {
let ann_rate = sigma_v / volume;
let tau_eq = 1.0 / (capture_rate * ann_rate).sqrt();
let t_ratio = age_seconds / tau_eq;
(capture_rate / ann_rate).sqrt() * t_ratio.tanh()
}
pub fn dm_density_inside_body(
n_captured: f64,
body_radius: f64,
thermal_radius_fraction: f64,
) -> f64 {
let r_th = body_radius * thermal_radius_fraction;
3.0 * n_captured / (4.0 * PI * r_th * r_th * r_th)
}
pub fn heating_from_capture(capture_rate: f64, m_dm_kg: f64) -> f64 {
capture_rate * m_dm_kg * C * C
}
pub fn neutrino_signal_from_capture(capture_rate: f64, branching_nu: f64, n_nu: f64) -> f64 {
capture_rate * branching_nu * n_nu
}
pub fn gravitational_focusing_factor(v_dm: f64, v_esc: f64) -> f64 {
1.0 + (v_esc / v_dm).powi(2)
}