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)/sigma*3 z - mu/psi*3 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) * sigma*2 + mu psi**2) 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)