//! Kepler orbits, and how long trips between them take.
//!
//! Distances are in AU, times in days, masses in suns. Orbits lie flat in
//! the plane of the system (no inclination) and run anticlockwise.

use std::f64::consts::TAU;

use super::Vec2;

/// Days in a year for a body 1 AU from a star of one solar mass.
const YEAR: f64 = 365.25;
/// An Earth mass, in suns.
pub const EARTH_MASS: f64 = 3.003e-6;

#[derive(Clone, Copy, Debug)]
pub struct Orbit {
    /// Semi-major axis.
    pub a: f64,
    /// Eccentricity: 0 is a circle.
    pub e: f64,
    /// Angle from the x axis to the closest point of the orbit.
    pub periapsis: f64,
    /// Where along the orbit the body is at day 0, as a mean anomaly.
    pub mean_at_epoch: f64,
    /// Days per revolution.
    pub period: f64,
}

impl Orbit {
    /// An orbit around a body of `parent_mass` suns, with its period from
    /// Kepler's third law.
    pub fn new(a: f64, e: f64, periapsis: f64, mean_at_epoch: f64, parent_mass: f64) -> Orbit {
        let period = YEAR * (a.powi(3) / parent_mass).sqrt();
        Orbit { a, e, periapsis, mean_at_epoch, period }
    }

    /// Position relative to the parent on `day`.
    pub fn position(&self, day: f64) -> Vec2 {
        let mean = self.mean_at_epoch + TAU * day / self.period;
        self.at_eccentric_anomaly(eccentric_anomaly(mean, self.e))
    }

    /// `n` points evenly spaced around the ellipse, for drawing it.
    pub fn outline(&self, n: usize) -> impl Iterator<Item = Vec2> + '_ {
        (0..n).map(move |i| self.at_eccentric_anomaly(TAU * i as f64 / n as f64))
    }

    /// Closest and farthest distances from the parent.
    pub fn periapsis_distance(&self) -> f64 {
        self.a * (1.0 - self.e)
    }

    pub fn apoapsis_distance(&self) -> f64 {
        self.a * (1.0 + self.e)
    }

    fn at_eccentric_anomaly(&self, ea: f64) -> Vec2 {
        let x = self.a * (ea.cos() - self.e);
        let y = self.a * (1.0 - self.e * self.e).sqrt() * ea.sin();
        Vec2::new(x, y).rotated(self.periapsis)
    }
}

/// Solves Kepler's equation, M = E - e sin E, for E by Newton's method.
fn eccentric_anomaly(mean: f64, e: f64) -> f64 {
    let m = mean.rem_euclid(TAU);
    let mut ea = if e > 0.8 { std::f64::consts::PI } else { m };
    for _ in 0..12 {
        let step = (ea - e * ea.sin() - m) / (1.0 - e * ea.cos());
        ea -= step;
        if step.abs() < 1e-12 {
            break;
        }
    }
    ea
}

/// Days a trip of 1 AU takes. Trips are skewed well away from realism, so
/// that a far planet costs only a few days more than a near one: days grow
/// with the cube root of the distance. A moon is hours away, most planets
/// two to five days, and the far edge of the widest systems about ten.
pub const AU_TRIP: f64 = 2.25;

/// The shortest trip, however close the destination: undocking, getting
/// clear and settling into the new place.
pub const MIN_TRIP: f64 = 0.25;

pub fn trip_days(distance: f64) -> f64 {
    (AU_TRIP * distance.cbrt()).max(MIN_TRIP)
}

/// How many days a trip from `from` (where you are now, on `day`) to a
/// destination that moves along `destination(day)` takes. Aims for where
/// the destination will be on arrival, not where it is now.
pub fn plan_trip(from: Vec2, day: f64, destination: impl Fn(f64) -> Vec2) -> f64 {
    let mut days = trip_days(from.distance(destination(day)));
    for _ in 0..8 {
        days = trip_days(from.distance(destination(day + days)));
    }
    days
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn keplers_third_law() {
        let earthlike = Orbit::new(1.0, 0.0, 0.0, 0.0, 1.0);
        assert!((earthlike.period - 365.25).abs() < 1e-9);
        let far = Orbit::new(4.0, 0.0, 0.0, 0.0, 1.0);
        assert!(
            (far.period - 8.0 * 365.25).abs() < 1e-9,
            "four times out takes eight times as long"
        );
    }

    #[test]
    fn orbits_come_round_again() {
        let o = Orbit::new(2.3, 0.2, 1.0, 0.4, 0.8);
        let (p, q) = (o.position(17.0), o.position(17.0 + o.period));
        assert!(p.distance(q) < 1e-9);
    }

    #[test]
    fn ellipses_swing_between_periapsis_and_apoapsis() {
        let o = Orbit::new(2.0, 0.25, 0.7, 0.0, 1.0);
        let at_start = o.position(0.0).length();
        assert!((at_start - o.periapsis_distance()).abs() < 1e-9, "mean anomaly 0 is periapsis");
        let half = o.position(o.period / 2.0).length();
        assert!((half - o.apoapsis_distance()).abs() < 1e-9);
        let p = o.position(0.0);
        let angle = p.y.atan2(p.x);
        assert!((angle - 0.7).abs() < 1e-9, "periapsis points where it says");
    }

    #[test]
    fn keplers_second_law_bodies_hurry_near_periapsis() {
        let o = Orbit::new(1.0, 0.3, 0.0, 0.0, 1.0);
        let near = o.position(0.0).distance(o.position(1.0));
        let far = o.position(o.period / 2.0).distance(o.position(o.period / 2.0 + 1.0));
        assert!(near > far * 1.5);
    }

    #[test]
    fn trips_aim_ahead_of_a_moving_target() {
        let target = Orbit::new(1.0, 0.0, 0.0, 0.0, 1.0);
        let from = Vec2::new(5.0, 0.0);
        let days = plan_trip(from, 0.0, |d| target.position(d));
        let arrival = target.position(days);
        assert!((trip_days(from.distance(arrival)) - days).abs() < 1e-6);
    }

    #[test]
    fn trips_take_days_not_weeks() {
        assert!(trip_days(0.003) < 0.5, "a moon is hours away");
        assert!((2.0..3.0).contains(&trip_days(1.0)));
        assert!(trip_days(10.0) < 5.0);
        assert!(trip_days(90.0) < 11.0, "even the far edge of a wide system");
        assert_eq!(trip_days(1e-6), MIN_TRIP);
        // Far places cost more, but not much more.
        assert!((trip_days(27.0) / trip_days(1.0) - 3.0).abs() < 1e-9);
    }
}