Hyper/parabolic Kepler orbits and "mean anomaly"

Question

In an elliptical kepler orbit there is an easy recipe to describe the motion/position of a satellite at time $t$. One just follows the following steps - an important detail for me is that the numerical part has always the same error, the error doesn't increase with time, nor is it based on the timestep. (Which it would if one would solve a differential equation).

Calculate mean motion $n$ (given a semi major axis $a$, and central mass $M$, this is basically the 2 pi divided by the period): $$n = \sqrt{\frac{GM}{a^3}}$$

Now the definition of mean anomaly:

$$M(t) = M_0 + n \cdot t$$

And one can (numerically) solve the following equation, to get the eccentric anomaly $E$:

$$M = E - \varepsilon \sin(E)$$

And then simple geometry allows one to find the true anomaly $\theta$.

---

Now this is all fine and good working. I'm at a loss as how to to apply this to non-elliptical orbits. The main problems I have are that there's no notion of semi major axis (or when using a negative for hyperbolic orbits- it is infinite for parabolae, which would result in a mean motion of 0.).

Nor is there an meaning for the eccentric anomaly, non elliptical orbits don't have a geometric center. (The true anomaly would still be defined as the focal point).

Answer 1

For the case of hyperbolic orbit, you use the hyperbolic trigonometric function. So your equation for Hyperbolic Mean anomaly becomes:

$$M_h = e\sinh F - F$$

where $F$ is your hyperbolic eccentric anomaly, which is analogous to eccentric anomaly for ellipse and $M_h$ is still given by the same expression. The angle measurement for hyperbolic trigonometric identities are not the same but analogous, so using the center of hyperbola which exists between the real hyperbola and it's image, you get the sinh function value. Assuming the object to be at distance $r$ from the primary focus and center of hyperbola being the reference origin $\sinh F$ will be given by:

$$\sinh F = \dfrac{y}{b}$$

where $y$ is the perpendicular distance from the focal axis of the object on trajectory, and $b$ is the semi-minor axis of the conic.

As for the semi-major axis it can be obtained by orbit equation at periapse for hyperbola which is $a = \dfrac{h^2}{\mu} \dfrac{1}{e^2-1}$, which is positive.

For parabolic orbit you do not define Mean Anomaly, because of rightly mentioned absence of geometric center, instead use Baker's equation $2\sqrt{\dfrac{\mu}{p^3}}\left(t-\tau\right) = \tan\dfrac{\theta}{2} +\dfrac{1}{3}\tan^3\dfrac{\theta}{2}$ where $\tau$ is your periapse passage time and $t$ is the elapsed time since periapse, and $p$ is the semi-latus rectum which is well defined for a parabola. There are solutions available to this equation that gives the position of the object at any given time given some input parameters.

Hope this helps.

Answer 2

Starting from this answer, Here are the forms for elliptical, parabolic, and hyperbolic equivalents just to have everything in one place. I've transcribed from this handy table and added alternate forms for $E$ and $F$ from Wikipedia. I did some spot checking with numerical integration just to make sure the results agreed.

Ellipse, Circle $(0 \le \epsilon < 1)$:

$$\tan \frac E2 = \sqrt{\frac{1-e}{1+e}} \tan\frac\theta2$$ $$\text{or}$$ $$\cos E = \frac{e+\cos\theta}{1+e\cos\theta} \tan\frac\theta2,$$ $$M = E - e\sin E,$$ $$\Delta M = M_2 - M_1,$$ $$\Delta t = \sqrt{\frac{a^3}{\mu}} \Delta M,$$

Hyperbola ($\epsilon > 1)$:

$$\tanh \frac F2 = \sqrt{\frac{e-1}{e+1}} \tan\frac\theta2$$ $$\text{or}$$ $$\cosh F = \frac{e+\cos\theta}{1+e\cos\theta} \tan\frac\theta2,$$ $$M = e\sinh F-F,$$ $$\Delta M = M_2 - M_1,$$ $$\Delta t = \sqrt{\frac{(-a)^3}{\mu}} \Delta M,$$

Parabola($\epsilon = 1)$:

$$D = \tan\frac\theta2,$$ $$M = D + \frac{D^3}{3},$$ $$\Delta M = M_2 - M_1,$$ $$\Delta t = \sqrt{\frac{q^3}{\mu}} \Delta M,$$

---

To get the semi-major axis $a$ or to get $q$, use the following (don't worry that $a$ is negative for the hyperbola):

Ellipse, Hyperbola:

$$a=\frac{r_{peri}}{1-e}$$

Ellipse:

$$a=\frac{r_{peri}+r_{apo}}{2}$$

Circle:

$$a=r$$

Parabola:

$$q=r_{peri}$$

A quick check with $\mu=1$:

 e      theta      a     v_peri    E/D/F       M         t
1.5   90.000000  -2.0   1.581139  55.14281  40.94513  2.021271
1.0   90.000000   n/a   1.414214  57.29578  76.39437  1.885618
0.5   90.000000   2.0   1.224745  60.00000  35.19020  1.737177
0.0   90.000000   1.0   1.000000  90.00000  90.00000  1.570796

If you want to try it in Python:

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc  = -mu * x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

def get_D(theta, e):
    if e == 1.0:
        D    = np.tan(0.5*theta)
    else:
        D    = np.nan
    return D

def get_E(theta, e):
    if e < 1.0:
        term = np.sqrt((1.-e)/(1.+e)) * np.tan(0.5*theta)
        E    = 2.*np.arctan(term)
    else:
        E    = np.nan
    return E

def get_E_alt(theta, e):
    if e < 1.0:
        term = (e + np.cos(theta)) / (1. + e*np.cos(theta))
        E    = np.arccos(term)
    else:
        E    = np.nan
    return E

def get_F(theta, e):
    if e > 1.0:
        term = np.sqrt((e-1.)/(e+1.)) * np.tan(0.5*theta)
        F    = 2.*np.arctanh(term)
    else:
        F    = np.nan
    return F

def get_F_alt(theta, e):
    if e > 1.0:
        term = (e + np.cos(theta)) / (1. + e*np.cos(theta))
        F    = np.arccosh(term)
    else:
        F    = np.nan
    return F

def get_M_from_E(E, e):
    if e < 1.0:
        M = E - e*np.sin(E)
    else: 
        M = np.nan
    return M

def get_M_from_F(F, e):
    if e > 1.0:
        M = e*np.sinh(F) - F
    else: 
        M = np.nan
    return M

def get_M_from_D(D, e):
    if e == 1.0:
        M = D + D**3/3.
    else: 
        M = np.nan
    return M

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# http://www.bogan.ca/orbits/kepler/orbteqtn.html

quarterpi, halfpi, pi, twopi = [f*np.pi for f in [0.25, 0.5, 1, 2]]
rads, degs = pi/180, 180/pi

mu = 1.0

th0, th1 = 0.0, halfpi
print "th0, th1 (degs): ", degs*th0, degs*th1

eccs = [1.5, 1.0, 0.5, 0.0]

for e in eccs:

    print "e: ", e

    rp =  1.0  # periapsis

    if e < 1.0:
        print "     is ellipse!"

        ra = rp * (1+e)/(1-e)
        print "rp, ra: ", rp, ra

        a0 = 0.5*(rp + ra)
        v0 = np.sqrt(mu * (2./rp - 1./a0))
        print "a0, v0: ", a0, v0

        E0,  E1  = [get_E(th, e) for th in [th0, th1]]
        M0,  M1  = [get_M_from_E(E, e)  for E  in [E0,  E1 ]]
        print "E0, E1 (degs): ", degs*E0, degs*E1
        print "M0, M1 (degs): ", degs*M0, degs*M1

        print "E0, E1: ", E0, E1
        print "M0, M1: ", M0, M1

        dt = np.sqrt(a0**3/mu) * (M1-M0)

        print "dt (sec): ", dt

    elif e > 1.0:
        print "     is hyperbola!"

        ra = rp * (1+e)/(1-e)
        print "rp, ra: ", rp, ra

        a0 = 0.5*(rp + ra)
        v0 = np.sqrt(mu * (2./rp - 1./a0))
        print "a0, v0: ", a0, v0

        F0,  F1  = [get_F(th, e) for th in [th0, th1]]
        M0,  M1  = [get_M_from_F(F, e)  for F  in [F0,  F1 ]]
        print "F0, F1 (degs): ", degs*F0, degs*F1
        print "M0, M1 (degs): ", degs*M0, degs*M1

        print "F0, F1: ", F0, F1
        print "M0, M1: ", M0, M1

        dt = np.sqrt((-a0)**3/mu) * (M1-M0)

        print "dt (sec): ", dt

    elif e == 1.0:
        print "     is parabola!"

        print "rp: ", rp

        v0 = np.sqrt(mu * (2./rp))
        print "v0: ", v0

        D0,  D1  = [get_D(th, e) for th in [th0, th1]]
        M0,  M1  = [get_M_from_D(D, e)  for D  in [D0,  D1 ]]
        print "D0, D1 (degs): ", degs*D0, degs*D1
        print "M0, M1 (degs): ", degs*M0, degs*M1

        print "D0, D1: ", D0, D1
        print "M0, M1: ", M0, M1

        q = rp

        dt = np.sqrt(2.*q**3/mu) * (M1-M0)

        print "dt (sec): ", dt

    time = np.array([0, dt])
    X0   = np.array([rp, 0, 0, v0])

    answer, info = ODEint(deriv, X0, time, atol=1E-13, rtol=1E-13, full_output=True)

    x, y, vx, vy = answer.T
    theta = np.arctan2(y, x)

    print degs*theta[0], degs*theta[-1], " should be ", degs*th0, degs*th1