rust_kepler_solver/
ellipse.rs

1use std::f64::consts::PI;
2
3use serde::{Deserialize, Serialize};
4
5const DELTA_THRESHOLD: f64 = 1.0e-10;
6
7fn laguerre_delta(f: f64, f_prime: f64, f_prime_prime: f64) -> f64 {
8    let n: f64 = 2.0; // N=2, modified Newton-Raphson
9    let a = (n-1.0).powi(2) * f_prime.powi(2) - n*(n-1.0)*f*f_prime_prime;
10    let mut b = f64::sqrt(a.abs());
11    b = b.abs() * f_prime.signum(); // prevent catastrophic cancellation
12    - (n*f) / (f_prime + b)
13}
14
15/// ## Example
16/// ```rs
17/// use std::f64::consts::PI;
18/// use rust_kepler_solver::ellipse::EllipseSolver;
19/// 
20/// fn example_ellipse() {
21///     let eccentricity = 1.0;
22///     let solver = EllipseSolver::new(eccentricity);
23///     println!("{}", solver.solve(1.2));
24///     println!("{}", solver.solve(PI / 4.0));
25/// }
26/// ```
27#[derive(Debug, Clone, Serialize, Deserialize)]
28pub struct EllipseSolver {
29    eccentricity: f64
30}
31
32impl EllipseSolver {
33    pub fn new(eccentricity: f64) -> Self {
34        Self { eccentricity }
35    }
36
37    /// Works for 0 < `mean_anomaly` < 2pi
38    pub fn solve(&self, mean_anomaly: f64) -> f64 {
39        // Choosing an initial seed: https://www.aanda.org/articles/aa/full_html/2022/02/aa41423-21/aa41423-21.html#S5
40        // Yes, they're actually serious about that 0.999999 thing (lmao)
41        let mut eccentric_anomaly = mean_anomaly
42            + (0.999_999 * 4.0 * self.eccentricity * mean_anomaly * (PI - mean_anomaly))
43            / (8.0 * self.eccentricity * mean_anomaly + 4.0 * self.eccentricity * (self.eccentricity - PI) + PI.powi(2));
44
45        // Iteration using laguerre method
46        // According to this 1985 paper laguerre should practially always converge (they tested it 500,000 times on different values)
47        // https://link.springer.com/content/pdf/10.1007/bf01230852.pdf
48        loop {
49            let sin_eccentric_anomaly = eccentric_anomaly.sin();
50            let cos_eccentric_anomaly = eccentric_anomaly.cos();
51            let f = mean_anomaly - eccentric_anomaly + self.eccentricity*sin_eccentric_anomaly;
52            let f_prime = -1.0 + self.eccentricity*cos_eccentric_anomaly;
53            let f_prime_prime = -self.eccentricity*sin_eccentric_anomaly;
54            let delta = laguerre_delta(f, f_prime, f_prime_prime);
55            if delta.abs() < DELTA_THRESHOLD {
56                break;
57            }
58            eccentric_anomaly += delta;
59        }
60        eccentric_anomaly
61    }
62}
63
64#[cfg(test)]
65mod test {
66    use crate::bisection::bisection;
67
68    use super::EllipseSolver;
69
70    fn solve_with_bisection(e: f64, m: f64) -> f64 {
71        // We don't care about speed here, so just use as wide a range as possible
72        let f = |eccentric_anomaly: f64| eccentric_anomaly - m - e*eccentric_anomaly.sin();
73        bisection(&f, -100000.0, 100000.0)
74    }
75
76    #[test]
77    fn test_ellipse() {
78        let eccentricites: Vec<f64> = (1..999)
79            .map(|x| x as f64 / 1000.0)
80            .collect();
81        let mean_anomalies: Vec<f64> = (0..6283) // about pi*2*100
82            .map(|x| x as f64 / 1000.0)
83            .collect();
84            
85        for e in &eccentricites {
86            let solver = EllipseSolver::new(*e);
87            for m in &mean_anomalies {
88                let expected = solve_with_bisection(*e, *m);
89                let actual = solver.solve(*m);
90                let difference = if actual != 0.0 { (expected - actual) / actual } else { expected - actual }.abs();
91                if difference > 1.0e-4 {
92                    dbg!(expected, actual, e, m);
93                    panic!()
94                }
95            }
96        }
97    }
98}
rust_kepler_solver/ellipse.rs

rust_kepler_solver/
ellipse.rs
