How to calculate the flight path angle, γ, from a state vector?

@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

for r, v in zip(rstarts, vstarts):
X0 = np.array([r, 0, 0, v])
answer, info = ODEint(deriv, X0, time, full_output= True)

gammas = []
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()

• should this question be TLDRed to indicate it was a coding error, as it still seems to asking what's wrong with the formula
– user20636
Jan 21 '20 at 9:03

This is a problem that has plagued groups of people very knowledgeable about orbital dynamics but who learned using different textbooks: there are two different definitions of "flight path angle"!!

In addition to $\gamma$, the angle between the tangential direction and the velocity vector, there is $\beta$, the angle between the radial direction and the velocity vector. People often say "flight path angle" without saying which definition they're using. Confusing! (I just noticed that the diagram in Julio's answer also shows $\beta$)

If you work with $\beta$ instead of $\gamma$, $\beta$ is given by

$$\arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{1}$$

which goes from 0 ("straight up") to $\pi$ ("straight down"). Using $\gamma$, "straight up" is $\pi/2$ and "straight down" is $-\pi/2$, so converting $\beta$ to $\gamma$ you just subtract $\beta$ from $\pi/2$:

$$\gamma = \pi/2 - \arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{2}$$

This is equivalent to

$$\gamma = \arcsin\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{3}$$

I'm not familiar with the language you used for your calculations and plots, so I haven't looked at your algorithm to see why there are "too many wiggles".

• Thanks! I've added tags (numbers) to equations. Would you say there are too many wiggles, or is that wiggling behavior in fact reasonable? Since your $\beta$ (eq 1) is the same as my erroneous $\gamma$ except for an offset of half pi, then the wiggles in my plot should be the same as those in a proper plot of your $\beta$ (eq 1).
– uhoh
Jul 15 '18 at 0:02
• Looks like too many wiggles to me. I'll check that later. Jul 15 '18 at 0:54
• @uhoh , Actually, my eq 1 is just the negative of your equation. Something else is wrong. Of course you know something is wrong because all the plotted $\gamma$s are negative or zero, which can't be except for an inward spiral. For a Keplerian eccentric orbit $\gamma$ should cross zero exactly twice, at periapsis and apoapsis, and be monotonic between the extrema, for both the short (extremum through periapsis to the other extremum) and long (extremum through apoapsis to the other extremum) segments. I'll see if I can draw an example of what the $\gamma$ curve should look like. Jul 15 '18 at 5:16
• Oops, I should have said above, "My eq. 2 is just the negative of yours." I should log off and go to bed! Jul 15 '18 at 5:48
• @uhoh "tangential" to a sphere centered at the primary, not to the orbit. Personally I'd prefer to say "lateral velocity" but my first Orbital Dynamics professor at Stanford used "tangential". Oct 9 '19 at 0:08

I found the error in the script, it was due to my "homebrew" dot product. I had an extra squaring:

dotted = ((xx*vv)**2).sum(axis=0)     # WRONG

dotted = (xx*vv).sum(axis=0)          # Correct


So using this plus @TomSpilker's excellent clarifications I've use the following two methods to calculate gamma:

Method 1:

$$\gamma_1 = \arcsin\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{3}$$

Method 2:

A brute-force alternate method to double check:

$$\theta_r = \arctan2(y, x)$$

$$\theta_v = \arctan2(vy, x)$$

$$\theta_{tanj} = \theta_r + \frac{\pi}{2}$$

$$\gamma_2 = \theta_{tanj} - \theta_v$$

$$\gamma_{2mod} = \mod(\gamma_2+ \pi, 2\pi) - \pi$$

The modulo operation is only really needed in the computer program since each theta comes from a separate arctan2 operation:

gammas_1, gammas_2 = [], []

xx, vv = a.reshape(2, 2, -1)

dotted = (xx*vv).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma_1 = np.arcsin(dotted/(rabs*vabs))   # Per Tom Spilker's answer Eq. 3

theta_r = np.arctan2(xx[1], xx[0])
theta_v = np.arctan2(vv[1], vv[0])
theta_tanj = theta_r + halfpi

gamma_2 = theta_tanj - theta_v
gamma_2 = np.mod(gamma_2 + pi, twopi) - pi

gammas_1.append(gamma_1)
gammas_2.append(gamma_2)

plt.figure()
plt.subplot(2, 1, 1)
for gamma_1 in gammas_1:
plt.plot(time, gamma_1)
plt.title('gammas_1', fontsize=16)

plt.subplot(2, 1, 2)
for gamma_2 in gammas_2:
plt.plot(time, gamma_2)
plt.title('gammas_2', fontsize=16)

• Indeed the new $\gamma$ plot is what I expected. Hooray! Good sleuthing. Jul 15 '18 at 17:34