Pretty simple question here. I'm calculating time of flight between two true anomalies by converting those true anomalies to mean anomalies and using the following equation:

$\Delta t = t_2-t_1 = \sqrt{\frac{a^3}{\mu}}(M_2-M_1)$

In comparing this to a simple inertial propagation method (Cowell's Method) I'm getting good agreement for eccentricities lower than 0.5, but massive diagreement (on the order of 50% plus) when I test eccentricities greater than 0.5.

I know that if one were using a method that goes from $M -> E$ that solving the transcendental equation $M = E - e\sin(E)$ using Newton's method breaks down for high eccentricities. However, I'm not doing this... I'm going from true anomaly ($\theta$) -> eccentric anomaly ($E$) -> mean anomaly ($M$). It's all algebraic with that particular path.

Just wanted to know if anybody had experience or insight into whether Kepler's method itself breaks down at high eccentricities, or if I should be looking at my propagator?

Thanks! -Dave

  • 3
    $\begingroup$ Yes it does. The problem ultimately traces back to geometry. The mean anomaly is defined in terms of the center of the orbit, but the center of a parabola does not exist. A full answer would describe how to work around this failure, hopefully I or someone else will find time to derive it. $\endgroup$ Jul 23, 2020 at 1:08
  • 1
    $\begingroup$ Are you sure the resuts of Numerov's method (also called Cowell's method) are correct? The method is derived using a fifth order taylor expansion. To get an analytical exact method an infinite number of taylor series terms should be used. So Kepler's method does not break down at high eccentricities but Cowell's Method fails. $\endgroup$
    – Uwe
    Jul 23, 2020 at 13:19
  • $\begingroup$ @Uwe The issue is that things start to diverge before the inertial propagator should fail. I've run side-by-side tests of the propagator, and for less than one revolution, I usually get an error from the propagator (compared to a perfectly unperturbed orbit) on the order of a few meters even for eccentricities close to 1; eccentricities that are greater than where I'm seeing divergence between Kepler's and the Numerical time of flight. $\endgroup$
    – D. Hodge
    Jul 23, 2020 at 18:32
  • $\begingroup$ What about numerical errors? 1 m compared to the Earth radius would require at least 23 bits, that is the resolution of single precision floating point. Do you use single or double precision floats? If time increments are too small, double would not do. $\endgroup$
    – Uwe
    Jul 23, 2020 at 18:54

1 Answer 1


As noted in the comments, both methods potentially break down. Kepler's method is mathematically exact, but as formulated here it becomes ill-conditioned as the eccentricity approaches 1. The mean and eccentric anomalies become undefined for a parabolic orbit because a parabola lacks a center. Numerov's/Cowell's method approximates the propagator with a Taylor series whose accuracy is not guaranteed for eccentricities all the way up to a limiting value of 1.

To separate these issues we should put Kepler's method into a form that remains well-conditioned all the way to the parabolic limit. This proper conditioning will then guarantee the accuracy of the Kepler result, against which the Numerov/Cowell result can be compared.

In this discussion the time from periapsis is used together with the anomaly conversions given below, the latter taken from Wikipedia:

$t=\sqrt{\dfrac{\alpha^3}{\mu}}M=\sqrt{\dfrac{r_p^3}{(1-e)^3\mu}}M \tag{1}\label{Eq 1}$

$M=E-e\sin E\tag{2}\label{Eq 2}$

$E=2\tan^{-1}\left(\sqrt{\dfrac{1-e}{1+e}}\tan\dfrac{\theta}{2}\right)\tag{3}\label{Eq 3}$

In Eq. 1 $r_p$ is the periapsis distance, which unlike the semi major axis remains bounded and well-defined all the way up to (and beyond) eccentricity $1$.

When $E$ is calculated from $\theta$ via Eq. 3 with $e$ approaching 1, we get $E$ proportional to $(1-e)^{1/2}$; but in order to have finite and bounded times the mean anomaly $M$ must be proportional to $(1-e)^{3/2}$. Therefore Eq. 2 which connects $M$ to $E$ is ill-conditioned, for we are inputting terms with the lower power proportionality to get a difference with the higher-power proportionality.

To get $M$ in terms of quantities having the proper proportionality for a well-conditioned operation, render

$M=E-e\sin E = (E-\sin E)+(1-e)\sin E = ((\sin^{-1}s)-s)+(1-e)s$

where $s=\sin E$. This equation holds for $|E|\le\pi/2$ corresponding to $\theta\le 2\tan^{-1}(\sqrt{(1+e)/(1-e))}$. For larger true anomalies we can use Eqs. 1-3 directly because the ill-conditioning does not arise. Hereafter we focus on the case $\theta\le 2\tan^{-1}(\sqrt{(1+e)/(1-e))}$.

From Eq. 3 and the trigonometric identity $\sin(2\tan^{-1}u)=2u/(1+u^2)$ we obtain

$s=\dfrac{2\sqrt{\dfrac{1-e}{1+e}}\tan\dfrac{\theta}{2}}{1+\dfrac{1-e}{1+e}\tan^2\dfrac{\theta}{2}}\tag{4}\label{Eq 4}$

Now we have to tackle $\sin^{-1}s-s$ which is proportional to $s^3$ whereas the terms are proportional to $s$. To eliminate this ill-conditioning, convert this transcendental function to an integral involving an algebraic one:

$\displaystyle{(\sin^{-1}s)-s=\int_0^s\left(\dfrac{1}{\sqrt{1-x^2}}-1\right) dx = \int_0^s\left(\dfrac{x^2}{\sqrt{1-x^2}(\sqrt{1-x^2}+1)}\right) dx\tag{5}\label{Eq 5}}$

We have invoked the difference of squares factorization to rationalize the numerator, which gets rid on the ill-conditioned subtraction. This may be integrated numerically to obtain a well-conditioned result at high eccentricity. Since the integration is independent of the orbit, values as a function of $s$ may be stored beforehand in a table for lookup.

Put this all together and we obtain the time from periapsis for a highly eccentric orbit from a modified version of Eq. 1:

$\displaystyle{t=\sqrt{\dfrac{r_p^3}{\mu}}\left(\dfrac{1}{(1-e)^{3/2}}\int_0^s\left(\dfrac{x^2}{\sqrt{1-x^2}(\sqrt{1-x^2}+1)}\right) dx+\dfrac{s}{(1-e)^{1/2}}\right)}\tag{6}\label{Eq 6}$

with $s$ determined from Eq. 4.

  • 2
    $\begingroup$ Thanks @uhoh. It took several attempts but I finally got it to work! $\endgroup$ Jul 26, 2020 at 16:28
  • 2
    $\begingroup$ "bonus points" for "continuous improvement" :-) $\endgroup$
    – uhoh
    Aug 31, 2021 at 23:03
  • $\begingroup$ To avoid confusion, it should be said that $p$ is usually the parameter (or semi latus rectum as the mathematician say) of the orbit. Periapsis distance is a(1-e) so p in eq 3 is periapsis distance. Parameter is a(1-e*2). $\endgroup$
    – W H G
    Sep 1, 2021 at 15:40
  • $\begingroup$ @whg Using the answer to another question as a guide, I changed the variable name to read $r_p$. $\endgroup$ Sep 1, 2021 at 18:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.