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

Question

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

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.