From cbce0f93d12bed8b73471e5c4ec9b25779b37563 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Tue, 24 Sep 2024 15:35:46 -0700 Subject: [PATCH] Improve numerical stability for near parabolic orbits (#119) * Improve numerical stability for near parabolics * make eccentric anom not wrap for hyper/parabolic * changelog --- CHANGELOG.md | 2 ++ src/kete_core/src/elements.rs | 38 ++++++++++++++++++------- src/kete_core/src/propagation/kepler.rs | 15 ++++++---- 3 files changed, 40 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ec4e15a..de60544 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -30,6 +30,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Orbital elements calculations for True Anomaly and Eccentric Anomaly were numerically + unstable with nearly parabolic orbits. - Time scaling bugs when loading times from both Horizons and the MPC were fixed, these were causing offsets of about a minute in the loaded states. However it was shifting both the perihelion time and the epoch time by the same amount, leading to only minor diff --git a/src/kete_core/src/elements.rs b/src/kete_core/src/elements.rs index 5eda62c..e687dc7 100644 --- a/src/kete_core/src/elements.rs +++ b/src/kete_core/src/elements.rs @@ -2,7 +2,7 @@ //! This allows conversion to and from cometary orbital elements to [`State`]. use crate::constants::GMS_SQRT; use crate::prelude::{Desig, Frame, NeosResult, State}; -use crate::propagation::{compute_eccentric_anomaly, compute_true_anomaly}; +use crate::propagation::{compute_eccentric_anomaly, compute_true_anomaly, PARABOLIC_ECC_LIMIT}; use nalgebra::Vector3; use std::f64::consts::TAU; @@ -117,7 +117,7 @@ impl CometElements { }; let peri_time: f64 = { - if (ecc - 1.0).abs() < 1e-8 { + if (ecc - 1.0).abs() < PARABOLIC_ECC_LIMIT { // Parabolic let mut true_anomaly = ecc_vec.angle(pos); if vp_mag.is_sign_negative() { @@ -126,7 +126,7 @@ impl CometElements { let d = (true_anomaly / 2.0).tan(); let dt = (2f64.sqrt() * peri_dist.powf(1.5) / GMS_SQRT) * (d + d.powi(3) / 3.0); epoch - dt - } else if ecc < 1e-8 { + } else if ecc < 1e-6 { let semi_major = (2.0 / p_mag - v_mag2).recip(); let mean_motion = semi_major.abs().powf(-1.5) * GMS_SQRT; // for circular cases mean_anomaly == true_anomaly @@ -188,8 +188,8 @@ impl CometElements { /// Convert orbital elements into a cartesian coordinate position and velocity. /// Units are in AU and AU/Day. fn to_pos_vel(&self) -> NeosResult<[[f64; 3]; 2]> { - let elliptical = self.eccentricity < 0.9999; - let hyperbolic = self.eccentricity > 1.0001; + let elliptical = self.eccentricity < 1.0 - PARABOLIC_ECC_LIMIT; + let hyperbolic = self.eccentricity > 1.0 + PARABOLIC_ECC_LIMIT; let parabolic = !elliptical & !hyperbolic; // these handle parabolic in a non-standard way which allows for the @@ -262,17 +262,22 @@ impl CometElements { Ok([pos, vel]) } + /// Compute the eccentric anomaly for the cometary elements. pub fn eccentric_anomaly(&self) -> NeosResult { - compute_eccentric_anomaly(self.eccentricity, self.mean_anomaly(), self.peri_dist) - .map(|x| x.rem_euclid(TAU)) + compute_eccentric_anomaly(self.eccentricity, self.mean_anomaly(), self.peri_dist).map(|x| { + match self.eccentricity { + ecc if ecc > 1.0 - PARABOLIC_ECC_LIMIT => x, + _ => x.rem_euclid(TAU), + } + }) } /// Compute the semi major axis in AU. /// NAN is returned if the orbit is parabolic. pub fn semi_major(&self) -> f64 { match self.eccentricity { - ecc if ((ecc - 1.0).abs() <= 1e-5) => f64::NAN, + ecc if ((ecc - 1.0).abs() <= PARABOLIC_ECC_LIMIT) => f64::NAN, ecc => self.peri_dist / (1.0 - ecc), } } @@ -290,7 +295,7 @@ impl CometElements { /// Compute the mean motion in radians per day. pub fn mean_motion(&self) -> f64 { match self.eccentricity { - ecc if ((ecc - 1.0).abs() <= 1e-5) => { + ecc if ((ecc - 1.0).abs() <= PARABOLIC_ECC_LIMIT) => { GMS_SQRT * 1.5 / 2f64.sqrt() / self.peri_dist.powf(1.5) } _ => GMS_SQRT / self.semi_major().abs().powf(1.5), @@ -302,7 +307,7 @@ impl CometElements { let mm = self.mean_motion(); let mean_anomaly = (self.epoch - self.peri_time) * mm; match self.eccentricity { - ecc if ecc < 0.9999 => mean_anomaly.rem_euclid(TAU), + ecc if ecc < 1.0 - PARABOLIC_ECC_LIMIT => mean_anomaly.rem_euclid(TAU), _ => mean_anomaly, } } @@ -349,6 +354,19 @@ mod tests { }; assert!((elem.true_anomaly().unwrap() - 1.198554792).abs() < 1e-6); assert!(elem.to_pos_vel().is_ok()); + + let elem = CometElements { + desig: Desig::Empty, + frame: Frame::Ecliptic, + epoch: 2455562.5, + eccentricity: 0.99999, + inclination: 2.792526803, + lon_of_ascending: 0.349065850, + peri_time: 2455369.7, + peri_arg: -0.8726646259, + peri_dist: 0.5, + }; + assert!((elem.true_anomaly().unwrap() - 2.6071638616282553).abs() < 1e-6); } #[test] diff --git a/src/kete_core/src/propagation/kepler.rs b/src/kete_core/src/propagation/kepler.rs index 191a054..b527992 100644 --- a/src/kete_core/src/propagation/kepler.rs +++ b/src/kete_core/src/propagation/kepler.rs @@ -13,6 +13,9 @@ use core::f64; use nalgebra::{ComplexField, Vector3}; use std::f64::consts::TAU; +/// How close to ecc=1 do we assume the orbit is parabolic +pub const PARABOLIC_ECC_LIMIT: f64 = 1e-3; + /// Compute the eccentric anomaly for all orbital classes. /// /// # Arguments @@ -30,13 +33,13 @@ pub fn compute_eccentric_anomaly(ecc: f64, mean_anom: f64, peri_dist: f64) -> Ne "Eccentricity must be greater than 0".into(), ))?, ecc if ecc < 1e-6 => Ok(mean_anom), - ecc if ecc < 0.9999 => { + ecc if ecc < 1.0 - PARABOLIC_ECC_LIMIT => { // Elliptical let f = |ecc_anom: f64| -ecc * ecc_anom.sin() + ecc_anom - mean_anom.rem_euclid(TAU); let d = |ecc_anom: f64| 1.0 - 1.0 * ecc * ecc_anom.cos(); Ok(newton_raphson(f, d, mean_anom.rem_euclid(TAU), 1e-11)?.rem_euclid(TAU)) } - ecc if ecc < 1.0001 => { + ecc if ecc < 1.0 + PARABOLIC_ECC_LIMIT => { // Parabolic // Find the zero point of -mean_anom + peri_dist * ecc_anom + ecc_anom.powi(3) / 6.0 // this is a simple cubic equation which can be solved analytically. @@ -68,7 +71,9 @@ pub fn compute_true_anomaly(ecc: f64, mean_anom: f64, peri_dist: f64) -> NeosRes let ecc_anom = compute_eccentric_anomaly(ecc, mean_anom, peri_dist)?; let anom = match ecc { - e if (e - 1.0).abs() < 1e-6 => (ecc_anom / (2.0 * peri_dist).sqrt()).atan() * 2.0, + e if (e - 1.0).abs() < PARABOLIC_ECC_LIMIT => { + (ecc_anom / (2.0 * peri_dist).sqrt()).atan() * 2.0 + } e if e < 1.0 => (((1.0 + ecc) / (1.0 - ecc)).sqrt() * (ecc_anom / 2.0).tan()).atan() * 2.0, e if e > 1.0 => { ((-(1.0 + ecc) / (1.0 - ecc)).sqrt() * (ecc_anom / 2.0).tanh()).atan() * 2.0 @@ -97,12 +102,12 @@ pub fn eccentric_anomaly_from_true(ecc: f64, true_anom: f64, peri_dist: f64) -> "Eccentricity must be greater than 0".into(), ))?, e if e < 1e-6 => true_anom, - e if e < 0.9999 => { + e if e < 1.0 - PARABOLIC_ECC_LIMIT => { let (sin_true, cos_true) = true_anom.sin_cos(); let t = (1.0 + ecc * cos_true).recip(); ((1.0 - ecc.powi(2)).sqrt() * sin_true * t).atan2((ecc + cos_true) * t) } - e if e < 1.0001 => (0.5 * true_anom).tan() * (2.0 * peri_dist).sqrt(), + e if e < 1.0 + PARABOLIC_ECC_LIMIT => (0.5 * true_anom).tan() * (2.0 * peri_dist).sqrt(), _ => { let v = (-(1.0 - ecc) / (1.0 + ecc)).sqrt(); ((true_anom * 0.5).tan() * v).atanh() * 2.0