planetsfactory 0.0.3

Planet factory — classify, build and catalogue planets for any star system: Solar System, TRAPPIST-1, Kepler-90, Proxima Centauri, or fully custom worlds.
Documentation
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
}