mercurys 0.0.2

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

pub const MERCURY_MASS: f64 = 3.3011e23;
pub const MERCURY_RADIUS: f64 = 2_439_700.0;
pub const SEMI_MAJOR_AXIS: f64 = 0.387098 * AU;
pub const ECCENTRICITY: f64 = 0.205630;
pub const INCLINATION_DEG: f64 = 7.005;
pub const LONGITUDE_ASCENDING_NODE_DEG: f64 = 48.331;
pub const ARGUMENT_PERIHELION_DEG: f64 = 29.124;
pub const MEAN_ANOMALY_J2000_DEG: f64 = 174.796;

const MU_SUN: f64 = G * SOLAR_MASS;
const MU_MERCURY: f64 = G * MERCURY_MASS;
const TWO_PI: f64 = 2.0 * std::f64::consts::PI;

pub struct MercuryOrbit {
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    pub mean_anomaly_rad: f64,
    pub eccentric_anomaly_rad: f64,
    pub true_anomaly_rad: f64,
}

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

impl MercuryOrbit {
    pub fn true_anomaly_deg(&self) -> f64 {
        self.true_anomaly_rad.to_degrees()
    }
    pub fn new() -> Self {
        let m0 = MEAN_ANOMALY_J2000_DEG.to_radians();
        let mut s = Self {
            semi_major_axis_m: SEMI_MAJOR_AXIS,
            eccentricity: ECCENTRICITY,
            mean_anomaly_rad: m0,
            eccentric_anomaly_rad: m0,
            true_anomaly_rad: 0.0,
        };
        s.solve_kepler();
        s
    }

    pub fn at_mean_anomaly(mean_anomaly_rad: f64) -> Self {
        let mut s = Self {
            semi_major_axis_m: SEMI_MAJOR_AXIS,
            eccentricity: ECCENTRICITY,
            mean_anomaly_rad,
            eccentric_anomaly_rad: mean_anomaly_rad,
            true_anomaly_rad: 0.0,
        };
        s.solve_kepler();
        s
    }

    pub fn orbital_period_s(&self) -> f64 {
        kepler_period(self.semi_major_axis_m, MU_SUN)
    }

    pub fn orbital_period_days(&self) -> f64 {
        self.orbital_period_s() / crate::SECONDS_PER_EARTH_DAY
    }

    pub fn mean_motion(&self) -> f64 {
        TWO_PI / self.orbital_period_s()
    }

    pub fn velocity_at_distance(&self, r: f64) -> f64 {
        vis_viva(MU_SUN, r, self.semi_major_axis_m).abs().sqrt()
    }

    pub fn perihelion_m(&self) -> f64 {
        periapsis(self.semi_major_axis_m, self.eccentricity)
    }

    pub fn aphelion_m(&self) -> f64 {
        apoapsis(self.semi_major_axis_m, self.eccentricity)
    }

    pub fn specific_orbital_energy(&self) -> f64 {
        orbital_energy(MU_SUN, self.semi_major_axis_m)
    }

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

    pub fn escape_velocity_at_surface() -> f64 {
        escape_velocity(MU_MERCURY, MERCURY_RADIUS)
    }

    pub fn mean_orbital_velocity(&self) -> f64 {
        TWO_PI * self.semi_major_axis_m / self.orbital_period_s()
    }

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

    pub fn current_radius(&self) -> f64 {
        true_anomaly_to_radius(
            self.semi_major_axis_m,
            self.eccentricity,
            self.true_anomaly_rad,
        )
    }

    pub fn solar_irradiance(&self) -> f64 {
        let r_au = self.current_radius() / AU;
        sciforge::hub::domain::meteorology::radiation::solar_constant() / (r_au * r_au)
    }

    pub fn position(&self) -> (f64, f64, f64) {
        let r = self.current_radius();
        let nu = self.true_anomaly_rad;
        let omega = ARGUMENT_PERIHELION_DEG.to_radians();
        let big_omega = LONGITUDE_ASCENDING_NODE_DEG.to_radians();
        let inc = INCLINATION_DEG.to_radians();

        let u = omega + nu;
        let x = r * (big_omega.cos() * u.cos() - big_omega.sin() * u.sin() * inc.cos());
        let y = r * (big_omega.sin() * u.cos() + big_omega.cos() * u.sin() * inc.cos());
        let z = r * (u.sin() * inc.sin());
        (x, y, z)
    }

    fn solve_kepler(&mut self) {
        let m = self.mean_anomaly_rad;
        let e = self.eccentricity;
        let mut ea = m;
        for _ in 0..15 {
            ea = ea - (ea - e * ea.sin() - m) / (1.0 - e * ea.cos());
        }
        self.eccentric_anomaly_rad = ea;
        let half_ea = ea / 2.0;
        self.true_anomaly_rad =
            2.0 * ((1.0 + e).sqrt() * half_ea.sin()).atan2((1.0 - e).sqrt() * half_ea.cos());
    }

    pub fn step(&mut self, dt_s: f64) {
        self.mean_anomaly_rad += self.mean_motion() * dt_s;
        self.mean_anomaly_rad %= TWO_PI;
        self.solve_kepler();
    }

    pub fn gr_perihelion_precession_rad_per_orbit(&self) -> f64 {
        let a = self.semi_major_axis_m;
        let e = self.eccentricity;
        let c = sciforge::hub::domain::common::constants::C;
        6.0 * std::f64::consts::PI * MU_SUN / (a * (1.0 - e * e) * c * c)
    }
}