1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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);
}
}
|
