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

Question

@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()

Answer 1

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".

Answer 2

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 = [], []
for a in answers:

    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)