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?
-
$\begingroup$ It looks like they just calculated the julian dates and subtracted them. What seems confusing about it? $\endgroup$– PaulSep 22, 2017 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$– user10509Sep 22, 2017 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 HammenSep 23, 2017 at 12:06
2 Answers
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.
-
$\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$– RustySep 23, 2017 at 4:05
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
-
-
$\begingroup$ There is a problem here! As noted in this comment some of these were mis-transcribed. I'll fix this today after coffee... $\endgroup$– uhohJul 26, 2020 at 0:49