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

fn rayleighbetargb() -> [f64; 3] {
    let mh2 = 2.0 * atomic_mass(1) * 1e-3;
    let mhe = atomic_mass(2) * 1e-3;
    let mch4 = (atomic_mass(6) + 4.0 * atomic_mass(1)) * 1e-3;
    let matm = 0.80 * mh2 + 0.19 * mhe + 0.015 * mch4;
    let ndensity = N_A * *crate::ONEBARDENSITY / matm;
    let nminus1 = 1.32e-4;
    let nsqminus1sq = (2.0 * nminus1) * (2.0 * nminus1);
    let coeff = 8.0 * std::f64::consts::PI.powi(3) * nsqminus1sq / (3.0 * ndensity);
    let lambdas = [680e-9f64, 550e-9, 440e-9];
    [
        coeff / lambdas[0].powi(4),
        coeff / lambdas[1].powi(4),
        coeff / lambdas[2].powi(4),
    ]
}

pub struct AtmosphereParams {
    pub planetradiusm: f64,
    pub atmosphereheightm: f64,
    pub rayleighcoefficients: [f64; 3],
    pub rayleighscaleheightm: f64,
    pub miecoefficient: f64,
    pub miescaleheightm: f64,
    pub mieasymmetry: f64,
    pub methaneabsorption: [f64; 3],
    pub sunintensity: f64,
    pub numsamples: u32,
    pub numlightsamples: u32,
}

impl Default for AtmosphereParams {
    fn default() -> Self {
        let beta = rayleighbetargb();
        Self {
            planetradiusm: crate::NEPTUNEEQUATORIALRADIUS,
            atmosphereheightm: 600000.0,
            rayleighcoefficients: beta,
            rayleighscaleheightm: *crate::SCALEHEIGHT,
            miecoefficient: 5e-6,
            miescaleheightm: 15000.0,
            mieasymmetry: 0.75,
            methaneabsorption: [0.0, 0.001, 0.01],
            sunintensity: 22.0 / (30.07 * 30.07),
            numsamples: 16,
            numlightsamples: 8,
        }
    }
}

impl AtmosphereParams {
    pub fn rayleighdensity(&self, altitudem: f64) -> f64 {
        (-altitudem / self.rayleighscaleheightm).exp()
    }

    pub fn miedensity(&self, altitudem: f64) -> f64 {
        (-altitudem / self.miescaleheightm).exp()
    }

    pub fn scatterrayleigh(&self, costheta: f64) -> f64 {
        3.0 / (16.0 * std::f64::consts::PI) * (1.0 + costheta * costheta)
    }

    pub fn scattermie(&self, costheta: f64) -> f64 {
        let g = self.mieasymmetry;
        let g2 = g * g;
        let num = 3.0 * (1.0 - g2) * (1.0 + costheta * costheta);
        let den =
            8.0 * std::f64::consts::PI * (2.0 + g2) * (1.0 + g2 - 2.0 * g * costheta).powf(1.5);
        num / den
    }

    pub fn skycolor(&self, sundir: [f64; 3], viewdir: [f64; 3], cameraaltitudem: f64) -> [f64; 3] {
        let costheta = sundir[0] * viewdir[0] + sundir[1] * viewdir[1] + sundir[2] * viewdir[2];
        let rayleighphase = self.scatterrayleigh(costheta);
        let miephase = self.scattermie(costheta);
        const LAMBDAR: f64 = 680.0;
        const LAMBDAG: f64 = 550.0;
        const LAMBDAB: f64 = 440.0;
        let inv4r = (LAMBDAG / LAMBDAR).powi(4);
        let inv4g = 1.0;
        let inv4b = (LAMBDAG / LAMBDAB).powi(4);
        let nsteps = self.numsamples as usize;
        let atmotop = self.planetradiusm + self.atmosphereheightm;
        let raylen = atmotop - (self.planetradiusm + cameraaltitudem);
        let steplen = raylen / nsteps as f64;
        let mut color = [0.0f64; 3];
        let mut opticalr = 0.0;
        let mut opticalg = 0.0;
        let mut opticalb = 0.0;
        for i in 0..nsteps {
            let h = cameraaltitudem + steplen * (i as f64 + 0.5);
            let rhoray = self.rayleighdensity(h);
            let rhomie = self.miedensity(h);
            opticalr += (self.rayleighcoefficients[0] * inv4r * rhoray
                + self.miecoefficient * rhomie
                + self.methaneabsorption[0])
                * steplen;
            opticalg += (self.rayleighcoefficients[1] * inv4g * rhoray
                + self.miecoefficient * rhomie
                + self.methaneabsorption[1])
                * steplen;
            opticalb += (self.rayleighcoefficients[2] * inv4b * rhoray
                + self.miecoefficient * rhomie
                + self.methaneabsorption[2])
                * steplen;
            let nlight = self.numlightsamples as usize;
            let lightstep = (atmotop - (self.planetradiusm + h)) / nlight as f64;
            let mut sunr = 0.0;
            let mut sung = 0.0;
            let mut sunb = 0.0;
            for j in 0..nlight {
                let lh = h + lightstep * (j as f64 + 0.5);
                let lr = self.rayleighdensity(lh);
                let lm = self.miedensity(lh);
                sunr += (self.rayleighcoefficients[0] * inv4r * lr
                    + self.miecoefficient * lm
                    + self.methaneabsorption[0])
                    * lightstep;
                sung += (self.rayleighcoefficients[1] * inv4g * lr
                    + self.miecoefficient * lm
                    + self.methaneabsorption[1])
                    * lightstep;
                sunb += (self.rayleighcoefficients[2] * inv4b * lr
                    + self.miecoefficient * lm
                    + self.methaneabsorption[2])
                    * lightstep;
            }
            let attenr = (-sunr - opticalr).exp();
            let atteng = (-sung - opticalg).exp();
            let attenb = (-sunb - opticalb).exp();
            color[0] += (self.rayleighcoefficients[0] * inv4r * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * attenr
                * steplen;
            color[1] += (self.rayleighcoefficients[1] * inv4g * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * atteng
                * steplen;
            color[2] += (self.rayleighcoefficients[2] * inv4b * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * attenb
                * steplen;
        }
        [
            self.sunintensity * color[0],
            self.sunintensity * color[1],
            self.sunintensity * color[2],
        ]
    }
}
