use crate::physics::orbit::{MERCURY_MASS, MERCURY_RADIUS};
use sciforge::hub::domain::common::constants::{G, R_GAS};
use sciforge::hub::domain::meteorology::atmosphere::scale_height;

pub const SURFACE_PRESSURE_PA: f64 = crate::EXOSPHERE_SURFACE_PRESSURE_PA;

const MU_MERCURY: f64 = G * MERCURY_MASS;
const G_SURFACE: f64 = MU_MERCURY / (MERCURY_RADIUS * MERCURY_RADIUS);

pub struct ExosphericSpecies {
    pub name: &'static str,
    pub molar_mass_kg: f64,
    pub column_density_cm2: f64,
    pub scale_height_km: f64,
}

impl ExosphericSpecies {
    pub fn number_density_at(&self, altitude_m: f64, surface_temp_k: f64) -> f64 {
        let h_s = scale_height(surface_temp_k, self.molar_mass_kg, G_SURFACE);
        let n0 = self.column_density_cm2 * 1e4 / h_s;
        n0 * (-altitude_m / h_s).exp()
    }

    pub fn jeans_escape_rate(&self, surface_temp_k: f64) -> f64 {
        let v_thermal = (2.0 * R_GAS * surface_temp_k / self.molar_mass_kg).sqrt();
        let v_esc = (2.0 * MU_MERCURY / MERCURY_RADIUS).sqrt();
        let lambda = (v_esc / v_thermal).powi(2);
        let surface_area = 4.0 * std::f64::consts::PI * MERCURY_RADIUS * MERCURY_RADIUS;
        let n0 = self.column_density_cm2 * 1e4 / (self.scale_height_km * 1000.0);
        n0 * v_thermal / (2.0 * std::f64::consts::PI.sqrt())
            * (1.0 + lambda)
            * (-lambda).exp()
            * surface_area
    }
}

pub fn sodium() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "Na",
        molar_mass_kg: 0.02299,
        column_density_cm2: crate::NA_COLUMN_DENSITY,
        scale_height_km: crate::NA_SCALE_HEIGHT_KM,
    }
}

pub fn potassium() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "K",
        molar_mass_kg: 0.03910,
        column_density_cm2: crate::K_COLUMN_DENSITY,
        scale_height_km: crate::K_SCALE_HEIGHT_KM,
    }
}

pub fn calcium() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "Ca",
        molar_mass_kg: 0.04008,
        column_density_cm2: crate::CA_COLUMN_DENSITY,
        scale_height_km: 60.0,
    }
}

pub fn magnesium() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "Mg",
        molar_mass_kg: 0.02431,
        column_density_cm2: crate::MG_COLUMN_DENSITY,
        scale_height_km: 100.0,
    }
}

pub fn helium() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "He",
        molar_mass_kg: 0.00400,
        column_density_cm2: crate::HE_COLUMN_DENSITY,
        scale_height_km: crate::HE_SCALE_HEIGHT_KM,
    }
}

pub fn hydrogen() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "H",
        molar_mass_kg: 0.00101,
        column_density_cm2: crate::H_COLUMN_DENSITY,
        scale_height_km: 500.0,
    }
}

pub fn oxygen() -> ExosphericSpecies {
    ExosphericSpecies {
        name: "O",
        molar_mass_kg: 0.01600,
        column_density_cm2: crate::O_COLUMN_DENSITY,
        scale_height_km: 150.0,
    }
}

pub fn all_species() -> Vec<ExosphericSpecies> {
    vec![
        sodium(),
        potassium(),
        calcium(),
        magnesium(),
        helium(),
        hydrogen(),
        oxygen(),
    ]
}

pub fn mean_scale_height(temp_k: f64) -> f64 {
    scale_height(temp_k, 0.02299, G_SURFACE)
}

pub fn mean_solar_irradiance() -> f64 {
    crate::SOLAR_CONSTANT_MEAN
}

pub fn effective_temperature_k() -> f64 {
    let s = crate::SOLAR_CONSTANT_MEAN;
    let albedo = crate::BOND_ALBEDO;
    let sigma = sciforge::hub::domain::common::constants::SIGMA_SB;
    ((s * (1.0 - albedo)) / (4.0 * sigma)).powf(0.25)
}
