use crate::temporal::calendar::SECONDSPERDAY;
use sciforge::hub::domain::astronomy::celestial::gravitational_force;
use sciforge::hub::domain::astronomy::orbits::{
    angular_momentum, apoapsis, escape_velocity, kepler_period, orbital_energy, periapsis, vis_viva,
};
use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};

pub const SEMIMAJORAXIS: f64 = 30.07 * AU;
pub const ECCENTRICITY: f64 = 0.008678;
pub const INCLINATIONDEG: f64 = 1.77;
pub const LONGITUDEASCENDINGNODEDEG: f64 = 131.784;
pub const ARGUMENTPERIHELIONDEG: f64 = 276.336;

const MUSUN: f64 = G * SOLAR_MASS;

pub struct NeptuneOrbit {
    pub semimajoraxism: f64,
    pub eccentricity: f64,
    pub trueanomalyrad: f64,
}

impl Default for NeptuneOrbit {
    fn default() -> Self {
        Self::new()
    }
}

impl NeptuneOrbit {
    pub fn new() -> Self {
        Self {
            semimajoraxism: SEMIMAJORAXIS,
            eccentricity: ECCENTRICITY,
            trueanomalyrad: 0.0,
        }
    }

    pub fn orbitalperiods(&self) -> f64 {
        kepler_period(self.semimajoraxism, MUSUN)
    }

    pub fn orbitalperioddays(&self) -> f64 {
        self.orbitalperiods() / SECONDSPERDAY
    }

    pub fn velocityatdistance(&self, r: f64) -> f64 {
        vis_viva(MUSUN, r, self.semimajoraxism)
    }

    pub fn perihelionm(&self) -> f64 {
        periapsis(self.semimajoraxism, self.eccentricity)
    }

    pub fn aphelionm(&self) -> f64 {
        apoapsis(self.semimajoraxism, self.eccentricity)
    }

    pub fn specific_orbital_energy(&self) -> f64 {
        orbital_energy(MUSUN, self.semimajoraxism)
    }

    pub fn specific_angular_momentum(&self) -> f64 {
        angular_momentum(MUSUN, self.semimajoraxism, self.eccentricity)
    }

    pub fn escape_velocity_at_surface() -> f64 {
        escape_velocity(G * crate::NEPTUNEMASS, crate::NEPTUNEEQUATORIALRADIUS)
    }

    pub fn gravitational_force_sun(&self) -> f64 {
        let r = self.current_radius();
        gravitational_force(SOLAR_MASS, crate::NEPTUNEMASS, r)
    }

    pub fn current_radius(&self) -> f64 {
        let a = self.semimajoraxism;
        let e = self.eccentricity;
        a * (1.0 - e * e) / (1.0 + e * self.trueanomalyrad.cos())
    }

    pub fn mean_orbital_velocity(&self) -> f64 {
        let period = self.orbitalperiods();
        2.0 * std::f64::consts::PI * self.semimajoraxism / period
    }
}
