use sciforge::hub::domain::common::constants::DEGREE_TO_RAD;

pub const SIDEREALDAYS: f64 = 57996.0;
pub const AXIALTILTDEG: f64 = 28.32;
pub const AXIALTILTRAD: f64 = AXIALTILTDEG * DEGREE_TO_RAD;
pub const PRECESSIONPERIODYEARS: f64 = 26000.0;

pub struct NeptuneRotation {
    pub angularvelocityrads: f64,
    pub axialtiltrad: f64,
}

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

impl NeptuneRotation {
    pub fn new() -> Self {
        Self {
            angularvelocityrads: 2.0 * std::f64::consts::PI / SIDEREALDAYS,
            axialtiltrad: AXIALTILTRAD,
        }
    }

    pub fn surfacevelocityatlatitude(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        self.angularvelocityrads * crate::NEPTUNEEQUATORIALRADIUS * latrad.cos()
    }

    pub fn centripetalaccelerationatlatitude(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        let omega = self.angularvelocityrads;
        omega * omega * crate::NEPTUNEEQUATORIALRADIUS * latrad.cos()
    }

    pub fn coriolisparameter(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        2.0 * self.angularvelocityrads * latrad.sin()
    }

    pub fn momentofinertia(&self) -> f64 {
        0.23 * crate::NEPTUNEMASS * crate::NEPTUNEEQUATORIALRADIUS * crate::NEPTUNEEQUATORIALRADIUS
    }

    pub fn rotationalkineticenergy(&self) -> f64 {
        0.5 * self.momentofinertia() * self.angularvelocityrads.powi(2)
    }

    pub fn angular_momentum(&self) -> f64 {
        self.momentofinertia() * self.angularvelocityrads
    }

    pub fn precession_rate_rad_per_year(&self) -> f64 {
        2.0 * std::f64::consts::PI / PRECESSIONPERIODYEARS
    }

    pub fn effective_tilt_rad(&self, juliancenturiessincej2000: f64) -> f64 {
        let rate = self.precession_rate_rad_per_year() / 100.0;
        self.axialtiltrad + rate * juliancenturiessincej2000 * 0.001
    }

    pub fn day_length_variation_due_to_tilt(&self, dayofyear: u32, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        let orbitalperioddays = 60190.0;
        let declination = self.axialtiltrad
            * (2.0 * std::f64::consts::PI * (dayofyear as f64 - 81.0) / orbitalperioddays).sin();
        let coshourangle = -(latrad.tan() * declination.tan());
        if coshourangle > 1.0 {
            return 0.0;
        }
        if coshourangle < -1.0 {
            return SIDEREALDAYS / 3600.0;
        }
        2.0 * coshourangle.acos() * (SIDEREALDAYS / 3600.0) / (2.0 * std::f64::consts::PI)
    }
}
