suns 0.0.1

use sciforge::hub::domain::common::constants::{G, K_B, SIGMA_SB, SOLAR_MASS};

use crate::{
    PROTON_MASS, SOLAR_CORE_DENSITY, SOLAR_CORE_TEMP, SOLAR_EFFECTIVE_TEMP,
    SOLAR_LUMINOSITY, SOLAR_MEAN_DENSITY, SOLAR_RADIUS,
};

pub const CORE_HYDROGEN_FRACTION: f64 = 0.34;
pub const ENVELOPE_HYDROGEN_FRACTION: f64 = 0.7346;
pub const ENVELOPE_HELIUM_FRACTION: f64 = 0.2485;
pub const ENVELOPE_METAL_FRACTION: f64 = 0.0169;
pub const MEAN_MOLECULAR_WEIGHT_CORE: f64 = 0.83;
pub const MEAN_MOLECULAR_WEIGHT_ENVELOPE: f64 = 0.61;
pub const ADIABATIC_INDEX: f64 = 5.0 / 3.0;
pub const EDDINGTON_LUMINOSITY_FRACTION: f64 = 5.1e-5;

pub struct StellarLayer {
    pub name: &'static str,
    pub inner_radius_frac: f64,
    pub outer_radius_frac: f64,
    pub temp_inner_k: f64,
    pub temp_outer_k: f64,
    pub density_inner: f64,
    pub density_outer: f64,
}

impl StellarLayer {
    pub fn mean_temperature(&self) -> f64 {
        (self.temp_inner_k + self.temp_outer_k) / 2.0
    }

    pub fn mean_density(&self) -> f64 {
        (self.density_inner + self.density_outer) / 2.0
    }

    pub fn thickness_m(&self) -> f64 {
        (self.outer_radius_frac - self.inner_radius_frac) * SOLAR_RADIUS
    }

    pub fn volume_m3(&self) -> f64 {
        let r_outer = self.outer_radius_frac * SOLAR_RADIUS;
        let r_inner = self.inner_radius_frac * SOLAR_RADIUS;
        (4.0 / 3.0) * std::f64::consts::PI * (r_outer.powi(3) - r_inner.powi(3))
    }

    pub fn mass_estimate_kg(&self) -> f64 {
        self.mean_density() * self.volume_m3()
    }

    pub fn pressure_at_base(&self) -> f64 {
        let rho = self.mean_density();
        let r = self.inner_radius_frac * SOLAR_RADIUS;
        let m_enclosed = (4.0 / 3.0) * std::f64::consts::PI * r.powi(3) * rho;
        G * m_enclosed * rho / r
    }

    pub fn sound_speed(&self) -> f64 {
        let t = self.mean_temperature();
        let mu = if self.inner_radius_frac < 0.3 {
            MEAN_MOLECULAR_WEIGHT_CORE
        } else {
            MEAN_MOLECULAR_WEIGHT_ENVELOPE
        };
        (ADIABATIC_INDEX * K_B * t / (mu * PROTON_MASS)).sqrt()
    }

    pub fn thermal_scale_height(&self) -> f64 {
        let t = self.mean_temperature();
        let mu = MEAN_MOLECULAR_WEIGHT_ENVELOPE;
        let r = ((self.inner_radius_frac + self.outer_radius_frac) / 2.0) * SOLAR_RADIUS;
        let g_local = G * SOLAR_MASS * (r / SOLAR_RADIUS).powi(3) / (r * r);
        K_B * t / (mu * PROTON_MASS * g_local.max(1e-10))
    }
}

pub fn core() -> StellarLayer {
    StellarLayer {
        name: "Core",
        inner_radius_frac: 0.0,
        outer_radius_frac: 0.25,
        temp_inner_k: SOLAR_CORE_TEMP,
        temp_outer_k: 7.0e6,
        density_inner: SOLAR_CORE_DENSITY,
        density_outer: 2.0e4,
    }
}

pub fn radiative_zone() -> StellarLayer {
    StellarLayer {
        name: "Radiative Zone",
        inner_radius_frac: 0.25,
        outer_radius_frac: 0.70,
        temp_inner_k: 7.0e6,
        temp_outer_k: 2.0e6,
        density_inner: 2.0e4,
        density_outer: 200.0,
    }
}

pub fn convective_zone() -> StellarLayer {
    StellarLayer {
        name: "Convective Zone",
        inner_radius_frac: 0.70,
        outer_radius_frac: 1.0,
        temp_inner_k: 2.0e6,
        temp_outer_k: SOLAR_EFFECTIVE_TEMP,
        density_inner: 200.0,
        density_outer: 2e-4,
    }
}

pub fn hydrostatic_pressure(depth_frac: f64) -> f64 {
    let rho = SOLAR_MEAN_DENSITY;
    let r = (1.0 - depth_frac) * SOLAR_RADIUS;
    let m_enclosed = SOLAR_MASS * (r / SOLAR_RADIUS).powi(3);
    G * m_enclosed * rho / r.max(1.0)
}

pub fn kelvin_helmholtz_timescale() -> f64 {
    G * SOLAR_MASS * SOLAR_MASS / (SOLAR_RADIUS * SOLAR_LUMINOSITY)
}

pub fn thermal_timescale() -> f64 {
    let e_thermal = 1.5 * K_B * SOLAR_CORE_TEMP / PROTON_MASS * SOLAR_MASS * 0.1;
    e_thermal / SOLAR_LUMINOSITY
}

pub fn nuclear_timescale() -> f64 {
    let efficiency = 0.007;
    let core_mass_frac = 0.1;
    let mc = SOLAR_MASS * core_mass_frac;
    efficiency * mc * crate::SPEED_OF_LIGHT * crate::SPEED_OF_LIGHT / SOLAR_LUMINOSITY
}

pub fn dynamical_timescale() -> f64 {
    1.0 / (G * SOLAR_MEAN_DENSITY).sqrt()
}

pub fn eddington_luminosity() -> f64 {
    let sigma_t = 6.6524e-29;
    4.0 * std::f64::consts::PI * G * SOLAR_MASS * crate::SPEED_OF_LIGHT / sigma_t
}

pub fn virial_temperature() -> f64 {
    let mu = MEAN_MOLECULAR_WEIGHT_CORE;
    G * SOLAR_MASS * mu * PROTON_MASS / (3.0 * K_B * SOLAR_RADIUS)
}

pub fn central_pressure_estimate() -> f64 {
    3.0 * G * SOLAR_MASS * SOLAR_MASS / (8.0 * std::f64::consts::PI * SOLAR_RADIUS.powi(4))
}

pub fn luminosity_from_temperature() -> f64 {
    4.0 * std::f64::consts::PI * SOLAR_RADIUS * SOLAR_RADIUS * SIGMA_SB
        * SOLAR_EFFECTIVE_TEMP.powi(4)
}

pub fn gravitational_binding_energy() -> f64 {
    3.0 * G * SOLAR_MASS * SOLAR_MASS / (5.0 * SOLAR_RADIUS)
}

pub fn surface_gravity() -> f64 {
    G * SOLAR_MASS / (SOLAR_RADIUS * SOLAR_RADIUS)
}

pub fn schwarzschild_radius() -> f64 {
    2.0 * G * SOLAR_MASS / (crate::SPEED_OF_LIGHT * crate::SPEED_OF_LIGHT)
}

pub fn roche_limit(body_density: f64) -> f64 {
    SOLAR_RADIUS * (2.0 * SOLAR_MEAN_DENSITY / body_density).powf(1.0 / 3.0)
}

pub fn jeans_mass(temperature: f64, density: f64) -> f64 {
    let cs = (K_B * temperature / (MEAN_MOLECULAR_WEIGHT_ENVELOPE * PROTON_MASS)).sqrt();
    let pi = std::f64::consts::PI;
    (pi / 6.0) * (pi * cs * cs / (G * density)).powf(1.5) * density
}