NRHO 9:2 orbit at Apolune [AWP State]

[rslippert]

#others using AWP might enjoy this NRHO 9:2 Halo orbit state to play with.

gatway92 = [1.021325, 0.0, -0.181619, 0.0 , -0.101736, 0.0 ]# semi-stable Apolune NRHO 9:2 orbit state

[rslippert]

Validation of the Purdue JNRHO 9:2 Gateway orbit in 64 lines of Python code
from orbital-mechanics.space/the-n-body-problem/Equations-of-Motion-CR3BP.html
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

m_1 = 5.974E24 # kg Earth mass 5.97219E24
m_2 = 7.348E22 # kg Moon mass 7.34767309E22
mu = m_2/(m_1 + m_2)
x_0 = 1 - mu
y_0 = .0455
z_0 = 0
vx_0 = -0.5
vy_0 = 0.5
vz_0 = 0
hstate = np.array([1-mu,.0455,0,-0.5,0.5,0])
NHRO92 = np.array([1.0242, 0.0, -0.1815,0.0, -0.105,0])
JNHRO = np.array([1.021335, 0.0, -0.181619,0.0, -0.101756,0])
Y_0 = JNHRO

-------------------------------------------------------------------------------

def nondim_cr3bp(t, state):
x, y, z, xv, yv, zv = state # Get the position and velocity from state
state_dot = np.zeros_like(state) # Define the derivative vector
state_dot[:3] = state[3:]
sigma = np.sqrt(np.sum(np.square([x + mu, y, z])))
psi = np.sqrt(np.sum(np.square([x - 1 + mu, y, z])))
state_dot[3] = 2 * yv + x - (1 - mu) * (x + mu)
/ sigma3 - mu * (x - 1 + mu) / psi3
state_dot[4] = -2 * xv + y - (1 - mu) * y / sigma3 - mu * y / psi3
state_dot[5] = -(1 - mu)/sigma3 * z - mu/psi3 * z
return state_dot

-------------------------------------------------------------------------------

t_0 = 0 # nondimensional time
t_f = 20 # nondimensional time (was 20)
t_points = np.linspace(t_0, t_f, 1000)
sol = solve_ivp(nondim_cr3bp, [t_0, t_f], Y_0, atol=1e-9, rtol=1e-6, t_eval=t_points)
Y = sol.y.T
r = Y[:, :3] # nondimensional distance
v = Y[:, 3:] # nondimensional velocity
x_2 = (1 - mu) * np.cos(np.linspace(0, np.pi, 100))
y_2 = (1 - mu) * np.sin(np.linspace(0, np.pi, 100))
fig, ax = plt.subplots(figsize=(5,5), dpi=96)
ax.plot(r[:, 0], r[:, 1], 'r', label="Trajectory") # Plot the orbits
ax.axhline(0, color='k')
ax.plot(np.hstack((x_2, x_2[::-1])), np.hstack((y_2, -y_2[::-1])))
ax.plot(-mu, 0, 'bo', label="$m_1$")
ax.plot(1 - mu, 0, 'go', label="$m_2$")
ax.plot(x_0, y_0, 'ro')
ax.set_aspect("equal")
plt.show()

Now Jacobi: orbital-mechanics.space/the-n-body-problem/jacobi-constant.html

speed_sq = np.sum(np.square(v), axis=1)
r[:, 0] += mu
sigma = np.sqrt(np.sum(np.square(r), axis=1))
r[:, 0] -= 1.0
psi = np.sqrt(np.sum(np.square(r), axis=1))
r[:, 0] = r[:, 0] + 1.0 - mu
J = 0.5 * speed_sq - (1 - mu) / sigma - mu / psi - 0.5 * ((1 - mu) * sigma2 + mu * psi2)
fig, ax = plt.subplots(figsize=(8,6), dpi=96)
ax.plot(sol.t, J, label="Jacobi Constant") # plot jacobi constant
ax.axhline(J[0], color="C1", label="Initial Jacobi Constant")
ax.legend(loc="center left")
ax.set_xlabel("$\tau$")
ax.set_ylabel("$J$")
plt.show()

[rslippert]

Purdue JNRHO 9:2 has a jacobi_constant=3.04719
To avoid confusion (I was confused by this), because this constant is sometimes a negative Jacobi constant
note that Purdue uses magnitude of the velocity, so see the below code for this calculation:

def jacobi_const( st ): # st is the state vector
global mu
r1 = math.sqrt( (st[0]+mu)**2 + st[1]**2 + st[2]**2 )
r2 = math.sqrt( (st[0]+mu-1)**2 + st[1]**2 + st[2]*2 )
j = 2((1-mu)/r1+mu/r2) + st[0]**2 + st[1]**2
- (st[3]**2 + st[4]**2 + st[5]**2)
return(j)