Skip to content

Commit

Permalink
Improve numerical stability for near parabolic orbits (#119)
Browse files Browse the repository at this point in the history
* Improve numerical stability for near parabolics

* make eccentric anom not wrap for hyper/parabolic

* changelog
dahlend authored Sep 24, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
1 parent 4374f5f commit cbce0f9
Showing 3 changed files with 40 additions and 15 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -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
38 changes: 28 additions & 10 deletions src/kete_core/src/elements.rs
Original file line number Diff line number Diff line change
@@ -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<f64> {
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]
15 changes: 10 additions & 5 deletions src/kete_core/src/propagation/kepler.rs
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit cbce0f9

Please sign in to comment.