use sciforge::hub::prelude::constants::N_A;
use sciforge::hub::prelude::constants::elements::atomic_mass;

#[derive(Debug, Clone)]
pub struct MolecularSpecies {
    pub name: &'static str,
    pub symbol: &'static str,
    pub molar_mass_kg_mol: f64,
    pub volume_fraction: f64,
    pub refractive_index_stp: f64,
    pub depolarization_factor: f64,
}

pub struct AtmosphereEndpoint {
    pub planet_radius_m: f64,
    pub atmosphere_height_m: f64,
    pub scale_height_m: f64,
    pub mie_coefficient: f64,
    pub mie_scale_height_m: f64,
    pub mie_asymmetry: f64,
    pub methane_absorption: [f64; 3],
    pub sun_intensity: f64,
    pub num_samples: u32,
    pub num_light_samples: u32,
    pub species: Vec<MolecularSpecies>,
    pub one_bar_density: f64,
    pub mean_molar_mass_kg_mol: f64,
    pub ammonia_absorption: f64,
}

impl AtmosphereEndpoint {
    pub fn neptune() -> Self {
        let h2_mass = 2.0 * atomic_mass(1) * 1e-3;
        let he_mass = atomic_mass(2) * 1e-3;
        let ch4_mass = (atomic_mass(6) + 4.0 * atomic_mass(1)) * 1e-3;
        let nh3_mass = (atomic_mass(7) + 3.0 * atomic_mass(1)) * 1e-3;
        let h2o_mass = (2.0 * atomic_mass(1) + atomic_mass(8)) * 1e-3;

        let species = vec![
            MolecularSpecies {
                name: "Molecular hydrogen",
                symbol: "H\u{2082}",
                molar_mass_kg_mol: h2_mass,
                volume_fraction: 0.80,
                refractive_index_stp: 1.000132,
                depolarization_factor: 0.0221,
            },
            MolecularSpecies {
                name: "Helium",
                symbol: "He",
                molar_mass_kg_mol: he_mass,
                volume_fraction: 0.19,
                refractive_index_stp: 1.000035,
                depolarization_factor: 0.025,
            },
            MolecularSpecies {
                name: "Methane",
                symbol: "CH\u{2084}",
                molar_mass_kg_mol: ch4_mass,
                volume_fraction: 0.015,
                refractive_index_stp: 1.000444,
                depolarization_factor: 0.0,
            },
            MolecularSpecies {
                name: "Ammonia",
                symbol: "NH\u{2083}",
                molar_mass_kg_mol: nh3_mass,
                volume_fraction: 1e-6,
                refractive_index_stp: 1.000376,
                depolarization_factor: 0.0,
            },
            MolecularSpecies {
                name: "Water",
                symbol: "H\u{2082}O",
                molar_mass_kg_mol: h2o_mass,
                volume_fraction: 1e-6,
                refractive_index_stp: 1.000256,
                depolarization_factor: 0.0,
            },
        ];

        let mean_molar = species
            .iter()
            .map(|s| s.molar_mass_kg_mol * s.volume_fraction)
            .sum::<f64>();

        let _ = N_A;

        Self {
            planet_radius_m: crate::NEPTUNEEQUATORIALRADIUS,
            atmosphere_height_m: 600_000.0,
            scale_height_m: *crate::SCALEHEIGHT,
            mie_coefficient: 5e-6,
            mie_scale_height_m: 15_000.0,
            mie_asymmetry: 0.75,
            methane_absorption: [0.0, 0.001, 0.01],
            sun_intensity: 22.0 / (30.07 * 30.07),
            num_samples: 16,
            num_light_samples: 8,
            species,
            one_bar_density: *crate::ONEBARDENSITY,
            mean_molar_mass_kg_mol: mean_molar,
            ammonia_absorption: 0.0,
        }
    }

    pub fn sky_luminance(&self, observer_alt: f64, cos_zenith: f64) -> f64 {
        let cos_z = cos_zenith.max(0.01);
        let rayleigh_weighted: f64 = self
            .species
            .iter()
            .map(|s| {
                let n_s =
                    self.one_bar_density * N_A / self.mean_molar_mass_kg_mol * s.volume_fraction;
                let delta = s.depolarization_factor;
                let king = (6.0 + 3.0 * delta) / (6.0 - 7.0 * delta);
                let ns = s.refractive_index_stp - 1.0;
                8.0 * std::f64::consts::PI.powi(3) * (2.0 * ns).powi(2) * king / (3.0 * n_s)
            })
            .sum();
        let altitude_factor = (-observer_alt / self.scale_height_m).exp();
        let mie_factor = (-observer_alt / self.mie_scale_height_m).exp() * self.mie_coefficient;
        let g = self.mie_asymmetry;
        let hg_phase = (1.0 - g * g)
            / (4.0 * std::f64::consts::PI * (1.0 + g * g - 2.0 * g * cos_z).powf(1.5));
        let methane_abs = self.methane_absorption[0] * 0.3
            + self.methane_absorption[1] * 0.59
            + self.methane_absorption[2] * 0.11;
        let scatter_550 = rayleigh_weighted / (550e-9_f64).powi(4) * altitude_factor;
        scatter_550 * 1e-20 * cos_z
            + mie_factor * 1e-6
            + methane_abs * 1e-8
            + hg_phase * mie_factor * 1e-4
            + self.ammonia_absorption * 1e-10
    }

    pub fn barometric_density(&self) -> f64 {
        self.one_bar_density
    }

    pub fn shell_volume(&self) -> f64 {
        4.0 / 3.0
            * std::f64::consts::PI
            * ((self.planet_radius_m + self.atmosphere_height_m).powi(3)
                - self.planet_radius_m.powi(3))
    }
}
