use crate::temporal::epoch::J2000EPOCH;
use sciforge::hub::prelude::constants::AU;

pub const NEPTUNEAXIALTILTDEG: f64 = 28.32;

pub struct SolarPosition {
    pub azimuthdeg: f64,
    pub elevationdeg: f64,
    pub distanceau: f64,
    pub direction: [f64; 3],
}

impl SolarPosition {
    pub fn compute(juliandate: f64, latdeg: f64, londeg: f64) -> Self {
        let d = juliandate - J2000EPOCH;
        let neptuneorbitalperioddays = 60182.0;
        let meananomaly = (256.228 + 360.0 * d / neptuneorbitalperioddays) % 360.0;
        let mrad = meananomaly.to_radians();
        let equationofcenter =
            0.6541 * mrad.sin() + 0.0033 * (2.0 * mrad).sin() + 0.00002 * (3.0 * mrad).sin();
        let eclipticlon = (meananomaly + equationofcenter + 180.0 + 44.971) % 360.0;
        let eclrad = eclipticlon.to_radians();
        let obliquityrad = NEPTUNEAXIALTILTDEG.to_radians();
        let declination = (obliquityrad.sin() * eclrad.sin()).asin();
        let rightascension = (obliquityrad.cos() * eclrad.sin()).atan2(eclrad.cos());
        let neptunesiderealrotperiodh = 16.11;
        let siderealtime =
            (280.16 + 360.0 / (neptunesiderealrotperiodh / 24.0) * d + londeg) % 360.0;
        let hourangle = siderealtime.to_radians() - rightascension;
        let latrad = latdeg.to_radians();
        let altitude = (latrad.sin() * declination.sin()
            + latrad.cos() * declination.cos() * hourangle.cos())
        .asin();
        let azimuth = (declination.sin() - altitude.sin() * latrad.sin())
            .atan2(altitude.cos() * latrad.cos());
        let dirx = altitude.cos() * azimuth.sin();
        let diry = altitude.sin();
        let dirz = altitude.cos() * azimuth.cos();
        let distanceau = 30.07 - 0.26 * mrad.cos();
        Self {
            azimuthdeg: azimuth.to_degrees(),
            elevationdeg: altitude.to_degrees(),
            distanceau,
            direction: [dirx, diry, dirz],
        }
    }

    pub fn distancem(&self) -> f64 {
        self.distanceau * AU
    }

    pub fn isabovehorizon(&self) -> bool {
        self.elevationdeg > 0.0
    }

    pub fn solarfluxwm2(&self) -> f64 {
        1361.0 / (self.distanceau * self.distanceau)
    }
}
