use sciforge::hub::domain::geology::tectonics::{
    airy_root, heat_flow, isostatic_equilibrium, thermal_subsidence,
};

pub struct MercuryInterior {
    pub core_radius_m: f64,
    pub inner_core_radius_m: f64,
    pub mantle_thickness_m: f64,
    pub crust_thickness_m: f64,
}

impl Default for MercuryInterior {
    fn default() -> Self {
        Self::new()
    }
}

impl MercuryInterior {
    pub fn new() -> Self {
        Self {
            core_radius_m: crate::CORE_RADIUS,
            inner_core_radius_m: crate::INNER_CORE_RADIUS,
            mantle_thickness_m: crate::MANTLE_THICKNESS,
            crust_thickness_m: crate::CRUST_THICKNESS,
        }
    }

    pub fn core_volume_m3(&self) -> f64 {
        (4.0 / 3.0) * std::f64::consts::PI * self.core_radius_m.powi(3)
    }

    pub fn inner_core_volume_m3(&self) -> f64 {
        (4.0 / 3.0) * std::f64::consts::PI * self.inner_core_radius_m.powi(3)
    }

    pub fn mantle_volume_m3(&self) -> f64 {
        let r_outer = self.core_radius_m + self.mantle_thickness_m;
        let r_inner = self.core_radius_m;
        (4.0 / 3.0) * std::f64::consts::PI * (r_outer.powi(3) - r_inner.powi(3))
    }

    pub fn crust_volume_m3(&self) -> f64 {
        let r_planet = crate::physics::orbit::MERCURY_RADIUS;
        let r_inner = r_planet - self.crust_thickness_m;
        (4.0 / 3.0) * std::f64::consts::PI * (r_planet.powi(3) - r_inner.powi(3))
    }

    pub fn core_mass_fraction(&self) -> f64 {
        let core_vol = self.core_volume_m3();
        let total_vol =
            (4.0 / 3.0) * std::f64::consts::PI * crate::physics::orbit::MERCURY_RADIUS.powi(3);
        (core_vol * crate::CORE_DENSITY) / (total_vol * crate::MEAN_DENSITY)
    }

    pub fn is_tectonically_active(&self) -> bool {
        false
    }

    pub fn contraction_estimate_km(&self) -> f64 {
        7.0
    }

    pub fn crustal_root_depth(&self, topography_m: f64) -> f64 {
        airy_root(topography_m, crate::CRUST_DENSITY, crate::MANTLE_DENSITY)
    }

    pub fn isostatic_elevation(&self) -> f64 {
        isostatic_equilibrium(
            self.crust_thickness_m,
            crate::CRUST_DENSITY,
            crate::MANTLE_DENSITY,
        )
    }
}

pub fn surface_heat_flow(conductivity_w_mk: f64, gradient_k_per_m: f64) -> f64 {
    heat_flow(conductivity_w_mk, gradient_k_per_m)
}

pub fn temperature_at_depth(depth_m: f64, gradient_k_per_m: f64) -> f64 {
    crate::MEAN_SURFACE_TEMP_K + depth_m * gradient_k_per_m
}

pub fn cmb_temperature() -> f64 {
    crate::CMB_TEMPERATURE
}

pub fn lithosphere_thermal_subsidence(
    initial_thickness_m: f64,
    time_s: f64,
    thermal_diffusivity: f64,
) -> f64 {
    thermal_subsidence(initial_thickness_m, time_s, thermal_diffusivity)
}

pub fn mean_heat_flux_mw_m2() -> f64 {
    crate::HEAT_FLUX_MW_M2
}
