Ctrl+K

Example: Hyperbolic Trajectory

Example: Hyperbolic Trajectory#

A geocentric trajectory has a perigee altitude of 300 km and a perigee velocity of 15 km/s. Calculate the time to fly from perigee to a true anomaly of \(\nu =\) 100°, and the position at that time. Then, calculate the true anomaly and speed 3 hr later.

Given True Anomaly, Find Time Since Perigee#

As for the elliptical case, the solution has three steps:

  1. Find the hyperbolic eccentric anomaly, \(F\), from the true anomaly, \(\nu\)

  2. Find the hyperbolic mean anomaly, \(M_h\), from the eccentric anomaly

  3. Find the time since perigee, \(t\), from the mean anomaly

Eq. (228) gives the eccentric anomaly in terms of the true anomaly. The only unknown parameter is the eccentricity of the hyperbola, which we need to find from the given orbital elements. Since we have \(v_p\) and \(r_p\) in the problem statement, we can calculate the angular momentum followed by the eccentricity.

import numpy as np
from scipy.optimize import newton

mu = 3.986004418E5  # km**3/s**2
nu_1 = np.radians(100)
r_p = 300 + 6378.1  # km
v_p = 15  # km/s
h = r_p * v_p  # km**2/s
e = h**2 / (r_p * mu) - 1

The eccentricity is \(e =\) 2.7696. Since \(e > 1\), this trajectory is a hyperbola.

We should find the true anomaly of the asymptote from Eq. (150), to ensure that our desired true anomaly is valid.

nu_infty = np.arccos(-1 / e)

The true anomaly of the asymptote is \(\nu_{\infty} =\) 111.17°. Therefore, our desired true anomaly is valid. Now we can calculate the eccentric anomaly, \(F\).

F_1 = 2 * np.arctanh(np.sqrt((e - 1)/(e + 1)) * np.tan(nu_1 / 2))

Then, the mean anomaly is found from Kepler’s equation, Eq. (229):

M_h1 = e * np.sinh(F_1) - F_1

Finally, calculating the time from the mean anomaly is done from the definition of the mean anomaly, Eq. (220):

t_1 = h**3 / mu**2 * 1 / (e**2 - 1)**(3/2) * M_h1

The total time is \(t_1 =\) 1.15 hr.

Given Time Since Perigee, Find True Anomaly#

Now, let’s calculate the true anomaly 3 hours later, after about 4 total hours since perigee have elapsed. Again, there are three steps:

  1. Given time since perigee, \(t\), calculate the hyperbolic mean anomaly, \(M_h\)

  2. Calculate the hyperbolic eccentric anomaly, \(F\), from the hyperbolic mean anomaly

  3. Calculate the true anomaly, \(\nu\), from the hyperbolic eccentric anomaly

Since we already have the orbital eccentricity and specific angular momentum, we can start by finding the mean anomaly at the time.

t_2 = 3 * 3600 + t_1  # sec
M_h2 = mu**2 / h**3 * (e**2 - 1)**(3/2) * t_2

Now, we need to solve Kepler’s equation to find the eccentric anomaly, \(F\). Since the equation is transcendental in \(F\), we need to use the Newton solver in SciPy. Since we know the derivative, we will define two Python functions:

  1. Kepler’s equation, \(f(F) = 0\)

  2. The derivative of Kepler’s equation with respect to \(F\), \(f'(F)\)

def kepler(F, M_h, e):
    """Kepler's equation, to be used in a Newton solver."""
    return e * np.sinh(F) - F - M_h

def d_kepler_d_F(F, M_h, e):
    """The derivative of Kepler's equation, to be used in a Newton solver.
    
    Note that the argument M_h is unused, but must be present so the function
    arguments are consistent with the kepler function.
    """
    return e * np.cosh(F) - 1

F_2 = newton(func=kepler, fprime=d_kepler_d_F, x0=np.pi, args=(M_h2, e))

With this value for \(F\), we can calculate the value for \(\nu\). To avoid quadrant ambiguity problems, we will use Eq. (228).

sqrt_e = np.sqrt((e + 1) / (e - 1)) 
nu_2 = (2 * np.arctan(sqrt_e * np.tanh(F_2 / 2))) % (2 * np.pi)

Like for the ellipse, to convert \(\nu\) to the range \([0, 2\pi)\), we take the modulus with \(2\pi\). In most programming languages, Python and MATLAB included, the arctan function returns a value between \(-\pi/2\) and \(\pi/2\). When the result is multiplied by 2, it gives the range from \(-\pi\) to \(\pi\). We need to transform this angle to be in the range of \(0\) to \(2\pi\). To do so, we can take the modulus of the angle with \(2\pi\).

The modulus is the remainder after division. In Python, the modulus operator is %, while in MATLAB, we have to use the function mod(numerator, denominator). This works for both positive and negative numbers, and ensures that we get the correct angle for the appropriate quadrant.

The true anomaly after 4.15 hr is \(\nu_2 =\) 107.78°.

Calculate the Speed of the Spacecraft#

To find the speed, we will calculate the velocity components. The radius at \(\nu_2 =\) 107.78° can be found from the orbit equation, Eq. (113).

r_2 = h**2 / mu / (1 + e * np.cos(nu_2))

The velocity components can be found from Eqs. (114) and (115).

v_perp = h / r_2
v_r = mu / h * e * np.sin(nu_2)
v_2 = np.sqrt(v_r**2 + v_perp**2)

The radius is \(r_2 =\) 1.6318E+05 km and the speed is \(v_2 =\) 10.51 km/s.

MATLAB Solution#

In MATLAB, the following code will give the same result:

function kepler
    mu = 3.986e5; % km^3/s^2
    r_p = 300 + 6378; % km
    v_p = 15; % km/s
    h = r_p * v_p;
    e = h^2 / (mu * r_p) - 1;
    
    nu_1 = deg2rad(100);
    F_1 = 2 * atanh(sqrt((e - 1) / (e + 1)) * tan(nu_1 / 2));
    M_h1 = e * sinh(F_1) - F_1;
    t_1 = h^3 / mu^2 * 1 / (e^2 - 1)^(3 / 2) * M_h1;
    
    t_2 = t_1 + 3 * 3600;
    M_h2 = mu^2 / h^3 * (e^2 - 1)^(3 / 2) * t_2;

    function x = fun(F, M_h, e)
        x = e * sinh(F) - F - M_h;
    end

    F_2 = fzero(@(x) fun(x, M_h2, e), [3, 4]);
    t2 = 2 * atan(sqrt((e + 1) / (e - 1)) * tanh(F_2 / 2));
    nu_2 = mod(t2, 2 * pi);
    disp(rad2deg(nu_2))
    
    r_2 = h^2 / mu * 1 / (1 + e * cos(nu_2));
    v_perp = h / r_2;
    v_r = mu / h * e * sin(nu_2);
    v = sqrt(v_perp^2 + v_r^2);
    disp([v_perp, v_r, v])
end

We are using fzero() again to solve Kepler’s equation. I’m not sure how sensitive fzero() will be to the initial guess.

Contents