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<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]
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