use crate::config::parameters::G;
use std::f64::consts::PI;
pub struct PlanetaryInterior {
pub total_mass: f64,
pub total_radius: f64,
pub core_mass_fraction: f64,
pub mantle_mass_fraction: f64,
pub core_density: f64,
pub mantle_density: f64,
}
impl PlanetaryInterior {
pub fn new(
total_mass: f64,
total_radius: f64,
core_frac: f64,
core_density: f64,
mantle_density: f64,
) -> Self {
Self {
total_mass,
total_radius,
core_mass_fraction: core_frac,
mantle_mass_fraction: 1.0 - core_frac,
core_density,
mantle_density,
}
}
pub fn core_radius(&self) -> f64 {
let core_mass = self.total_mass * self.core_mass_fraction;
(3.0 * core_mass / (4.0 * PI * self.core_density)).cbrt()
}
pub fn mantle_thickness(&self) -> f64 {
self.total_radius - self.core_radius()
}
pub fn central_pressure(&self) -> f64 {
let rho_mean = self.mean_density();
3.0 * G * self.total_mass * rho_mean / (4.0 * self.total_radius)
}
pub fn mean_density(&self) -> f64 {
self.total_mass / (4.0 / 3.0 * PI * self.total_radius.powi(3))
}
pub fn moment_of_inertia_factor(&self) -> f64 {
let r_core = self.core_radius();
let r_total = self.total_radius;
let rho_c = self.core_density;
let rho_m = self.mantle_density;
let i_core = 8.0 / 15.0 * PI * rho_c * r_core.powi(5);
let i_mantle = 8.0 / 15.0 * PI * rho_m * (r_total.powi(5) - r_core.powi(5));
(i_core + i_mantle) / (self.total_mass * r_total * r_total)
}
pub fn gravitational_binding_energy(&self) -> f64 {
3.0 * G * self.total_mass * self.total_mass / (5.0 * self.total_radius)
}
pub fn adiabatic_temperature_gradient(&self, surface_temp: f64) -> f64 {
let g = G * self.total_mass / (self.total_radius * self.total_radius);
let cp = 1200.0;
surface_temp + g / cp * self.total_radius * 0.5
}
}
pub fn adams_williamson_density(rho_0: f64, depth: f64, bulk_modulus: f64, g: f64) -> f64 {
rho_0 * (rho_0 * g * depth / bulk_modulus).exp()
}
pub fn core_mass_fraction_from_density(mean_density: f64) -> f64 {
let earth_mean = 5515.0;
let ratio = mean_density / earth_mean;
(0.32 * ratio).clamp(0.0, 0.8)
}
pub fn differentiation_energy(total_mass: f64, radius: f64, density_contrast: f64) -> f64 {
0.1 * density_contrast * G * total_mass * total_mass / radius
}