use sciforge::hub::domain::astronomy::orbits::escape_velocity;
use sciforge::hub::domain::common::constants::{EARTH_MASS, EARTH_RADIUS, G};
pub struct Impactor {
pub mass_kg: f64,
pub velocity_m_s: f64,
pub density_kg_m3: f64,
pub angle_deg: f64,
}
impl Impactor {
pub fn asteroid(diameter_m: f64, velocity_km_s: f64) -> Self {
let density = crate::STONY_ASTEROID_DENSITY;
let radius = diameter_m / 2.0;
let volume = (4.0 / 3.0) * std::f64::consts::PI * radius.powi(3);
Self {
mass_kg: density * volume,
velocity_m_s: velocity_km_s * 1000.0,
density_kg_m3: density,
angle_deg: 45.0,
}
}
pub fn kinetic_energy_j(&self) -> f64 {
0.5 * self.mass_kg * self.velocity_m_s.powi(2)
}
pub fn kinetic_energy_mt(&self) -> f64 {
self.kinetic_energy_j() / crate::MT_TNT_TO_JOULE
}
pub fn impact_velocity(&self) -> f64 {
let v_esc = escape_velocity(G * EARTH_MASS, EARTH_RADIUS);
(self.velocity_m_s.powi(2) + v_esc.powi(2)).sqrt()
}
pub fn crater_diameter_m(&self, target_density: f64) -> f64 {
let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
let angle_rad = self.angle_deg.to_radians();
let d_projectile =
(6.0 * self.mass_kg / (std::f64::consts::PI * self.density_kg_m3)).cbrt();
1.161
* (self.density_kg_m3 / target_density).powf(1.0 / 3.0)
* d_projectile.powf(0.78)
* self.velocity_m_s.powf(0.44)
* g.powf(-0.22)
* angle_rad.sin().powf(1.0 / 3.0)
}
pub fn fireball_radius_m(&self) -> f64 {
let energy_kt = self.kinetic_energy_j() / crate::KT_TNT_TO_JOULE;
55.0 * energy_kt.powf(0.4)
}
pub fn ejecta_volume_m3(&self, target_density: f64) -> f64 {
let crater_d = self.crater_diameter_m(target_density);
let crater_depth = crater_d / 5.0;
std::f64::consts::PI / 6.0 * crater_d * crater_d * crater_depth
}
}
pub fn chicxulub_equivalent() -> Impactor {
Impactor::asteroid(10_000.0, 20.0)
}
pub fn tunguska_equivalent() -> Impactor {
Impactor::asteroid(60.0, 15.0)
}