@Julio's excellent answer describes a flight path angle, and explains that it is the angle between the tangential direction (perpendicular to the radial vector to the central body) and the current velocity vector.
I've first tried to get the angle from this expression, but it's obviously wrong, since $\arccos$ is an even function and the angle can go from $-\pi/2$ to $\pi/2$:
$$\arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) - \frac{\pi}{2} \ \ \ \text{ (incorrect!)}$$
I've integrated orbits for GM ($\mu$) and SMA ($a$) of unity and starting distances from 0.2 to 1.8. That makes the period always $2 \pi$. When I plot the result of my function, I get too many wiggles.
What expression can I use to get the correct flight path angle gamma starting from state vectors?
Revised python for the erroneous part would be appreciated, but certainly not necessary for an answer.
def deriv(X, t):
x, v = X.reshape(2, -1)
acc = -x * ((x**2).sum())**-1.5
return np.hstack((v, acc))
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
T = twopi
time = np.linspace(0, twopi, 201)
a = 1.0
rstarts = 0.2 * np.arange(1, 10)
vstarts = np.sqrt(2./rstarts - 1./a) # from vis-viva equation
answers = []
for r, v in zip(rstarts, vstarts):
X0 = np.array([r, 0, 0, v])
answer, info = ODEint(deriv, X0, time, full_output= True)
answers.append(answer.T)
gammas = []
for a in answers:
xx, vv = a.reshape(2, 2, -1)
dotted = ((xx*vv)**2).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma = np.arccos(dotted/(rabs*vabs)) - halfpi
gammas.append(gamma)
if True:
plt.figure()
plt.subplot(4, 1, 1)
for x, y, vx, vy in answers:
plt.plot(x, y)
plt.plot(x[:1], y[:1], '.k')
plt.plot([0], [0], 'ok')
plt.title('y vs x')
plt.subplot(4, 1, 2)
for x, y, vx, vy in answers:
plt.plot(time, x, '-b')
plt.plot(time, y, '--r')
plt.title('x (blue) y (red, dashed)')
plt.xlim(0, twopi)
plt.subplot(4, 1, 3)
for x, y, vx, vy in answers:
plt.plot(time, vx, '-b')
plt.plot(time, vy, '--r')
plt.title('vx (blue) vy (red), dashed')
plt.xlim(0, twopi)
plt.subplot(4, 1, 4)
for gamma in gammas:
plt.plot(time, gamma)
plt.title('gamma?')
plt.xlim(0, twopi)
plt.show()