3
$\begingroup$

I was working an example in Charles Brown's "Elements of Spacecraft Design" but I can not figure out how the Time of Flight was calculated. Here is the problem Example 6.5. I got the correct answer when calculating a, rP, and e, but after that I haven't been successful. Could anyone help?

$\endgroup$
  • $\begingroup$ It looks like they just calculated the julian dates and subtracted them. What seems confusing about it? $\endgroup$ – Paul Sep 22 '17 at 2:58
  • $\begingroup$ Can you edit and include the problem in here? We shouldn't have to go somewhere else, and it will help against link rot. $\endgroup$ – Jan Doggen Sep 22 '17 at 7:30
  • $\begingroup$ I rolled the question back to the original, Rusty. Having two pages from the referenced book as images violates both copyright and accessibility concerns. $\endgroup$ – David Hammen Sep 23 '17 at 12:06
2
$\begingroup$

How to compute time of flight

It's easy if you know the gravitational parameter $\mu\equiv GM$, the semi-major axis $a$, the eccentricity $e$, the initial true anomaly $\theta_0$, and the final true anomaly $\theta_1$ (or equivalently, the change in true anomaly $\Delta \theta \equiv \theta_1 - \theta_0$).

Step 1: Compute the initial and final eccentric anomalies from the eccentricity and the initial and final true anomalies via $$\tan \frac E2 = \sqrt{\frac{1-e}{1+e}} \tan\frac\theta2$$ Step 2: Compute the initial and final mean anomalies from the eccentricity and the initial and final eccentric anomalies via Kepler's equation $$M = E - e\sin E$$ Step 3: Compute the change in mean anomaly via $$\Delta M = M_2 - M_1$$ Step 4: Compute the mean motion from the gravitational parameter and semi-major axis length via $$n = \sqrt{\frac \mu {a^3}}$$ Step 5: Compute the time of flight from the mean motion and the change in mean anomaly via $$\Delta t = \frac {\Delta M} n$$


How to design a transfer orbit

This is considerably more challenging. Here the goal is to find the transfer orbit given the initial and final position of the object. This is Lambert's problem. The underlying equations are transcendental; there is no nice closed form solution. Instead you need to make a guess and then iteratively refine the guess. There are a number of ways to iteratively solve Lambert's problem.

$\endgroup$
  • $\begingroup$ In the example above for the 'Theta at Earth = 180', Theta_1 = 180 and Theta_2 = 312.99. I would then calculate E and M for both and use those from then on right? I am doing that but I don't get 116.15 days. I am getting 168. $\endgroup$ – Rusty Sep 23 '17 at 4:05
2
$\begingroup$

edit: There's more about the math for hyperbolic orbits in @DavidHammon's other answer.

I'll add to @DavidHammon's clear answer by including the 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$ and $r_{peri}=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
$\endgroup$
  • $\begingroup$ edit suggestions welcomed! $\endgroup$ – uhoh Sep 22 '17 at 6:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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