\( \DeclareMathOperator{\sgn}{sgn} \DeclareMathOperator{\d}{d} \DeclareMathOperator{\artanh}{artanh} \DeclareMathOperator{\arsinh}{arsinh} \newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}} \)
Kepler's Equation is an equation that describes the relationship between time and place for objects in an elliptic orbit under the influence of gravity. There are similar equations for parabolic and hyperbolic orbits, which we for convenience here also call Kepler's Equations.
Using Kepler's Equation, you can calculate the time at which a certain position is reached. If you want to calculate the position at a certain time, then for elliptic and hyperbolic orbits this cannot be done as a calculation involving a small number of evaluations of well-known functions, but only by first estimating a solution and then successively approaching the true solution ever closer until the solution is close enough to the true one.
Using the methods described below, you can calculate \( τ = \tan\left( \frac{1}{2}ν \right) \), where \( ν \) is the true anomaly, which describes the position of the planet (or other celestial object) in its orbit around the Sun (or other celestial object). The true anomaly is the angle between the direction to the perifocus and the direction to the planet, as seen from the focus of the orbit.
For an elliptic orbit, the eccentricity \( e \) is greater or equal to 0 and less than 1. The mean anomaly \( M \), eccentric anomaly \( E \), and true anomaly \( ν \) (Greek letter "nu") of an object in an elliptic orbit are determined by
\begin{align} M & = t\sqrt{\frac{Γ}{a^3}} \\ M & = E - e \sin(E) \label{eq:ekepler} \\ τ = \tan\left( \frac{1}{2}ν \right) & = \sqrt{\frac{1 + e}{1 - e}} \tan\left( \frac{1}{2}E \right) \label{eq:eEtoν} \end{align}
where \( a \) is the length of the semimajor axis of the orbit, \( Γ \) is the gravity parameter of the system, \( t \) is the time since the (last) passage through the perifocus, \( e \) is the eccentricity of the orbit, and all angles are measured in radians (2π radians equals 360 degrees). \( M \) increases at a constant rate (at 2π radians per orbital period), so it is a measure of time. \( ν \) is the angular distance of the object from its perifocus, as seen from the focus.
Equation \eqref{eq:ekepler} is Kepler's Equation (for elliptic orbits). The problem is that we want to know \( E \) given \( M \) and \( e \), and that equation cannot be rearranged into a form \( E = … \) in which the right-hand side is a finite collection of terms that do not have \( E \) in them.
The eccentricity \( e \) of a parabolic orbit is equal to exactly 1. For such an orbit, there is a direct method to find the position as a function of time, without approximations:
\begin{align} W & = t\sqrt{\frac{9Γ}{8q^3}} \label{eq:pkepler} \\ u & = \left( W + \sqrt{W^2 + 1} \right)^{1/3} \\ \tan\left( \frac{1}{2}ν \right) & = u - \frac{1}{u} \end{align}
For a hyperbolic orbit, the eccentricity \( e \) is greater than 1.
\begin{align} M & = t\sqrt{\frac{Γ}{|a|^3}} \\ M & = e \sinh E - E \label{eq:hkepler} \\ \tan\left( \frac{1}{2}ν \right) & = \sqrt{\frac{e + 1}{e - 1}} \tanh\left( \frac{1}{2}E \right) \label{eq:hEtoν} \end{align}
Equation \eqref{eq:hkepler} has the same problem that equation \eqref{eq:ekepler} has.
To find the solution, you can use the following methods, depending on the type of orbit.
If the eccentricity \( e \) is less than 1, then the orbit is a circle (\( e = 0 \)) or ellipse. The method is then as follows, for mean anomaly \( M \) measured in radians:
Replace \( M \) if possible by an equivalent smaller value that is definitely between −2π and +2π, and preferably between −π and +π, by adding or subtracting the appropriate multiple of 2π. How you can do this depends a lot on the possibilities of your calculating machine or programming language.
\begin{equation} M ← M \bmod 2π \end{equation}
Calculate the initial estimate \( E_1 \) of the solution.
\begin{align} s_1 & = \left| \frac{M}{δ} \right| \\ s_2 & = (6|M|)^{1/3} \\ E_1 & = \sgn(M) \min(s_1, s_2) \end{align}
If you want it to be less complicated, then you can take always
\begin{equation} E_1 = (6M)^{1/3} \end{equation}
Calculate a correction from the previous estimate:
\begin{equation} d_i = E_i + \frac{M + e \sin(E_i) - E_i}{1 - e \cos(E_i)} \end{equation}
If \( d_i \) is exactly equal to 0, then you have found \( E \). Stop the procedure.
If \( |d_i| \) is greater than \( |d_{i-1}| \) and also is less than \( |ε E_i| \), then you cannot find a better value for \( E \). Stop the procedure.
If \( |d_i| \) is greater than \( |d_{i-1}| \) and \( i \) is greater than \( n \), then give up. The current method (and especially the initial estimate) is not suitable to find the value of \( E \) for the current \( M \) and \( e \). (Because of the well-chosen initial estimate, you should never end up here, but one never knows.)
If you haven't quit yet, then return to step 3.
The value of \( ε \) must fit the relative precision with which you can calculate. If your calculator or programming language can just tell the difference between 1 and \( 1 + ε_0 \) (so if you add \( ε_0 \) to 1, and then subtract 1 from the result, then you find a result that is not 0, but if you do the same with a value smaller than \( ε_0 \) then you do end up with 0), then \( ε = \sqrt{ε_0} \) seems a reasonable value to use. I used 10^{−7}, in a computer program that calculated using double precision floating-point values.
A reasonable value for \( n \) is 10.
The last-found \( E_i \) is your best estimate for \( E \). Use Equation \eqref{eq:eEtoν} to calculate \( \tan\left( \frac{1}{2}ν \right) \).
If the eccentricity \( e \) is equal to 1, then the orbit is a parabola. Then use Equations \eqref{eq:pkepler}ff to calculate \( \tan\left( \frac{1}{2}ν \right) \).
If the eccentricity \( e \) is greater than 1, then the orbit is a hyperbola. The method for hyperbolic orbits is in many but not all details the same as that for elliptic orbits.
Do not replace \( M \) by a smaller value (which did get asked for elliptic orbits).
Calculate the initial estimate \( E_1 \) for the solution. If \( |M| \lt 3e \), then
\begin{align} s_1 & = \left| \frac{M}{δ} \right| \\ s_2 & = (6|M|)^{1/3} \\ E_1 & = \sgn(M) \min(s_1, s_2) \end{align}
Otherwise,
\begin{equation} E_1 = \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) \end{equation}
If you want it to be less complicated, then you can take always
\begin{equation} E_1 = (6M)^{1/3} \end{equation}
Calculate a correction from the previous estimate:
\begin{align} N_i & = \frac{M + E_i}{e \cosh(E_i)} - \tanh(E_i) \\ D_i & = 1 - \frac{1}{e \cosh(E_i)} \\ d_i & = \frac{N_i}{D_i} \end{align}
The termination conditions are the same as those for elliptic orbits. If \( d_i \) is exactly equal to 0, then you have found \( E \). Stop the procedure.
If \( |d_i| \) is greater than \( |d_{i-1}| \) and also is less than \( |ε E_i| \), then you cannot find a better value for \( E \). Stop the procedure.
If \( |d_i| \) is greater than \( |d_{i-1}| \) and \( i \) is greater than \( n \), then give up. The current method (and especially the initial estimate) is not suitable to find the value of \( E \) for the current \( M \) and \( e \). (Because of the well-chosen initial estimate, you should never end up here, but one never knows.)
If you haven't quit yet, then return to step 3.
The values of \( ε \) and \( n \) are the same as for elliptic orbits.
The last-found \( E_i \) is your best estimate for \( E \). Use Equation \eqref{eq:hEtoν} to calculate \( \tan\left( \frac{1}{2}ν \right) \).
Given \( τ = \tan\left( \frac{1}{2}ν \right) \), the eccentricity \( e \), and the perifocal distance \( q = a |1 - e| \) between the focus and the perifocus, the distance \( r \) between the focus and the target object is, for all circular, elliptic, parabolic, and hyperbolic orbits:
\begin{align} ρ & = \frac{1 + e}{1 + e + (1 - e)τ^2} \\ r & = qρ (1 + τ^2) \notag \\ & = q \frac{(1 + e) (1 + τ^2)}{1 + e + (1 - e) τ^2} \\ & = q \frac{1 + e}{1 + e \cos ν} \end{align}
For elliptic orbits:
\begin{equation} r = q \frac{1 - e \cos E}{1 - e} \end{equation}
For parabolic orbits:
\begin{align} r & = q \frac{2}{1 + \cos ν} \notag \\ & = q (1 + τ^2) \end{align}
For hyperbolic orbits:
\begin{equation} r = q\frac{e \cosh(E) - 1}{e - 1} \end{equation}
If you know \( τ = \tan\left( \frac{1}{2}ν \right) \), then you can calculate \( \sin(ν) \) and \( \cos(ν) \):
\begin{align} \sin(ν) & = \frac{2τ}{1 + τ^2} \\ \cos(ν) & = \frac{1 - τ^2}{1 + τ^2} \end{align}
and then the coordinates of the target object (in the orbital plane):
\begin{align} x & = r \cos(ν) = qρ (1 - τ^2) = q \frac{(1 + e) (1 - τ^2)}{1 + e + (1 - e) τ^2} \label{eq:coords} \\ y & = r \sin(ν) = 2qρτ = q \frac{2 (1 + e) τ}{1 + e + (1 - e) τ^2} \end{align}
For a parabolic orbit, \( e = 1 \) so \( δ = 0 \). For a nearly-parabolic orbit (\( |δ| \ll 1 \)) we have
\begin{align} ρ & ≈ 1 + \frac{1}{2}δτ^2 \\ r & ≈ q\left( 1 + τ^2 + \frac{1}{2}δτ^2 (1 + τ^2) \right) \\ x & ≈ q\left( 1 - τ^2 + \frac{1}{2}δτ^2 (1 - τ^2) \right) \\ y & ≈ 2qτ \left( 1 + \frac{1}{2}δτ^2 \right) \end{align}
I define the perifocal anomaly \( m \) as follows:
\begin{equation} m = \frac{M}{|δ|^{3/2}} \end{equation}
The maximum perifocal anomaly for an elliptic orbit is \( μ/|δ|^{3/2} \).
For parabolic orbits,
\begin{align} W & = \frac{3m}{2^{3/2}} \\ u & = (W + \sqrt{W^2 + 1})^{1/3} \\ τ & = u - \frac{1}{u} \\ m \ll 1 & ⇒ τ ≈ \frac{m}{\sqrt{2}} \\ m \gg 1 & ⇒ τ ≈ \frac{(6m)^{1/3}}{\sqrt{2}} \end{align}
For circular, elliptic, and hyperbolic orbits, with \( q = a|1 - e| = a|δ| \), we have
\begin{align} M & = t \sqrt{\frac{Γ}{a^3}} = t \sqrt{\frac{Γ|δ|^3}{q^3}} \\ m & = t \sqrt{\frac{Γ}{q^3}} \end{align}
so the perifocal anomaly \( m \) for nearly-parabolic orbits with fixed perifocal distance \( q \) is independent of how close the shape of the orbit is to a parabola (i.e., of how small \( δ \) is), while the mean anomaly \( M \) depends very strongly on how small \( δ \) is.
For a nearly-parabolic orbit (with \( |δ| \ll 1 \) we have:
\begin{align} |m| \ll 1 & ⇒ τ ≈ \frac{m}{\sqrt{2}} \\ 1 \ll |m| \ll \frac{1}{|δ|^{3/2}} & ⇒ τ ≈ \frac{(6m)^{1/3}}{\sqrt{2}} \\ |m| \gg \frac{1}{|δ|^{3/2}} & ⇒ τ ≈ \sgn(m)\sqrt{\frac{2}{|δ|}} \left( 1 - \frac{1}{|m||δ|^{3/2}} \right) \end{align}
so for \( |m| \ll |δ|^{−3/2} \) an object in such an orbit behaves itself as if it is in a parabolic orbit.
The next few tables show a some specific solutions to Kepler's Equation, which we'll use later to check the performance of approximation formulas. All angles are given in radians. One radian is 180/π ≈ 57,2957795131 degrees, and
\[ τ = \tan\left( \frac{1}{2} ν \right) \]
For cases where \( e \) is close to 1, the following definitions are useful:
\begin{align} δ & ≡ e - 1 \\ m & ≡ \frac{M}{|δ|^{3/2}} \label{eq:m} \end{align}
\( m \) is the perifocal anomaly.
Table 1 is for an anomaly of 0.0001 radians.
Table 1: Equation of Kepler: Solutions for 0.0001 Radian
# | \(M\) | \(m\) | \(e\) | \(E\) | \(τ\) | \(ν\) |
---|---|---|---|---|---|---|
1 | 0.0001 | 0.000100000000 | 0 | 0.000100000000 | 5.00000000 × 10^{−5} | 0.000100000000 |
2 | 0.0001 | 0.000101518971 | 0.01 | 0.000101010101 | 5.10126517 × 10^{−5} | 0.000102025303 |
3 | 0.0001 | 0.00316227766 | 0.9 | 0.000999998500 | 0.00217944638 | 0.00435888587 |
4 | 0.0001 | 0.100000000 | 0.99 | 0.00998358122 | 0.0704184571 | 0.140604812 |
5 | 0.0001 | 3.16227766 | 0.999 | 0.0614230944 | 1.37355061 | 1.88299657 |
6 | 0.0001 | 100.000000 | 0.9999 | 0.0819842185 | 5.80026395 | 2.80013747 |
7 | 0.0001 | 100.000000 | 1.0001 | 0.0819610818 | 5.79242631 | 2.79968440 |
8 | 0.0001 | 3.16227766 | 1.001 | 0.0613913007 | 1.37266327 | 1.88238152 |
9 | 0.0001 | 0.100000000 | 1.01 | 0.00998325102 | 0.0707679178 | 0.141300268 |
10 | 0.0001 | 0.00316227766 | 1.1 | 0.000999998167 | 0.00229128346 | 0.00458255889 |
11 | 0.0001 | 1.01518971 × 10^{−7} | 100 | 1.01010101 × 10^{−6} | 5.10126517 × 10^{−7} | 1.02025303 × 10^{−6} |
12 | 0.0001 | 1.00000150 × 10^{−13} | 1000000 | 1.00000100 × 10^{−10} | 5.00001000 × 10^{−11} | 1.00000200 × 10^{−10} |
13 | 9.85037563 × 10^{−5} | 0.0001 | 0.01 | 9.94987437 × 10^{−5} | 5.02493781 × 10^{−5} | 0.000100498756 |
14 | 3.16227766 × 10^{−6} | 0.0001 | 0.9 | 3.16227766 × 10^{−5} | 6.89202437 × 10^{−5} | 0.000137840487 |
15 | 1.00000000 × 10^{−7} | 0.0001 | 0.99 | 9.99999998 × 10^{−6} | 7.05336798 × 10^{−5} | 0.000141067359 |
16 | 3.16227766 × 10^{−9} | 0.0001 | 0.999 | 3.16227765 × 10^{−6} | 7.06929981 × 10^{−5} | 0.000141385996 |
17 | 1.00000000 × 10^{−10} | 0.0001 | 0.9999 | 9.99999998 × 10^{−7} | 7.07089102 × 10^{−5} | 0.000141417820 |
18 | 0 | 0.0001 | 1 | 0 | 7.07106780 × 10^{−5} | 0.000141421356 |
19 | 1.00000000 × 10^{−10} | 0.0001 | 1.0001 | 9.99999998 × 10^{−7} | 7.07124457 × 10^{−5} | 0.000141424891 |
20 | 3.16227766 × 10^{−9} | 0.0001 | 1.001 | 3.16227765 × 10^{−6} | 7.07283535 × 10^{−5} | 0.000141456707 |
21 | 1.00000000 × 10^{−7} | 0.0001 | 1.01 | 9.99999998 × 10^{−6} | 7.08872343 × 10^{−5} | 0.000141774468 |
22 | 3.16227766 × 10^{−6} | 0.0001 | 1.1 | 3.16227765 × 10^{−5} | 7.24568836 × 10^{−5} | 0.000144913767 |
23 | 0.0985037563 | 0.0001 | 100 | 0.000994987271 | 0.000502493656 | 0.00100498723 |
24 | 99999.8500 | 0.0001 | 1000000 | 0.0998340290 | 0.0498756461 | 0.0996687023 |
Line 4 in Table 1 says that for \( M = 0.0001 \) radians and \( e = 0.99 \), the eccentric anomaly is equal to \( E = 0.00998358122 \) radians, and the true anomaly is equal to \( ν = 0.140604809 \) radians = 8.05606213°. Is that correct?
Equation \eqref{eq:ekepler} says that
\[ M = E - e \sin(E) \]
With the current values of \( E \) and \( e \) we find
\[ E - e \sin(E) = 0.00998358122 - 0.99 × \sin(0.00998358122) = 0.0001 \]
which is equal to \( M \). Equation \eqref{eq:eEtoν} then yields
\begin{align*} \tan\left( \frac{1}{2}ν \right) & = \sqrt{\frac{1 + e}{1 - e}} \tan\left( \frac{1}{2}E \right) \\ & = \sqrt{\frac{1.99}{0.01}} × \tan(0.5×0.00998358122) \\ & = 14.106736 × \tan(0.00499179061) \\ & = 14.106736 × 0.00499183207 \\ & = 0.0704184571 \end{align*}
which fits what is listed in the τ column.
Then
\[ ν = 2\arctan(0.0704184571) = 0.140604809 \]
which also fits what it says in the table.
Table 2 is for an anomaly of 1 radian.
Table 2: Equation of Kepler: Solutions for 1 Radian
# | \(M\) | \(m\) | \(e\) | \(E\) | \(τ\) | \(ν\) |
---|---|---|---|---|---|---|
1 | 1 | 1.00000000 | 0 | 1.00000000 | 0.546302490 | 1.00000000 |
2 | 1 | 1.01518971 | 0.01 | 1.00846012 | 0.557353696 | 1.01694301 |
3 | 1 | 31.6227766 | 0.9 | 1.86208669 | 5.85747591 | 2.80340907 |
4 | 1 | 1000.00000 | 0.99 | 1.92763555 | 20.3140949 | 3.04321826 |
5 | 1 | 31622.7766 | 0.999 | 1.93387356 | 64.8144720 | 3.11073780 |
6 | 1 | 1000000.00 | 0.9999 | 1.93449428 | 205.143679 | 3.13184347 |
7 | 1 | 1000000.00 | 1.0001 | 1.72897376 | 98.7940852 | 3.12134922 |
8 | 1 | 31622.7766 | 1.001 | 1.72768618 | 31.2337093 | 3.07758114 |
9 | 1 | 1000.00000 | 1.01 | 1.71487376 | 9.85240023 | 2.93928924 |
10 | 1 | 31.6227766 | 1.1 | 1.59281168 | 3.03376885 | 2.50477756 |
11 | 1 | 0.00101518971 | 100 | 0.0101008366 | 0.00510113418 | 0.0102021799 |
12 | 1 | 1.00000150 × 10^{−9} | 1000000 | 1.00000100 × 10^{−6} | 5.00001000 × 10^{−7} | 1.00000200 × 10^{−6} |
13 | 0.985037563 | 1 | 0.01 | 0.993416520 | 0.547483734 | 1.00181857 |
14 | 0.0316227766 | 1 | 0.9 | 0.282532839 | 0.619895127 | 1.10983994 |
15 | 0.00100000000 | 1 | 0.99 | 0.0885485963 | 0.624974249 | 1.11716160 |
16 | 3.16227766 × 10^{−5} | 1 | 0.999 | 0.0279769359 | 0.625467687 | 1.11787112 |
17 | 1.00000000 × 10^{−6} | 1 | 0.9999 | 0.00884630818 | 0.625516891 | 1.11794185 |
18 | 0 | 1 | 1 | 0 | 0.625522357 | 1.11794971 |
19 | 1.00000000 × 10^{−6} | 1 | 1.0001 | 0.00884613583 | 0.625527822 | 1.11795757 |
20 | 3.16227766 × 10^{−5} | 1 | 1.001 | 0.0279714858 | 0.625576995 | 1.11802825 |
21 | 0.00100000000 | 1 | 1.01 | 0.0883762467 | 0.626067340 | 1.11873295 |
22 | 0.0316227766 | 1 | 1.1 | 0.277078928 | 0.630836813 | 1.12557114 |
23 | 985.037563 | 1 | 100 | 2.98623497 | 0.912981379 | 1.47988203 |
24 | 999998500. | 1 | 1000000 | 7.60090122 | 0.999001498 | 1.56979733 |
Line 18 from Table 2 says that for \( m = 1 \) and \( e = 1 \) (i.e., for a parabolic orbit) we have dat \( \tan\left( \frac{1}{2}ν \right) = 0.625522357 \). Is that correct?
According to Equations \eqref{eq:pkepler}ff and \eqref{eq:m} we have
\[ W = \sqrt{\frac98}m = 1.06066017×1 = 1.06066017 \]
and then
\begin{align*} u & = \left( W + \sqrt{W^2 + 1} \right)^{1/3} \\ & = \left( 1.06066017 + \sqrt{1.06066017^2 + 1} \right)^{1/3} \\ & = (1.06066017 + 1.45773797)^{1/3} \\ & = 2.51839815^{1/3} \\ & = 1.36053002 \end{align*}
and then
\begin{align*} τ = \tan\left( \frac{1}{2}ν \right) & = u - \frac{1}{u} \\ & = 1.36053002 - \frac{1}{1.36053002} \\ & = 0.625522357 \end{align*}
which fits what it says in the τ column.
Table 3 is for an anomaly of 10000 radians. That is relevant only for parabolic and hyperbolic orbits.
Table 3: Equation of Kepler: Solutions for 10000 Radians
# | \(M\) | \(m\) | \(e\) | \(E\) | \(τ\) | \(ν\) |
---|---|---|---|---|---|---|
1 | 10000 | 1.00000000 × 10^{+10} | 1.0001 | 9.90437751 | 141.410763 | 3.12744969 |
2 | 10000 | 316227766. | 1.001 | 9.90347791 | 44.7280654 | 3.09688545 |
3 | 10000 | 10000000.0 | 1.01 | 9.89452619 | 14.1760164 | 3.00074262 |
4 | 10000 | 316227.766 | 1.1 | 9.80915781 | 4.58207213 | 2.71184720 |
5 | 10000 | 10.1518971 | 100 | 5.29887209 | 1.00000580 | 1.57080212 |
6 | 10000 | 1.00000150 × 10^{−5} | 1000000 | 0.00999984334 | 0.00499988501 | 0.00999968669 |
7 | 0 | 10000 | 1 | 0 | 27.6461704 | 3.06928143 |
8 | 0.0100000000 | 10000 | 1.0001 | 0.389974639 | 27.2318138 | 3.06818213 |
9 | 0.316227766 | 10000 | 1.001 | 1.20643179 | 24.1257778 | 3.05874120 |
10 | 10.0000000 | 10000 | 1.01 | 3.27015981 | 13.1393971 | 2.98967154 |
11 | 316.227766 | 10000 | 1.1 | 6.37425935 | 4.56697679 | 2.71047028 |
12 | 9850375.63 | 10000 | 100 | 12.1909984 | 1.01004025 | 1.58078634 |
13 | 9.99998500 × 10^{+12} | 10000 | 1000000 | 16.8112413 | 1.00000090 | 1.57079723 |
Line 3 of Table 3 says that for \( M = 10000 \) radians and \( e = 1.01 \), the eccentric anomaly is equal to 9.89452619, and the true anomaly is equal to \( ν = 2×\arctan(14.1760164) = 3.00074267 \) radians = 171.929891°. Is that correct?
According to Equation \eqref{eq:hkepler} we must then have
\[ M = e \sinh(E) - E \]
With the current values of \( E \) and \( e \) we find
\[ e \sinh(E) - E = 1.01 × \sinh(9.89452619) - 9.89452619 = 10000 \]
which is exactly the value of \( M \). Equation \eqref{eq:hEtoν} then leads to
\begin{align*} τ = \tan\left( \frac{1}{2}ν \right) & = \sqrt{\frac{e + 1}{e - 1}} \tanh\left( \frac{1}{2}E \right) \\ & = \sqrt{\frac{2.01}{0.01}} × \tanh(0.5 × 9.89452619) \\ & = 14.1774469 × \tanh(4.9472631) \\ & = 14.1774469 × 0.999899105 \\ & = 14.1760164 \end{align*}
which fits the value from the τ column.
The above tables show that near \( e = 1 \) the values of \( E \) and \( τ \) vary much faster with \( e \) for fixed \( M \) than for fixed \( m \). A parabolic orbit (with \( e = 1 \)) cannot be described at all with \( M \) or \( E \) (which are exactly 0 then), but they can be described with \( m \).
In extreme cases, when \( M \) or \( e \) or a related quantity is very small or rather very large, then we can find reasonable approximate solutions for \( E \).
If \( |M| ≪ \sqrt{\frac{|δ|^3}{e}} \) and also \( |M| ≪ |δ| \), so the mean anomaly is very small, then
\begin{equation} E ≈ \frac{M}{|δ|} \label{eq:E1} \end{equation}
For example, when \( M = 0.0001 \) and \( e = 0.01 \), then \( M ≪ \sqrt{\frac{|δ|^3}{e}} ≈ 0.001 \) and \( |M| ≪ δ ≈ 0.99 \), so we've met the requirements for using the approximation, so
\[ E ≈ \frac{M}{|δ|} = \frac{0.001}{0.99} ≈ 0.000101010101 \]
identical to what it says on line 2 of Table 1.
If \( M \) is very small and also \( e ≪ 1 \), so the orbit is nearly circular, then
\begin{align} τ & ≈ \frac{1}{2} M ( 1 + 2e ) \\ ν & ≈ (1 + 2e) M \end{align}
For example, if \( M = 0.0001 \) and \( e = 0.01 \), then we've met the requirements for using the approximation (see above), so
\[ τ ≈ \frac{1}{2} M (1 + 2e) = 0.5 × 0.0001 × 1.02 ≈ 0.00005100 \]
compared to 0.00005101 from the table, and
\[ ν ≈ (1 + 2e) M = 1.02 × 0.0001 = 0.000102 \]
compared to 0.00010203 from the table.
If \( M \) is very small and also \( |δ| ≪ 1 \), so the orbit is nearly parabolic, then
\begin{align} τ & ≈ \frac{M}{\sqrt{2|δ|^3}} = \frac{m}{\sqrt{2}} \\ ν & ≈ \frac{M\sqrt{2}}{\sqrt{|δ|^3}} = m\sqrt{2} \end{align}
For example, when \( M = 0.0001 \) and \( e = 0.99 \), then \( M ≪ \sqrt{\frac{|δ|^3}{e}} ≈ 10^{−6} \) and \( M ≪ |δ| ≈ 0.01 \), then we've met the requirements for using the approximation, so
\[ τ ≈ \frac{M}{\sqrt{2|δ|^3}} = \frac{0.0001}{\sqrt{2×0.01^3}} ≈ 0.0707 \]
compared to 0.0704 on line 4 in Table 1, and
\[ ν ≈ \frac{0.0001 × \sqrt{2}}{0.01^3} = 0.1414 \]
compared to 0.1406 from the table.
If \( |δ|^{3/2} ≪ M ≪ 1 \) and also \( |δ| ≪ 1 \), so the mean anomaly is mid-sized and the orbit is nearly parabolic, then
\begin{equation} E ≈ (6M)^{1/3} \label{eq:E2} \end{equation}
and
\begin{align} τ & ≈ \frac{(6M)^{1/3}}{\sqrt{2|δ|}} = \frac{(6m)^{1/3}}{\sqrt{2}} \\ ν & ≈ \sgn(M) \left( π - \frac{\sqrt{8|δ|}}{(6M)^{1/3}} \right) = \sgn(m) \left( π - \frac{\sqrt{8}}{(6m)^{1/3}} \right) \end{align}
For example, when \( M = 0.0001 \) and \( e = 0.9999 \), then \( 10^{−6} = |δ|^{3/2} ≪ M ≪ 1 \), so we've met the requirements for using the approximation, so
\[ E ≈ (6M)^{1/3} = 0.0006^{1/3} ≈ 0.0843 \]
and
\[ τ ≈ \frac{(6M)^{1/3}}{\sqrt{2|δ|}} ≈ \frac{(6×0.0001)^{1/3}}{\sqrt{2×0.0001}} ≈ 5.96 \]
compared to 5.80 on line 6 in Table 1. Then
\[ ν ≈ \sgn(0.0001) × \left( π - \frac{\sqrt{8×0.0001}}{(6×0.0001)^{1/3}} \right) ≈ π - 0.3353 ≈ 2.806 \]
compared with 2.800 from the table.
If \( M ≫ e \gt 1 \), so it is a hyperbolic orbit and the mean anomaly is very large, then
\begin{equation} E ≈ \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) \label{eq:E3} \end{equation}
For example, when \( M = 10000 \) and \( e = 100 \) then \( M ≫ e \), so we've met the requirements for using the approximation, so
\[ E ≈ \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) = \sgn(10000)×\log\left( 1 + \frac{2×10000}{100} \right) ≈ 5.303 \]
compared with 5.299 on line 5 in Table 3.
If also \( |δ| ≪ 1 \), so the orbit is nearly parabolic, then
\begin{align} τ & ≈ \sgn(M) \sqrt{\frac{2}{|δ|}} \left(1 - \frac{1}{|M|} \right) \\ ν & ≈ \sgn(M) \left( π - \frac{\sqrt{2|δ|}|M|}{|M| - 1} \right) \end{align}
For example, when \( M = 10000 \) and \( e = 1.01 \) then \( M ≫ e \) and \( |δ| = 0.01 ≪ 1 \), so we've met the requirements for using the approximation, so
\[ τ ≈ \sgn(M) \sqrt{\frac{2}{|δ|}} \left( 1 - \frac{1}{|M|} \right) = \sgn(10000)×\sqrt{\frac{2}{0.01}}×\left( 1 - \frac{1}{10000} \right) ≈ 14.14 \]
compared to 14.18 in line 3 in Table 3. And
\[ ν ≈ \sgn(10000) × \left( π - \frac{\sqrt{2×0.01} × 10000}{9999} \right) ≈ π - 0.1414 ≈ 3.0002 \]
compared to 3.0007 in the table.
If the eccentricity is very large and also \( e ≫ 1 \), so the orbit is strongly hyperbolic, then
\begin{align} τ & ≈ \sgn(M) \left( 1 + \frac{1}{e} - \frac{e}{|M|} \right) \\ ν_∞ & = \arctan\left( \sqrt{\frac{e + 1}{e - 1}} \right) = \arccos\left( -\frac{1}{e} \right) \\ ν & ≈ ν_∞ \sgn(M) - \frac{\sqrt{e^2 - 1}}{M} \end{align}
For example, when \( M = 10000 \) and \( e = 100 \) then \( |M| ≫ e \), so we've met the requirements for using the approximation, so
\begin{align*} τ & ≈ \sgn(M) \left( 1 + \frac{1}{e} - \frac{e + 1}{|M|} \right) \\ & = \sgn(10000)×\left( 1 + \frac{1}{100} - \frac{101}{10000} \right) \\ & ≈ 0.9999 \end{align*}
compared with 1.000006 on line 5 in Table 3. And
\begin{align*} ν_∞ & = \arccos\left( -\frac{1}{e} \right) = \arccos\left( -\frac{1}{100} \right) ≈ 1.580796 \\ ν & ≈ ν_∞ - \frac{\sqrt{100^2 - 1}}{10000} ≈ ν_∞ - 0.0099995 ≈ 1.570797 \end{align*}
compared to 1.570802 in the table.
With Kepler's Equation you can calculate the time (\( t \) or \( M \)) when you know the location in the orbit (\( ν \) or \( E \)) and the shape of the orbit (\( e \)), but not easily calculate the location when you know the time and shape ― except for a parabolic orbit. To calculate the location from the time and the shape you need two things:
You repeat the second step with each next approximation, until the approximate value no longer changes: then you have found the solution to the precision of your calculations.
How quickly you find the location depends on how accurate your initial estimate is, how much better each next approximation is, and how little effort it takes to calculate the initial estimate and the successive approximations.
Many different formulas have been invented for the initial estimate and the next approximations, usually based on the behavior of the solutions to Kepler's Equation.
Suppose that we can write Kepler's Equation like
\begin{equation} M = f(E) \label{eq:f} \end{equation}
and can rewrite it to
\begin{equation} E = g(M,E) \label{eq:g} \end{equation}
for a suitable function \( g \). Then we could construct an approximation formula
\begin{equation} E_{i+1} = g(M,E_i) \end{equation}
so
\begin{align} E_2 & = g(M,E_1) \\ E_3 & = g(M,E_2) \\ E_4 & = g(M,E_3) \end{align}
and so on. Do we then get ever closer to the solution?
The difference between estimate \( E_i \) and the true value \( E \) is
\begin{equation} d_i ≡ E_i - E \end{equation}
Then
\begin{equation} E_{i+1} = g(M,E_i) = g(M,E + d_i) ≈ g(M,E) + g'(M,E) d_i \end{equation}
where \( g' = \pdd{g}{E} \) is the derivative of \( g \) with respect to \( E \), and where the last near-equality holds when \( d_i \) is very small, which is when the approximation is already close to the true solution. Then
\begin{equation} d_{i+1} ≈ g'(M,E) d_i \end{equation}
so the approximations only improve when \( |g'(M,E)| \lt 1 \), and the solution is found quickly only if \( |g'(M,E)| ≪ 1 \).
We get further with the method of Newton-Raphson. Equation \eqref{eq:f} leads to
\begin{equation} \d M = f'(E) \d E \end{equation}
where \( f' \) is the derivative of \( f \) with respect to \( E \). So, for small deviations we have
\begin{equation} ∆M ≈ f'(E) ∆E \end{equation}
We know
\begin{equation} ∆M = f(E_i) - M \end{equation}
so
\begin{equation} ∆E ≈ \frac{∆M}{f'(E_i)} ≈ \frac{f(E_i) - M}{f'(E_i)} \end{equation}
so the next approximation is
\begin{equation} E_{i+1} = E_i + \frac{M - f(E_i)}{f'(E_i)} \end{equation}
How much better is that one than the previous one? We have
\begin{align} E + d_{i+1} & = E + d_i + \frac{M - f(E + d_i)}{f'(E + d_i)} \notag \\ & ≈ E + \frac{1}{2} d_i^2 \frac{f''(E)}{f'(E)} \end{align}
where we've used \( M = f(E) \) and have assumed that \( d_i f'''(E) ≪ f''(E) \) and \( d_i f''(E) ≪ f'(E) \)
so
\begin{equation} d_{i+1} ≈ \frac{1}{2} d_i^2 \frac{f''(E)}{f'(E)} \end{equation}
so the last estimate is better than the previous one if \( |d_{i+1}| \lt |d_i| \) which means \( |d_i f''(E)| \lt 2|f'(E)| \) and then \( |d_i| \lt β ≡ 2\left| \frac{f'(E)}{f''(E)} \right|| \). The boundary indicated by \( β \) is vague (because found after neglecting certain terms), and at the beginning we don't know the value of \( E \) anyway and have to calculate \( β \) based on \( E_1 \) instead, which is not accurate, either. So if \( |d_i| \) isn't comfortably less than \( β \) then that may still not be small enough to locate the solution quickly. But if \( |d_i| ≪ β \) then the number of correct digits after the decimal point roughly doubles for each next estimate.
This method takes more calculation effort for each repeat than the previous method did, but take far fewer repeats to find the true solution than the previous method did.
A good initial estimate is important, but good stopping conditions are important, too. How do you know that you have found the true solution, or at least cannot get any closer to it anymore?
The most desirable outcome is of course that the correction that you calculate (based on how well the previous estimate did) is exactly zero, because then you cannot make the estimate more precise anymore.
But sometimes you do not find a correction of exactly zero, but end up in a cycle of successive small corrections, for example first a small positive correction, and then a small negative correction, and then a small positive one again, without ever getting to zero. To break out of such a cycle, one can adopt the rule that the iteration should stop if the new correction has a greater magnitude than the previous correction.
But that doesn't always work. If the initial estimate is far from the true solution, then the first few corrections can be quite large, and may sometimes be greater in magnitude than the previous correction, until you get close enough to the true solution, and then the corrections quickly decrease in magnitude. If you follow the rule "stop when the correction is bigger than before" right from the start, then you may stop too soon, when the estimate is still far from the true solution.
All in all, I use the following rules:
We need reasonable values for \( ε \) and \( n \). I've chosen \( ε = 10^{−7} \) and \( n = 10 \).
For elliptic orbits
\[ M = E - e \sin E \]
so for the first-order approximation
\begin{align} g(M,E) & = M + e \sin E \\ g'(M,E) & = e \cos E \end{align}
so the newest estimate is better than the previous one if \( |e \cos E| \lt 1 \), which is always true for elliptical orbits (with \( e \lt 1 \)). However, if \( e ≈ 1 \) and \( E ≈ 0 \) then the estimate improves only marginally in each step, so then it takes very many steps before you approach the solution. A very rough estimate of that number of steps is
\[ N ≈ \frac{1}{1 - |g'(M,E)|} = \frac{1}{1 - |e \cos E|} \]
which for \( e = 1 + δ ≈ 1 \) and \( |E| ≪ 1 \) is approximately equal to \( \frac{1}{|δ|} \), which for near-parabolic orbits is a very large number.
For the second-order approximation, we have
\begin{align} f(E) & = E - e \sin E \\ f'(E) & = 1 - e \cos E \\ f''(E) & = e \sin E \end{align}
and
\begin{equation} E_{i+1} = E_i + \frac{M + e \sin(E_i) - E_i}{1 - e \cos(E_i)} \end{equation}
\begin{equation} β = \left| \frac{2(1 - e \cos E)}{e \sin E} \right| \end{equation}
\( E_{i+1} \) is more accurate than \( E_i \) if \( |d_i| \lt β \). The least value of \( β \) is attained when \( \cos(E) = e \); then
\[ β = β_{\text{min}} = 2\frac{\sqrt{1 - e^2}}{e} \]
For near-parabolic orbits with \( e = 1 + δ \) and \( |δ| ≪ 1 \) we find \( β_{\text{min}} ≈ \sqrt{2|δ|} \) for \( |E| ≈ \sqrt{2|δ|} \), so then the initial estimate must be quite accurate already, to be able to locate the solution quickly.
Let \( M = 1 \), \( e = 0.5 \), and \( E_1 = 1.5 \). Then
\begin{align*} β & ≈ \left| 2×\frac{1 - 0.5×\cos(1.5)}{0.5×\sin(1.5)} \right| = 3.87 \\ E_2 & = E_1 + \frac{M + e \sin(E_1) - E_1}{1 - e \cos(E_1)} = 1.5 - 0.0012984304 = 1.4987015696 \\ E_3 & = E_2 + \frac{M + e \sin(E_2) - E_2}{1 - e \cos(E_2)} = 1.4987015696 - 0.0000004361 = 1.4987011335 \\ E_4 & = E_3 + \frac{M + e \sin(E_3) - E_3}{1 - e \cos(E_3)} = 1.4987011335 - 0.0000000000 = 1.4987011335 \end{align*}
So with this initial estimate, the true solution is already found after 3 repeats, to 10 decimals after the decimal point. The "capture area" reaches out to about 3.87 radians from the solution, so pretty much any initial estimate quickly leads to the solution.
Now a different example. Let \( M = 0.0001 \) and \( e = 0.999999 \) and \( E_1 = 0.0002 \). Then
\begin{align*} β & ≈ 2×\frac{1 - 0.999999×\cos(0.0002)}{0.999999×\sin(0.0002)} = 0.0102 \\ E_i & = E_{i-1} + \frac{N_i}{D_i} = E_{i-1} + ∆_i \end{align*}
\(i\) | \(N_i\) | \(D_i\) | \(∆_i\) | \(E_i\) |
---|---|---|---|---|
0 | 0.0002 | |||
1 | 0.0000999998 | 0.0000102 | 98.03902 | 98.03922 |
2 | −98.64418 | 1.796175 | −54.91903 | 3.305431 |
3 | −43.87931 | 0.3491598 | −125.6711 | −82.55095 |
4 | 81.78702 | 0.3548208 | 230.5023 | 147.9513 |
5 | −148.2434 | 1.956366 | −75.77488 | 72.17647 |
6 | −72.0963 | 1.996788 | −36.10614 | 36.07033 |
7 | −37.06855 | 1.057949 | −35.03813 | 1.032199 |
8 | −0.173671 | 0.4870684 | −0.3565639 | 0.6756354 |
9 | −0.0514278 | 0.2196911 | −0.2282422 | 0.4473932 |
10 | −0.01467686 | 0.09842299 | −0.1491202 | 0.298273 |
11 | −0.004303392 | 0.04415552 | −0.09745988 | 0.2008131 |
12 | −0.001247142 | 0.02009626 | −0.06205841 | 0.1387547 |
13 | −0.0003449473 | 0.009611985 | −0.0358872 | 0.1028675 |
14 | −8.142579 × 10^{−5} | 0.005287189 | −0.01540058 | 0.08746689 |
15 | −1.157164 × 10^{−5} | 0.003823786 | −0.003026225 | 0.08444066 |
16 | −3.953999 × 10^{−7} | 0.003563991 | −0.000110943 | 0.08432972 |
17 | −5.188181 × 10^{−10} | 0.003554641 | −1.459552 × 10^{−7} | 0.08432957 |
18 | −8.881784 × 10^{−16} | 0.003554628 | −2.498653 × 10^{−13} | 0.08432957 |
19 | 0 | 0.003554628 | 0 | 0.08432957 |
Now the capture area reaches only about 0.01 radians from the solution, so the initial estimate must be a lot more accurate than that to get to the solution quickly. Our initial estimate turns out to have been about that inaccurate, and jumped around seemingly randomly for the first 7 iterations until it happened to get close enough to the solution to home in on that solution after 19 iterations. With a bit of bad luck, it would have taken many more random jumps to get close enough to the solution.
I've tried various initial estimates \( E_1 \), mostly based on the study of extreme approximations described earlier.
In all cases I first replaced \( M \) by \( M \bmod 2π \), so that \( M \) ends up between 0 and 2π radians. In an elliptic orbit, the motion repeats itself when the mean anomaly \( M \) has increased by 2π, so \( M + 2π \) indicates the same location as \( M \) does, but some of the initial estimates do not automatically provide that. By first reducing \( M \) to lie between 0 and 2π we make even those initial estimates repeat themselves.
The formulas are symmetric around \( M = 0 \), so the calculations for \( -M \) are equally easy or difficult as the calculations for \( M \) are, so it is probably better to reduce \( M \) to between −π and +π rather than to between 0 and 2π, but it is interesting to see how large this effect is in practice: How much more effort is it to do the calculations for \( π ≤ M ≤ 2π \) than for \( -π ≤ M ≤ π \), which yield the exact same planet locations?
I tried the following initial estimates for elliptic orbits:
An approximation for small mean anomaly (Equation. \eqref{eq:E1})
\[ E_1 = \frac{M}{|δ|} \]
An approximation for mid-size mean anomaly and near-parabolic orbits (Equation \eqref{eq:E2})
\[ E_1 = (6M)^{1/3} \]
An approximation for large mean anomaly (Equation \eqref{eq:E3})
\[ E_1 = \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) \]
Based on a (hopefully nearby) parabolic orbit
\begin{align*} W & = \frac{3M\sqrt{2}}{4|δ|^{3/2}} \\ u & = (W + \sqrt{W^2 + 1})^{1/3} \\ τ & = u - \frac{1}{u} \\ E_1 & = 2 \arctan\left( τ\sqrt{\frac{1-e}{1+e}} \right) \end{align*}
Based on a (hopefully nearby) parabolic orbit, but with an adjustment for the eccentricity difference
\begin{align*} W & = \frac{3M\sqrt{2}}{4|δ|^{3/2}}\left( 1 - \frac{|δ|}{4} \right) \\ u & = (W + \sqrt{W^2 + 1})^{1/3} \\ τ & = u - \frac{1}{u} \\ E_1 & = 2 \arctan\left(τ \sqrt{\frac{1-e}{1+e}}\right) \end{align*}
Take one step with the first-order method
\[ E_1 = M + e \sin(M) \]
A simple (zeroth-order) estimate for small eccentricity
\[ E_1 = M \]
A formula estimated from the appearance of the graph of \( E \) against \( M \)
\[ E_1 = M + e \left(1 + \frac{1}{2}e^2 \right) \frac{\sqrt{M} \sqrt{2π - M}}{2\sqrt{π}} \]
Another formula estimated from the appearance of the graph of \( E \) against \( M \)
\[ E_1 = M + e \left( 1 + \frac{1}{5}e^2 \right) \frac{\sqrt{M} \sqrt{2π - M}}{2\sqrt{π}} \]
The initial estimate that was used in Section 2.1 is a combination of initial estimates 1 and 2. We refer to that combination as "initial estimate 0".
For hyperbolic orbits
\[ M = e \sinh E - E \]
so for the first-order approximation
\begin{align} g(M,E) & = \arsinh\left( \frac{M + E}{e} \right) \\ g'(M,e) & = \frac{1}{\sqrt{e^2 + (M + E)^2}} \end{align}
so the newest estimate is better than the previous one if \( \frac{1}{\sqrt{e^2 + (M + E)^2}} \lt 1 \), which for hyperbolic orbits (with \( e \gt 1 \)) is always the case. However, if \( e ≈ 1 \) and \( M ≈ E ≈ 0 \), then the improvement that each step brings is very small, so then it takes very many steps before you get near the solution, just like for near-parabolic elliptic orbits.
For the second-order approximation we have
\begin{align} f(E) & = e \sinh E - E \\ f'(E) & = e \cosh E - 1 \\ f''(E) & = e \sinh E \end{align}
and
\begin{align} E_{i+1} & = E_i + \frac{M + E_i - e \sinh E_i}{e \cosh(E_i) - 1} \label{eq:hyp1} \\ E_{i+1} & = E_i + \frac{\frac{M + E_i}{e \cosh(E_i)} - \tanh(E_i)}{1 - \frac{1}{e \cosh(E_i)}} \label{eq:hyp2} \\ β & = \left| \frac{2(e \cosh E - 1)}{e \sinh E} \right| \end{align}
where Equation \eqref{eq:hyp2} is better to use than Equation \eqref{eq:hyp1} is, because Equation \eqref{eq:hyp1} yields denominators and numerators that are too great to fit in calculating machines and most computer programs already for quite modest values of \( E \), and then you cannot calculate a useful correction anymore. Equation \eqref{eq:hyp2} does not have that problem, because there only the denominators feature a fast-rising function.
The least value of \( β \) is attained for \( \cosh(E) = e \); then
\[ β = β_{\text{min}} = 2\frac{\sqrt{e^2 - 1}}{e} \]
For near-parabolic orbits with \( e = 1 + δ \) and \( |δ| ≪ 1 \) we find \( β_{\text{min}} ≈ \sqrt{2δ} \) for \( |E| ≈ \sqrt{2δ} \), so then the initial estimate must be quite accurate already, to be able to locate the solution quickly.
For example, \( \cosh(1) ≈ 1.5 \) and \( \sinh(1) ≈ 1.2 \) but \( \cosh(700) ≈ \sinh(700) ≈ 10^{304} \), so already for a modest rise in input value the output values have become enormous.
If \( M = 3 \), \( e = 2 \), \( E_1 = 700 \), then Equation \eqref{eq:hyp1} yields
\begin{equation*} \begin{split} E_2 & = 700 + \frac{3 + 700 - 2×\sinh(700)}{2×\cosh(700) - 1} \\ & = 700 + \frac{−1.014232×10^{304}}{1.014232×10^{304}} \\ & = 700 - 1 \\ & = 699 \end{split} \end{equation*}
and Equation \eqref{eq:hyp2} yields
\begin{equation*} \begin{split} E_2 & = 700 + \frac{\frac{3 + 700}{2×\cosh(700)} - \tanh(700)}{1 - \frac{1}{2×\cosh(700)}} \\ & = 700 + \frac{\frac{703}{1.014232×10^{304}} - 1}{1 - \frac{1}{1.014232×10^{304}}} \\ & = 700 + \frac{6.93×10^{−302} - 1}{1 - 10^{−304}} \\ & = 700 + \frac{−1}{1} \\ & = 699 \end{split} \end{equation*}
Although all input values have at most 3 digits, the greatest intermediate result has no less than 300 digits. Equation \eqref{eq:hyp1} divides two enormous numbers, but Equarion \eqref{eq:hyp2} divides two numbers near 1. If your calculator or computer thinks that \( \cosh(E_i) \) or \( \sinh(E_i) \) are infinitely large, then Equation \eqref{eq:hyp2} can still yield useful results, because \( \tanh(E) \) never gets greater in magnitude than 1.
I've tried various initial estimates \( E_1 \), mostly based on the study of extreme approximations described earlier. For these hyperbolic orbits I did not force \( M \) to be between 0 and 2π, because hyperbolic orbits do not repeat themselves, so \( M + 2π \) indicates a different position than \( M \) does.
An estimate for small mean anomaly (Equation \eqref{eq:E1}), just like for elliptical orbits
\[ E_1 = \frac{M}{δ} \]
An estimate for mid-size mean anomaly and near-parabolic orbits (Equation \eqref{eq:E2}), just like for elliptical orbits
\[ E_1 = (6M)^{1/3} \]
An estimate for large mean anomaly (Equation \eqref{eq:E3}), just like for elliptical orbits
\[ E_1 = \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) \]
Based on a (hopefully nearby) parabolic orbit, just like for elliptical orbits, but with some extra manipulations because for a hyperbolic orbit some values of \( ν \) are not possible
\begin{align*} W & = \sgn(M) \frac{3\sqrt{2}}{4δ^{3/2}} \min\left( |M|, \frac{4}{3} \right) \\ u & = (W + \sqrt{W^2 + 1})^{1/3} \\ τ & = u - \frac{1}{u} \\ τ_2 & = τ \sqrt{\frac{e - 1}{e + 1}} \\ τ_3 & = \sgn(τ_2) \min(|τ_2|, 1) \\ E_1 & = 2 \artanh(τ_3) \end{align*}
Based on a (hopefully nearby) parabolic orbit, but with an adjustment for the eccentricity difference
\begin{align*} W & = \sgn(M) \frac{3\sqrt{2}}{4δ^{3/2}} \min\left( |M|\left( 1 + \frac{1}{4}δ \right), \frac{4}{3} \right) \\ u & = (W + \sqrt{W^2 + 1})^{1/3} \\ τ & = u - \frac{1}{u} \\ τ_2 & = τ \sqrt{\frac{e - 1}{e + 1}} \\ τ_3 & = \sgn(τ_2) \min(|τ_2|, 1) \\ E_1 & = 2 \artanh(τ_3) \end{align*}
Take one step with the first-order method, from the estimate using Equation \eqref{eq:E1}
\[ E_1 = \arsinh\left( \frac{M + \frac{M}{e - 1}}{e} \right) \]
The initial estimate that was used in Section 2.3 is a combination of initial estimates 1, 2, and 3. We refer to that combination as "initial estimate 0".
I've calculated \( τ \) for a large number of combinations of \( M \) or \( m \) and \( e \), with all initial estimates mentioned above, and have measured how many repeats were necessary to find the true solution to the highest attainable accuracy, and how much calculation time it took.
For \( M \) or \( m \) I've used the following values: 0, 10^{−9}, 10^{−8}, 10^{−7}, 10^{−6}, 10^{−5}, 0.0001, 0.001, 0.01, 0.02π (3.6°) through 1.98π radians (356.4°) in steps of 0.02π (3.6°), and also 10, 100, 1000, 10000, 10^{5}, and 10^{6} radians. For \( e \) I've used the following values: 0, 10^{−6}, 10^{−5}, 0.0001, 0.001, 0.01 through 0.99 in steps of 0.01, then 0.999, 0.9999, 1 − 10^{−5}, 1 − 10^{−6}, 1 − 10^{−7}, 1 − 10^{−8}, 1 − 10^{−9}, 1, 1 + 10^{−9}, 1 + 10^{−8}, 1 + 10^{−7}, 1 + 10^{−6}, 1 + 10^{−5}, 1.0001, 1.001, then 1.01 through 2.00 in steps of 0.01, then 3, 5, 10, 100, 1000, 10000, 10^{5}, 10^{6}.
The measured number of repeats is roughly the same on different computers, because it depends mostly on the chosen data type ― the more precision the data type allows, the more repeats are needed to attain that greater precision. I used "double precision floating point" numbers, i.e., "binary64" according to the IEEE 754 standard. The measured calculation time can be very different on different computers, because it depends also on properties of the computer, such as the CPU's speed.
Table 4 shows a summary of the repeat counts results, if the mentioned anomalies are mean anomalies \( M \). Table 5 shows a similar summary, if the mentioned anomalies are perifocal anomalies \( m \). Table 6 combines both previous tables.
Table 4: Number of Repeats Needed to Solve Kepler's Equation (mean anomaly)
Ellipse (2π) | Ellipse (π) | Hyperbola | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
# | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % |
0 | 10 | 5.1 | 1.4 | 38 | 9 | 4.5 | 1.4 | 54 | 8 | 4.6 | 1.2 | 53 |
1 | (∞) 64 | 6.4 | 3.2 | 17 | (∞) 64 | 5.6 | 3.1 | 29 | 22 | 5.3 | 1.5 | 19 |
2 | 10 | 5.4 | 1.1 | 26 | 9 | 5.0 | 1.1 | 30 | 9 | 5.0 | 0.94 | 34 |
3 | (∞) 22 | 5.5 | 1.4 | 18 | (∞) 22 | 5.3 | 1.6 | 16 | 23 | 4.6 | 1.2 | 58 |
4 | 9 | 5.6 | 1.4 | 13 | 8 | 5.1 | 1.4 | 18 | 8 | 4.7 | 1.0 | 45 |
5 | 9 | 5.7 | 1.4 | 13 | 8 | 5.1 | 1.5 | 19 | 8 | 4.7 | 1.0 | 44 |
6 | (∞) 26 | 4.6 | 1.5 | 76 | (∞) 26 | 4.5 | 1.7 | 63 | 20 | 4.7 | 1.4 | 45 |
7 | (∞) 23 | 5.2 | 1.7 | 35 | (∞) 22 | 5.0 | 1.9 | 30 | ||||
8 | (∞) 21 | 5.4 | 2.0 | 32 | 21 | 4.7 | 1.5 | 46 | ||||
9 | (∞) 26 | 5.3 | 1.9 | 32 | 22 | 4.7 | 1.5 | 46 | ||||
* | 7 | 4.3 | 1.3 | 7 | 3.9 | 1.3 | 7 | 4.0 | 0.99 |
Table 5: Number of Repeats Needed to Solve Kepler's Equation (perifocal anomaly)
Ellipse (2π) | Ellipse (π) | Hyperbola | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
# | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % |
0 | 10 | 4.9 | 1.4 | 35 | 9 | 4.5 | 1.6 | 44 | 10 | 5.0 | 1.4 | 39 |
1 | (∞) 23 | 5.2 | 1.8 | 24 | 9 | 4.5 | 1.6 | 42 | 22 | 5.3 | 1.5 | 26 |
2 | 10 | 5.4 | 1.1 | 16 | 10 | 5.4 | 1.3 | 8 | 10 | 5.4 | 1.1 | 22 |
3 | (∞) 23 | 5.4 | 1.5 | 18 | 9 | 5.2 | 1.4 | 15 | 22 | 5.2 | 1.6 | 39 |
4 | 9 | 4.8 | 1.3 | 39 | 8 | 4.3 | 1.2 | 48 | 8 | 4.6 | 1.2 | 57 |
5 | 9 | 4.9 | 1.4 | 35 | 9 | 4.3 | 1.3 | 47 | 8 | 4.6 | 1.2 | 51 |
6 | (∞) 30 | 5.0 | 1.8 | 47 | 9 | 4.7 | 1.6 | 38 | 20 | 4.9 | 1.5 | 46 |
7 | (∞) 23 | 5.5 | 1.8 | 18 | 10 | 5.1 | 1.6 | 15 | 22 | 5.2 | 1.6 | 39 |
8 | (∞) 21 | 5.0 | 1.5 | 34 | 9 | 4.7 | 1.4 | 33 | 22 | 5.2 | 1.6 | 39 |
9 | (∞) 22 | 5.1 | 1.5 | 31 | 9 | 4.8 | 1.4 | 30 | 22 | 5.2 | 1.6 | 39 |
* | 7 | 3.9 | 1.2 | 6 | 3.6 | 1.2 | 7 | 4.1 | 1.1 |
Table 6: Number of Repeats Needed to Solve Kepler's Equation (combined)
Ellipse (2π) | Ellipse (π) | Hyperbola | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
# | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % | \(X\) | \(μ\) | \(σ\) | % |
0 | 10 | 5.0 | 1.4 | 37 | 9 | 4.5 | 1.5 | 49 | 10 | 4.8 | 1.3 | 46 |
1 | (∞) 64 | 5.8 | 2.6 | 21 | (∞) 64 | 5.0 | 2.5 | 35 | 22 | 5.3 | 1.5 | 22 |
2 | 10 | 5.4 | 1.1 | 21 | 10 | 5.2 | 1.2 | 19 | 10 | 5.2 | 1.1 | 28 |
3 | (∞) 23 | 5.5 | 1.5 | 18 | (∞) 22 | 5.2 | 1.5 | 15 | 23 | 4.9 | 1.5 | 48 |
4 | 9 | 5.2 | 1.4 | 26 | 8 | 4.7 | 1.4 | 33 | 8 | 4.6 | 1.1 | 51 |
5 | 9 | 5.3 | 1.5 | 24 | 9 | 4.7 | 1.5 | 33 | 8 | 4.7 | 1.1 | 48 |
6 | (∞) 30 | 4.8 | 1.6 | 61 | (∞) 26 | 4.6 | 1.6 | 51 | 20 | 4.8 | 1.4 | 45 |
7 | (∞) 23 | 5.3 | 1.7 | 27 | (∞) 22 | 5.1 | 1.8 | 22 | ||||
8 | (∞) 21 | 5.2 | 1.8 | 33 | 21 | 4.7 | 1.4 | 40 | ||||
9 | (∞) 26 | 5.2 | 1.7 | 31 | 22 | 4.8 | 1.5 | 38 | ||||
* | 7 | 4.1 | 1.2 | 7 | 3.8 | 1.3 | 7 | 4.0 | 1.1 |
The columns in the preceding tables provide the following information:
After the first column there follow three parts that each have the same number of columns. The first part "Ellipse (2π)" shows results for elliptic orbits (eccentricity less than 1) based on all anomalies mentioned before, to beyond 2π radians (360°).
The second part "Ellipse (π)" shows results for elliptic orbits based on anomalies from 0 through π radians (180°). For mean anomalies, the range from 0 through 2π radians is equivalent to the range from -π to +π radians, and (also for perifocal anomalies) a calculation for a positive anomaly is as easy or difficult (and yields the same absolute values) as the calculation for the corresponding positive anomaly, so for mean anomalies (but not perifocal anomalies) the anomalies up to π radians suffice.
The third part "Hyperbola" shows results for hyperbolic orbits (eccentricity greater than 1) for all anomalies mentioned before.
For example: Initial estimate number 4 from Table 4 needed at most 9 repeats for elliptic orbits, if all anomalies are taken into account. The average number of repeats was 5.78, with a spread of 1.2. In 12% of cases there was no initial estimate that needed fewer repeats than initial estimate number 4.
Figure 1 shows how many repeats were needed to find the true solution for fixed mean anomaly, for the various initial estimates (identified by the number to the lower right of each picture), and is associated with Table 4. Black means 1 repeat, white means 15 or more repeats, and shades of grey indicate values in between. In each picture, the eccentricity increases from left to right, and the mean anomaly increases from bottom to top. The black vertical line near the middle is for eccentricity equal to 1 (the parabolic orbit), for which always only a single repeat is needed. The left-hand part of each picture is for elliptic orbits (eccentricity less than 1), and the right-hand part is for hyperbolic orbits (eccentricity greater than 1).
Likewise, Figure 2 shows the number of repeats for fixed perifocal anomaly (like Table 5).
A picture like Figures 1 and 2 but then for the calculation time instead of for the number of repeats looks like the corresponding picture for the number of repeats, but with a lot more noise on top of it.
The calculation time that was needed to find the solutions is in general proportional to the number of repeats that was needed to find the solutions: the more repeats, the more calculation time. However, it turns out that on my computer system, with my source of the \( \tan \) function (glibc 2.17), calculating \( \tan(x) \) for \( x \) equal to \( \frac{1}{2}π \) radians (90°) takes hundreds of times more calculation time than for other values of \( x \). The typical calculation time was near 0.5 μs, but the greatest calculation time was near 200 μs. There are similar but less great deviations for certain other combinations of \( M \) and \( e \). If the calculation time is very important to you, then you are advised to run some tests like I did, to see if there are strange cases that you would not predict based on just the repeat counts.
We see that initial estimates numbers 0, 2, 4, and 5 for elliptic orbits (for the chosen anomalies) always provide a solution, but the other initial estimates do not. That is not very strange, because these initial estimates were designed especially for near-parabolic orbits, and those are the most difficult cases for which it is important that the initial estimate be already quite accurate.
If we limit ourselves to anomalies not exceeding π radians (180°), then initial estimate number 1 always provides a solution for \( e ≤ 0.49 \), initial estimate 3 for \( e ≤ 0.99 \) or \( |M| ≥ 10.8° \), initial estimate 6 for \( e ≤ 0.99 \) or \( |M| ≥ 7.2° \), and initial estimate 7 for \( e ≤ 0.97 \) or \( |M| ≥ 28.8° \). Initial estimates numbers 8 and 9 always gave a solution for the chosen anomalies (not exceeding π radians), but those initial estimates need significantly more repeats for the combinations of \( M \) and \( e \) for which initial estimates numbers 3, 6, and 7 find no solution at all, so it seems wise to avoid those combinations also for initial estimates numbers 8 and 9.
None of initial estimates 2, 4, and 5 is obviously better than the others, when judged by the number of repeats, but initial estimate 0 takes on average about 10% fewer repeats. If we look at the calculation time, then initial estimate numbers 0 and 2 are better than numbers 4 and 5, because the calculation of initial estimate numbers 0 and 2 is less complicated than that of numbers 4 and 5 and therefore takes less time.
Judged by the number of repeats, initial estimate number 4 is on average the best for hyperbolic orbits, followed closely by initial estimate numbers 5 and 0. Judged by the calculation time, initial estimate numbers 0 and 2 are better than 4 and 5, just like for elliptic orbits.
We define
\begin{align} δ & ≡ e - 1 \\ f & ≡ \sqrt{\left| \frac{1 + e}{1 - e} \right|} \\ τ & ≡ \tan\left( \frac{1}{2}ν \right) \\ σ & ≡ \frac{τ}{f} \end{align}
where \( e \) is the eccentricity and \( ν \) is the true anomaly. \( δ \) indicates how far away the shape of the orbit is from being a parabola. \( τ \) is a measure for the position of the object in its orbit, just like \( ν \). When \( ν \) grows from −π to +π radians, then \( τ \) grows from −∞ to +∞.
For \( f \) we find the approximations
\begin{align*} e \ll 1 & ⇒ f ≈ 1 + e \\ |δ| \ll 1 & ⇒ f ≈ \sqrt{\frac{2}{|δ|}} \\ e \gg 1 & ⇒ f ≈ 1 + \frac{1}{e} \end{align*}
For \( τ \) we find approximations
\begin{align} |ν| \ll 1 & ⇒ τ ≈ \frac{1}{2}ν \\ ν ≈ ±π + ∆ν & ⇒ τ ≈ -\frac{2}{∆ν} \end{align}
assuming \( ν \) is measured in radians. The other way around,
\begin{align} |τ| \ll 1 & ⇒ ν ≈ 2τ \\ |τ| \gg 1 & ⇒ ν ≈ π\sgn(τ) - \frac{2}{τ} \end{align}
More generally,
\begin{equation} \tan\left( \frac{1}{2}(ν + ∆ν) \right) ≈ τ + \frac{1}{2}(1 + τ^2) ∆ν \end{equation}
Given \( e \), \( τ \), and perifocal distance \( q \), we have for the distance \( r \) to the barycenter of the orbit, and for the rectangular coordinates \( x \) and \( y \)
\begin{align} ρ & ≡ \frac{1 + e}{1 + e + (1 - e)τ^2} = \frac{1}{1 - \left(\frac{τ}{f}\right)^2\sgn(δ)} = \frac{1}{1 - σ^2\sgn(δ)} \\ r & = (1 + τ^2) qρ = q \frac{1 + τ^2}{1 - σ^2\sgn(δ)} \\ x & = (1 - τ^2) qρ = q \frac{1 - τ^2}{1 - σ^2\sgn(δ)} \\ y & = 2τ qρ = q \frac{2τ}{1 - σ^2\sgn(δ)} \end{align}
If \( |δ| \ll 1 \), then \( e ≈ 1 \) and then we find approximations
\begin{align} ρ & ≈ 1 + \frac{1}{2}δτ^2 \\ r & ≈ q\left( 1 + τ^2 + \frac{1}{2}δτ^2(1 + τ^2) \right) \\ x & ≈ q\left( 1 - τ^2 + \frac{1}{2}δτ^2(1 - τ^2) \right) \\ y & ≈ 2qτ\left( 1 + \frac{1}{2}δτ^2 \right) \end{align}
If \( |σ| \ll 1 \), then we find approximations
\begin{align} ρ & ≈ 1 + σ^2\sgn(δ) \\ r & ≈ q(1 + σ^2(\sgn(δ) + f^2)) \\ x & ≈ q(1 + σ^2(\sgn(δ) - f^2)) \\ y & ≈ 2qfσ(1 + σ^2\sgn(δ)) \end{align}
If \( |σ| = 1 - ∆σ ≈ 1 \) and \( 0 ≤ ∆σ \ll 1 \) and \( δ \gt 0 \) (so we're dealing with hyperbolic orbits) then we find approximations
\begin{align} ρ & ≈ \frac{1}{2∆σ} \\ r & ≈ q\frac{1 + f^2}{2∆σ} \\ x & ≈ q\frac{1 - f^2}{2∆σ} \\ y & ≈ q\frac{f}{∆σ}\sgn(σ) \end{align}
We derive approximation formulas by making assumptions about how large \( E \) and \( e \) are, but eventually want to find formulas for calculating \( E \) from \( m \) or \( M \) rather than the other way around. When we have found an approximation formula for a certain range of values of \( E \) and \( e \), then we convert that range of \( E \) into the corresponding range of \( m \) or \( M \). Recalling:
\begin{align} m & = \frac{M}{|δ|^{3/2}} \\ M & = m |δ|^{3/2} \end{align}
The Equation of Kepler for elliptic orbits is
\begin{equation} M = E - e \sin(E) \end{equation}
If \( |E| ≪ 1 \) then
\begin{equation} \sin(E) ≈ E - \frac{1}{6}E^3 \end{equation}
so
\begin{equation} M ≈ E - e \left( E - \frac{1}{6} E^3 \right) = E |δ| + \frac{1}{6}eE^3 \end{equation}
Kepler's Equation for hyperbolic orbits is
\begin{equation} M = e \sinh(E) - E \end{equation}
If \( |E| ≪ 1 \) then
\begin{equation} \sinh(E) ≈ E + \frac{1}{6}E^3 \end{equation}
so
\begin{equation} M ≈ e \left( E + \frac{1}{6} E^3 \right) - E = E |δ| + \frac{1}{6}eE^3 \label{eq:keplersmallE}\end{equation}
just like for elliptic orbits.
On the right-hand side of the last equation, the first term can be greater than, equal to, or less than the second term. The first term dominates if \( E \) is sufficiently small, and otherwise the second term dominates. The turning point is at
\begin{equation} E|δ| = \frac{1}{6}eE^3 ⇔ \frac{6|δ|}{e} = E^2 ⇔ |E| = \sqrt{\frac{6|δ|}{e}} ≈ \sqrt{\frac{|δ|}{e}} \end{equation}
If \( E ≪ \sqrt{\frac{|δ|}{e}} \) then the first term dominates on the right-hand side of Equation \eqref{eq:keplersmallE}, so
\begin{equation} M ≈ E |δ| ⇔ E ≈ \frac{M}{|δ|} \label{eq:MEverysmall} \end{equation}
Which range of \( M \) corresponds to this? If \( E ≪ \sqrt{\frac{|δ|}{e}} \) then
\begin{equation} M ≈ E |δ| ≪ \sqrt{\frac{|δ|}{e}} |δ| = \sqrt{\frac{|δ|^3}{e}} \end{equation}
We also assumed that \( |E| ≪ 1 \), so also
\begin{equation} M ≈ E |δ| ≪ |δ| \end{equation}
So equation \eqref{eq:MEverysmall} is a reasonable approximation if \( M \ll |δ| \) and \( M \ll \sqrt{\frac{|δ|^3}{e}} \). For \( e \lt \frac{1}{2} \) the first restriction is strongest, and otherwise the second one. For a nearly-circular orbit (\(e \ll 1 \)) we must have \( |M| \ll 1 \). For a nearly-parabolic orbit (\( |δ| \ll 1 \)) we must have \( |M| \ll |δ|^{3/2} \). For a strongly hyperbolic orbit (\( e \gg 1 \)) we must have \( |M| \ll e \). A reasonable alternative for all \( e \) is
\begin{equation} |M| \ll \frac{|δ|^{3/2} \left( 1 + |δ|^{3/2} \right) }{1 + δ^2} = \frac{|δ|^{3/2} + |δ|^3}{1 + δ^2} \end{equation}
For elliptic orbits
\begin{equation} τ = f\tan\left( \frac{1}{2}E \right) \label{eq:τe} \end{equation}
and for hyperbolic orbits
\begin{equation} τ = f\tanh\left( \frac{1}{2}E \right) \label{eq:τh} \end{equation}
Because \( |E| \ll 1 \), also
\[ \tan\left( \frac{1}{2}E \right) ≈ \tanh\left( \frac{1}{2}E \right) ≈ \frac{1}{2}E \]
which is much less than 1. We find
\begin{align} τ & ≈ \frac{fM}{2|δ|} \\ σ & ≈ \frac{M}{2|δ|} \end{align}
and \( |τ| \ll 1 \) and \( |σ| \ll \frac{1}{f} \).
If \( \sqrt{\frac{6|δ|}{e}} ≪ |E| ≪ 1 \), then the second term dominates on the right-hand side of Equation \eqref{eq:keplersmallE}, so
\begin{equation} M ≈ \frac{1}{6} eE^3 ⇔ E ≈ \left( \frac{6M}{e} \right)^{1/3} \end{equation}
If \( \sqrt{\frac{|δ|}{e}} ≪ E ≪ 1 \), then also
\begin{equation} \sqrt{\frac{|δ|}{e}} ≪ 1 ⇔ |δ| ≪ e \end{equation}
and that is only possible when \( e ≈ 1 \), so this case is only relevant for near-parabolic orbits, and then
\begin{equation} E ≈ (6M)^{1/3} \end{equation}
Which range of \( M \) corresponds to this? If \( \sqrt{\frac{|δ|}{e}} ≪ |E| ≪ 1 \) and \( |δ| ≪ 1 \) and \( e ≈ 1 \), then
\begin{equation} |δ|^{3/2} \ll |M| \ll 1 \end{equation}
Because \( |E| \ll 1 \)
\[ \tan\left( \frac{1}{2}E \right) ≈ \tanh\left( \frac{1}{2}E \right) ≈ \frac{1}{2} E \]
so
\begin{align} τ & ≈ \frac{1}{2}fE ≈ \frac{1}{2}f(6M)^{1/3} = f \left( \frac{3}{4}M \right)^{1/3} \\ σ & = \frac{τ}{f} ≈ \left( \frac{3}{4}M \right)^{1/3} \end{align}
and \( \sqrt{|δ|} \ll |σ| \ll 1 \) and \( f\sqrt{|δ|} \ll |τ| \ll f \).
For elliptic orbits, the case with \( |E| \gg 1 \) is not relevant, because for an elliptic orbit the positions keep repeating themselves, and you end up with the same position if you add to \( E \) an arbitrary multiple of 2π radians. For hyperbolic orbits, case \( |E| \gg 1 \) is relevant, because for a hyperbolic orbits the positions do not repeat themselves. Here we investigate only hyperbolic orbits (with \( e \gt 1 \)).
If \( |E| \gg 1 \) then
\begin{equation} \sinh(E) ≈ \frac{1}{2} \sgn(E) \exp(|E|) \end{equation}
which is far greater than \( E \) itself, so
\begin{equation} M ≈ \frac{1}{2}e\sgn(E)\exp(|E|) ⇔ E ≈ \sgn(M)\log\left( \frac{2|M|}{e} \right) \end{equation}
This approximating formula is valid when \( |E| \gg 1 \) so when \( |M| \gg e \), but we may want to use it for small \( E \) as well. Then it is undesirable that for \( |M| \lt \frac{1}{2}e \) the calculated value of \( E \) does not have the same sign as \( M \), while in reality those signs are always equal. Where that is useful we fix this by adding 1 to the value that we put into the log-function:
\[ E ≈ \sgn(M) \log\left( 1 + \frac{2|M|}{e} \right) \]
This makes no noticeable difference when \( |E| \gg 1 \) (for which we derived the approximation), but ensures that we always get the same sign for \( E \) as for \( M \).
For hyperbolic orbits
\begin{equation} τ = f\tanh\left( \frac{1}{2}E \right) \end{equation}
If \( |E| \gg 1 \) then
\begin{equation} \tanh\left( \frac{1}{2}E \right) ≈ \sgn(E)(1 - 2\exp(-|E|)) \end{equation}
so
\begin{align} τ & = f\tanh\left( \frac{1}{2}E \right) \notag \\ & ≈ f \sgn(E) (1 - 2\exp(-|E|)) \notag \\ & ≈ f \sgn(M) \left( 1 - 2\exp\left( -\log\left( \frac{2|M|}{e} \right) \right) \right) \notag \\ & = f \sgn(M) \left( 1 - \frac{e}{|M|} \right) \\ σ & ≈ \sgn(M) \left( 1 - \frac{e}{|M|} \right) \end{align}
so when \( M \) gets ever greater then \( |τ| \) gets ever closer to \( f \) and \( |σ| \) gets ever closer to 1. This corresponds to the asymptote of the hyperbola.
(\( ∧ \) means "and"; \( ⇒ \) means "then".)
We distinguish between three areas of validity:
\begin{align} Ⅰ & : |M| \ll |δ| ∧ |M| \ll |δ|^{3/2} \\ Ⅱ & : |δ|^{3/2} \ll |M| \ll e ≈ 1 \\ Ⅲ & : |M| \gg e \gt 1 \end{align}
Area Ⅰ is the circular orbit area: there the orbit looks like a circular one. Similarly, area Ⅱ is the parabolic orbit area, and area Ⅲ is the hyperbolic orbit area.
We found for \( σ \) in the three areas
\begin{align} Ⅰ & ⇒ σ ≈ \frac{M}{2|δ|} & ∧ & |σ| \ll 1 ∧ |σ| \ll \sqrt{|δ|} \\ Ⅱ & ⇒ σ ≈ \left( \frac{3}{4}M \right)^{1/3} & ∧ & \sqrt{|δ|} \ll |σ| \ll 1 \\ Ⅲ & ⇒ σ ≈ \sgn(M) \left( 1 - \frac{e}{|M|} \right) & ∧ & |σ| ≈ 1 \end{align}
and \( τ = fσ \) by definition.
For a nearly-circular orbit, \( e \ll 1 \). Then \( δ ≈ 1 \) and \( f ≈ 1 + e \). Then only area of validity Ⅰ applies, and then \( |M| \ll 1 \) and
\begin{align} σ & ≈ \frac{1}{2}M(1+e) & ∧ & |σ| \ll 1 \\ τ & ≈ \frac{1}{2}M(1+2e) & ∧ & |τ| \ll 1 \\ ν & ≈ M(1+2e) & ∧ & |ν| \ll 1 \end{align}
For a nearly-parabolic orbit we have \( e ≈ 1 \) and \( |δ| \ll 1 \) and \( f ≈ \sqrt{\frac{2}{|δ|}} \). Then areas Ⅰ and Ⅱ apply. For a hyperbolic nearly-parabolic orbit (\( e ≈ 1 ∧ e \gt 1 \)) area Ⅲ applies as well.
For area Ⅰ we find
\begin{align} σ & ≈ \frac{M}{2|δ|} & ∧ & |σ| \ll \sqrt{|δ|} \\ τ & ≈ \frac{M}{\sqrt{2}|δ|^{3/2}} & ∧ & |τ| \ll 1 \\ ν & = \frac{M\sqrt{2}}{|δ|^{3/2}} & ∧ & |ν| \ll 1 \end{align}
For area Ⅱ we find
\begin{align} σ & ≈ \left( \frac{3}{4}M \right)^{1/3} & ∧ & \sqrt{|δ|} \ll |σ| \ll 1 \\ τ & ≈ \sqrt{\frac{2}{|δ|}}\left( \frac{3}{4}M \right)^{1/3} & ∧ & 1 \ll |τ| \ll \sqrt{\frac{2}{|δ|}} \\ ν & ≈ \sgn(τ) \left( π - \frac{2}{|τ|} \right) & ∧ & π - \sqrt{|δ|} \ll |ν| \ll π \end{align}
For area Ⅲ (only for hyperbolic orbits) we find
\begin{align} σ & ≈ \sgn(M) \left( 1 - \frac{1}{|M|} \right) & ∧ & |δ| \ll |σ| \lt 1 \\ τ & ≈ \sgn(M) \sqrt{\frac{2}{|δ|}} \left( 1 - \frac{1}{|M|} \right) & ≈ & \sqrt{\frac{2}{|δ|}} \\ ν & ≈ ν_∞ \sgn(τ) \left( 1 - \frac{2}{|τ|} \right) & ≈ & ν_∞ \sgn(τ) \\ ν_∞ & ≈ \sgn(τ) \left( π - \sqrt{2|δ|} \right) & ≈ & π \end{align}
For a given value of \( M \) (\( |M| \ll 1 \)) the value of \( τ \) increases when \( δ \) gets closer to zero; if \( δ = 0 \) (for a parabolic orbit) then \( τ \) has become infinitely large. This means that for a parabolic orbit \( M \) has no meaning, and for a nearly-parabolic orbit it has little meaning. There is a good alternative.
The relationship between time \( t \) and mean anomaly \( M \) is
\begin{equation} M = t \sqrt{\frac{Γ}{a^3}} \end{equation}
Because \( q = a|1 - e| = a|δ| \), also \( a = \frac{q}{|δ|} \), so
\begin{align} M & = t \sqrt{\frac{Γ}{q^3}} |δ|^{3/2} \\ m & ≡ \frac{M}{|δ|^{3/2}} = t \sqrt{\frac{Γ}{q^3}} \end{align}
I call \( m \) the perifocal anomaly. For an orbit with a given perifocal distance \( q \), the perifocal anomaly \( m \) depends on the time but not on \( δ \), on how much the orbit looks like a parabola. In terms of the perifocal anomaly \( m \), the areas of validity for nearly-parabolic orbits are:
\begin{align} Ⅰ & : |m| \ll 1 \\ Ⅱ & : 1 \ll |m| \ll \frac{1}{|δ|^{3/2}} \\ Ⅲ & : |m| \gg \frac{1}{|δ|^{3/2}} ∧ e \gt 1 \end{align}
and the approximations for \( τ \) are:
\begin{align} Ⅰ & : \frac{m}{\sqrt{2}} & ∧ & |τ| \ll 1 \\ Ⅱ & : \frac{(6m)^{1/3}}{\sqrt{2}} & ∧ & |τ| \gg 1 \\ Ⅲ & : \sgn(m) \sqrt{\frac{2}{|δ|}} \left( 1 - \frac{1}{|m||δ|^{3/2}} \right) & ∧ & |τ| \gg 1 \end{align}
If \( δ \) gets ever closer to zero, then the boundary between validity areas Ⅱ and Ⅲ shifts to ever greater values of \( m \), finally disappearing into infinity for \( δ \) equal to zero.
This fits nicely with the precise formulas for the parabolic orbit (with \(e = 1 \)). With those formulas we can calculate \( τ \) precisely, but it is not easy to see how those formulas behave for very small or very large anomaly. We derive approximations for those cases.
\begin{align} W & ≡ t \sqrt{\frac{9Γ}{8q^3}} = \sqrt{\frac{9}{8}} m = \frac{3\sqrt{2}}{4} m \\ |W| \ll 1 & ⇒ τ ≈ \frac{2}{3}W = \frac{m}{\sqrt{2}} \\ |W| \gg 1 & ⇒ τ ≈ (2W)^{1/3} = \frac{(6m)^{1/3}}{\sqrt{2}} \end{align}
which are the same approximation formulas that we found also for near-parabolic orbits. Near-parabolic orbits therefore morph seamlessly into a true parabolic orbit when \( δ \) gets ever smaller and finally becomes equal to 0.
For a strongly hyperbolic orbit, \( e \gg 1 \) and \( f ≈ 1 + \frac{1}{e} \). Then only validity areas Ⅰ and Ⅲ apply. For area Ⅰ we find
\begin{align} σ & ≈ \frac{M}{2|δ|} ≈ \frac{M}{2e} & ∧ & |σ| \ll 1 \\ τ & ≈ \frac{M}{2e} \left( 1 + \frac{1}{e} \right) & ∧ & |τ| \ll 1 \\ ν & ≈ \frac{M}{e} \left( 1 + \frac{1}{e} \right) & ∧ & |ν| \ll 1 \end{align}
For area Ⅲ we find
\begin{align} σ & ≈ \sgn(M) \left( 1 - \frac{e}{|M|} \right) & ∧ & |σ| ≈ 1 \\ τ & ≈ \sgn(M) \left( 1 + \frac{1}{e} - \frac{e}{|M|} \right) & ∧ & |τ| ≈ 1 \\ ν & ≈ \sgn(M) \left( π - \frac{2}{1 + \frac{1}{e} - \frac{e}{|M|}} \right) & ∧ & |ν| ≈ π \end{align}
http://aa.quae.nl/en/reken/kepler.html;
Last updated: 2016−02−07