1 files changed, 164 insertions, 0 deletions
diff --git a/src/space/orbit.rs b/src/space/orbit.rs
new file mode 100644
index 0000000..0bcea5e
--- /dev/null
+++ b/src/space/orbit.rs
@@ -0,0 +1,164 @@
+//! 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);
+    }
+}