rust_kepler_solver/
hyperbola.rs

1use lazy_static::lazy_static;
2use serde::{Deserialize, Serialize};
3
4const CUBIC_DELTA_THRESHOLD: f64 = 1.0e-6;
5
6lazy_static! {
7    // From eq. 4 in the B. Wu et all paper
8    // Max value in each interval
9    static ref PADE_ECCENTRIC_ANOMALY_THRESHOLDS: [f64; 15] = [
10        40.0 / 8.0,
11        38.0 / 8.0,
12        34.0 / 8.0,
13        30.0 / 8.0,
14        26.0 / 8.0,
15        22.0 / 8.0,
16        18.0 / 8.0,
17        15.0 / 8.0,
18        13.0 / 8.0,
19        11.0 / 8.0,
20        9.0  / 8.0,
21        7.0  / 8.0,
22        5.0  / 8.0,
23        3.0  / 8.0,
24        29.0 / 200.0];
25
26    static ref PADE_ORDERS: [f64; 15] = [
27        10.0 / 2.0,
28        9.0 / 2.0,
29        8.0 / 2.0,
30        7.0 / 2.0,
31        6.0 / 2.0,
32        5.0 / 2.0,
33        8.0 / 4.0,
34        7.0 / 4.0,
35        6.0 / 4.0,
36        5.0 / 4.0,
37        4.0 / 4.0,
38        3.0 / 4.0,
39        2.0 / 4.0,
40        1.0 / 4.0,
41        0.0 / 4.0];
42}
43
44fn hyperbolic_kepler_equation(eccentricity: f64, eccentric_anomaly: f64) -> f64 {
45    eccentricity * f64::sinh(eccentric_anomaly) - eccentric_anomaly
46}
47
48fn pade_approximation(ec: f64, mh: f64, a: f64) -> [f64; 4] {
49    let ex = f64::exp(a);
50    let enx = f64::exp(-a);
51    let sa = (ex-enx) / 2.0; // sinh(a)
52    let ca = (ex+enx) / 2.0; // cosh(a)
53    let d1 = ca.powi(2) + 3.0;
54    let d2 = sa.powi(2) + 4.0;
55    let p1 = ca * (3.0 * ca.powi(2) + 17.0) / (5.0 * d1);
56    let p2 = sa * (3.0 * sa.powi(2) + 28.0) / (20.0 * d2);
57    let p3 = ca * (ca.powi(2) + 27.0) / (60.0 * d1);
58    let q1 = -2.0 * ca * sa / (5.0 * d1);
59    let q2 = (sa.powi(2) - 4.0) / (20.0 * d2);
60    let coefficient_f3 = ec * p3 - q2;
61    let coefficient_f2 = ec * p2 - (mh + a) * q2 - q1;
62    let coefficient_f1 = ec * p1 - (mh + a) * q1 - 1.0;
63    let coefficient_f0 = ec * sa - mh - a;
64    [coefficient_f3,coefficient_f2, coefficient_f1, coefficient_f0]
65}
66
67fn solve_cubic(coefficients: [f64; 4], mh: f64, ec: f64) -> f64 {
68    let mut x = mh / (ec - 1.0); // starting value from series expansion of HKE
69    loop {
70        // halley's method
71        let f = ((coefficients[0]*x + coefficients[1])*x + coefficients[2])*x + coefficients[3];
72        let f_prime = (3.0*coefficients[0]*x + 2.0*coefficients[1])*x + coefficients[2];
73        let f_prime_prime = 6.0*coefficients[0]*x + 2.0*coefficients[1];
74        let delta = -2.0*f*f_prime / (2.0*f_prime.powi(2) - f*f_prime_prime);
75        if delta.abs() < CUBIC_DELTA_THRESHOLD {
76            break;
77        }
78        x += delta;
79    }
80    x
81}
82
83/// ## Example
84/// ```rs
85/// use rust_kepler_solver::hyperbola::HyperbolaSolver;
86///
87/// fn example_hyperbola() {
88///     let eccentricity = 1.0;
89///     let solver = HyperbolaSolver::new(eccentricity);
90///     println!("{}", solver.solve(1.2));
91///     println!("{}", solver.solve(100.0));
92/// }
93/// ```
94#[derive(Debug, Clone, Serialize, Deserialize)]
95pub struct HyperbolaSolver {
96    eccentricity: f64,
97    pade_mean_anomaly_thresholds: [f64; 15],
98}
99
100impl HyperbolaSolver {
101    pub fn new(eccentricity: f64) -> Self {
102        let pade_mean_anomaly_thresholds = PADE_ECCENTRIC_ANOMALY_THRESHOLDS.map(|eccentric_anomaly| hyperbolic_kepler_equation(eccentricity, eccentric_anomaly));
103        Self { eccentricity, pade_mean_anomaly_thresholds }
104    }
105
106    /// Works with all values of mean anomaly 0 to infinity
107    pub fn solve(&self, mean_anomaly: f64) -> f64 {
108        // Solver assumes mean anomaly > 0
109        // The equation is symmetric, so for mean anomaly < 0, we just flip the sign o the output
110        let ec = self.eccentricity;
111        let mh = mean_anomaly.abs();
112
113        let f0 = if mh <= self.pade_mean_anomaly_thresholds[0] {
114            // For mh < 5 we use a 'piecewise pade approximation' to get the starting estimate
115            let mut i = 0;
116            while i < self.pade_mean_anomaly_thresholds.len()-1 && mh < self.pade_mean_anomaly_thresholds[i+1] {
117                i += 1;
118            }
119            let a = PADE_ORDERS[i];
120            let coefficients = pade_approximation(ec, mh, a);
121
122            solve_cubic(coefficients, mh, ec) + a
123
124        } else {
125            // For mh >= 5, we can use this... thing that I copied from the above paper
126            // I have no idea how it works, but it works very very very well
127            let fa = f64::ln(2.0 * mh / ec);
128            let ca = 0.5 * (2.0 * mh / ec + ec / (2.0 * mh));
129            let sa = 0.5 * (2.0 * mh / ec - ec / (2.0 * mh));
130            let top = 6.0 * (ec.powi(2) / (4.0 * mh) + fa) / (ec * ca - 1.0) 
131                + 3.0 * (ec * sa / (ec * ca - 1.0)) * ((ec.powi(2) / (4.0 * mh) + fa) / (ec * ca - 1.0)).powi(2);
132            let bottom = 6.0 + 6.0 * (ec * sa / (ec * ca - 1.0)) * ((ec.powi(2) / (4.0 * mh) + fa) / (ec * ca - 1.0))
133                + (ec * ca / (ec * ca - 1.0)) * ((ec.powi(2) / (4.0 * mh) + fa) / (ec * ca - 1.0)).powi(2);
134            let delta = top / bottom;
135            fa + delta
136        };
137
138        // Halley method
139        let f = ec * f0.sinh() - f0 - mh;
140        let f_prime = ec * f0.cosh() - 1.0;
141        let f_prime_prime = f_prime + 1.0;
142        let f1 = f0 - (2.0 * f / f_prime) / (2.0 - f * f_prime_prime / f_prime.powi(2));
143
144        f1 * mean_anomaly.signum()
145    }
146}
147
148#[cfg(test)]
149mod test {
150    use crate::bisection::bisection;
151
152    use super::HyperbolaSolver;
153
154    fn solve_with_bisection(e: f64, m: f64) -> f64 {
155        // We don't care about speed here, so just use as wide a range as possible
156        let f = |eccentric_anomaly: f64| e * f64::sinh(eccentric_anomaly) - eccentric_anomaly - m;
157        bisection(&f, -100000.0, 100000.0)
158    }
159
160    #[test]
161    fn test_hyperbola() {
162        let eccentricites: Vec<f64> = (1..999)
163            .map(|x| 1.0 + f64::powi(x as f64, 2) / 1000.0)
164            .collect();
165        let mean_anomalies: Vec<f64> = (0..10000)
166            .map(|x| f64::powi(x as f64, 2) / 10000.0)
167            .collect();
168            
169        for e in &eccentricites {
170            let solver = HyperbolaSolver::new(*e);
171            for m in &mean_anomalies {
172                let expected = solve_with_bisection(*e, *m);
173                let actual = solver.solve(*m);
174                let difference = if actual.abs() < 1.0e-5 { expected - actual } else { (expected - actual) / actual }.abs();
175                if difference > 1.0e-4 {
176                    dbg!(expected, actual, difference, e, m);
177                    panic!()
178                }
179            }
180        }
181    }
182}
rust_kepler_solver/hyperbola.rs

rust_kepler_solver/
hyperbola.rs
