# Does Kepler's Method for Determining Time of flight between two True Anomalies break down with eccentricities approaching 1?

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

• 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. – Oscar Lanzi Jul 23 at 1:08
• 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. – Uwe Jul 23 at 13:19
• @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. – D. Hodge Jul 23 at 18:32
• 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. – Uwe Jul 23 at 18:54

As noted in the comments, both methods potentially break down. Kepler's method is mathematically exact, but as formulated here it breaks down 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 us 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{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 $$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:

$$t=\sqrt{\dfrac{p^3}{\mu}}\left(\dfrac{(\sin^{-1}s)-s}{(1-e)^{3/2}}+\dfrac{s}{(1-e)^{1/2}}\right)\tag{6}\label{Eq 6}$$

with $$s$$ determined from Eq. 4 and $$(\sin^{-1}s)-s$$ determined from Eq. 5.

• I started to make a MathJax edit but then realized it's kind-of tricky with the integrals. Parenthesis that grow with the size of the fraction can be gotten as follows; $$\left(\frac{1}{1+\frac{1}{2}}\right$$ renders as $$\left(\frac{1}{1+\frac{1}{2}}\right).$$ If you scroll down on the MathJax help page to Contents you can see a link to the section on [equation numbering]() as well. – uhoh Jul 26 at 15:53
• Thanks @uhoh. It took several attempts but I finally got it to work! – Oscar Lanzi Jul 26 at 16:28
• Beautiful! btw I see that if you insert \displaystyle in the beginning of an equation that contains integral signs (right after the dollar sign), they grow to full size like this: $\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}$ – uhoh Jul 26 at 17:51