The Motion Equations for CR3BP

The equations of motion for the negligible mass (denoted as m₃) in the rotating frame (x, y, z-hat frame) are as follows:

These equations may seem daunting, but they become manageable once you understand the meaning of each variable. Here, x, y, and z (ρ vector) represent the position coordinates of m₃ relative to the center of mass of the primaries, in a rotating frame that moves with the primaries. The dots over these variables indicate time derivatives: one dot represents velocity, while two dots indicate acceleration. The parameter µ denotes the non-dimensional mass ratio of the primaries, while r₁ and r₂ are the distances from m₃ to each primary.

The aforementioned equations employ non-dimensional variables. The article referenced earlier elaborates on the use of these non-dimensional parameters, which is vital for grasping the intricacies of the three-body problem. For convenience, we use M* for mass units (kg), L* for length units (km), and T* for time units (s). More details will be provided in the coding section.

To solve the CR3BP, we will utilize an ordinary differential equation (ODE) solver to numerically integrate the motion equations for m₃. This requires defining a state vector, which includes the position and velocity vectors in the rotating frame, along with a time-derivative state vector.

The first video, "3 Body Simulation in Python: Rare Orbits," offers a visual demonstration of the three-body problem and its solutions in Python.

Setting Up the Code

To implement our code, we need to import several packages:

• NumPy: For creating and manipulating arrays.
• odeint from SciPy: For numerical integration.
• pyplot from Matplotlib: For visualizing the numerical results.

# Importing Packages

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

Creating the Numerical Integration Function

Next, we will define a Python function called model_CR3BP, which will be utilized by odeint to carry out the numerical integration of our state vector. This function will extract the position and velocity components from the state vector to compute the time-derivative state vector.

# CR3BP Model

def model_CR3BP(state, t):

x = state[0]

y = state[1]

z = state[2]

x_dot = state[3]

y_dot = state[4]

z_dot = state[5]

x_ddot = x + 2*y_dot - ((1 - mu)*(x + mu)) / ((x + mu)**2 + y**2 + z**2)**(3/2)

• (mu*(x - (1 - mu))) / ((x - (1 - mu))**2 + y**2 + z**2)**(3/2)

y_ddot = y - 2*x_dot - ((1 - mu)*y) / ((x + mu)**2 + y**2 + z**2)**(3/2)

• (mu*y) / ((x - (1 - mu))**2 + y**2 + z**2)**(3/2)

z_ddot = -((1 - mu)*z) / ((x + mu)**2 + y**2 + z**2)**(3/2)

• (mu*z) / ((x - (1 - mu))**2 + y**2 + z**2)**(3/2)

dstate_dt = [x_dot, y_dot, z_dot, x_ddot, y_ddot, z_ddot]

return dstate_dt

Defining Non-Dimensional Parameters

For our simulation, we will focus on the Earth-Moon system. However, you can adapt the parameters for different systems, such as Pluto-Charon or Sun-Jupiter. Below are the relevant parameters:

# Defining ND Parameters

G = 6.67408E-20 # Universal Gravitational Constant [km^3 kg^-1 s^-2]

mEarth = 5.97219E+24 # Mass of the Earth [kg]

mMoon = 7.34767E+22 # Mass of the Moon [kg]

a = 3.844E+5 # Semi-major axis of Earth and Moon [km]

m1 = mEarth

m2 = mMoon

Mstar = m1 + m2 # Non-Dimensional Mass Parameter

Lstar = a # Non-Dimensional Length Parameter

Tstar = (Lstar**3 / (G * Mstar))**(1/2) # Non-Dimensional Time Parameter

mu = m2 / Mstar

print('u03BC = ' + str(mu))

Initial Conditions and Time Interval

Before we can run the simulation, we need to set initial conditions and define a time interval for the integration. The following example shows arbitrary initial conditions that will generate the orbits seen earlier.

# Initial Conditions [Initial State Vector]

X_0 = 50000 / Lstar # Non-Dimensional x

Y_0 = 0 # Non-Dimensional y

Z_0 = 0 # Non-Dimensional z

VX_0 = 1.08 * Tstar / Lstar # Non-Dimensional x_dot

VY_0 = 3.18 * Tstar / Lstar # Non-Dimensional y_dot

VZ_0 = 0.68 * Tstar / Lstar # Non-Dimensional z_dot

state_0 = [X_0, Y_0, Z_0, VX_0, VY_0, VZ_0] # Non-Dimensional Initial Conditions

# Time Array

t = np.linspace(0, 15, 1000) # Non-Dimensional Time

Numerical Integration of the Model

Now we will combine our previous code into the odeint function to perform the numerical integration, yielding a state vector for each time step defined in our array.

# Numerically Integrating

sol = odeint(model_CR3BP, state_0, t)

# Rotational Frame Position Time History

X_rot = sol[:, 0]

Y_rot = sol[:, 1]

Z_rot = sol[:, 2]

# Inertial Frame Position Time History

X_Iner = sol[:, 0] * np.cos(t) - sol[:, 1] * np.sin(t)

Y_Iner = sol[:, 0] * np.sin(t) + sol[:, 1] * np.cos(t)

Z_Iner = sol[:, 2]

Incorporating Primary Masses

While optional, it can be useful to track the motion of the two primaries to better understand m₃'s trajectory. The following code calculates their positions in both the rotating and inertial frames.

# Constant m1 and m2 Rotational Frame Locations for CR3BP Primaries

m1_loc = [-mu, 0, 0]

m2_loc = [(1 - mu), 0, 0]

# Moving m1 and m2 Inertial Locations for CR3BP Primaries

X_m1 = m1_loc[0] * np.cos(t) - m1_loc[1] * np.sin(t)

Y_m1 = m1_loc[0] * np.sin(t) + m1_loc[1] * np.cos(t)

Z_m1 = m1_loc[2] * np.ones(len(t))

X_m2 = m2_loc[0] * np.cos(t) - m2_loc[1] * np.sin(t)

Y_m2 = m2_loc[0] * np.sin(t) + m2_loc[1] * np.cos(t)

Z_m2 = m2_loc[2] * np.ones(len(t))

Visualizing the Results

Finally, we will create plots to visualize the motion of each mass in both the rotating and inertial frames. The colors used will help distinguish between the spacecraft (green), the Moon (black), and the Earth (blue).

# Rotating Frame Plot

fig = plt.figure()

ax = plt.axes(projection='3d')

# Adding Figure Title and Labels

ax.set_title('Rotating Frame CR3BP Orbit (u03BC = ' + str(round(mu, 6)) + ')')

ax.set_xlabel('x [ND]')

ax.set_ylabel('y [ND]')

ax.set_zlabel('z [ND]')

# Plotting Rotating Frame Positions

ax.plot3D(X_rot, Y_rot, Z_rot, c='green')

ax.plot3D(m1_loc[0], m1_loc[1], m1_loc[2], c='blue', marker='o')

ax.plot3D(m2_loc[0], m2_loc[1], m2_loc[2], c='black', marker='o')

# Setting Axis Limits

xyzlim = np.array([ax.get_xlim3d(), ax.get_ylim3d(), ax.get_zlim3d()]).T

XYZlim = np.asarray([min(xyzlim[0]), max(xyzlim[1])])

ax.set_xlim3d(XYZlim)

ax.set_ylim3d(XYZlim)

ax.set_zlim3d(XYZlim * 3 / 4)

# Inertial Frame Plot

fig = plt.figure()

ax = plt.axes(projection='3d')

# Adding Figure Title and Labels

ax.set_title('Inertial Frame CR3BP Orbit (u03BC = ' + str(round(mu, 6)) + ')')

ax.set_xlabel('X [ND]')

ax.set_ylabel('Y [ND]')

ax.set_zlabel('Z [ND]')

# Plotting Inertial Frame Positions

ax.plot3D(X_Iner, Y_Iner, Z_Iner, c='green')

ax.plot3D(X_m1, Y_m1, Z_m1, c='blue')

ax.plot3D(X_m2, Y_m2, Z_m2, c='black')

# Setting Axis Limits

ax.set_xlim3d(XYZlim)

ax.set_ylim3d(XYZlim)

ax.set_zlim3d(XYZlim * 3 / 4)

plt.show()

The output displays the orbits of the spacecraft and the primary masses in both frames of reference. As illustrated, the spacecraft's orbit can be chaotic with certain initial conditions, leading to various potential orbital paths, including periodic ones. I encourage experimentation with different initial conditions to discover unique trajectories.

The second video, "Three Body Problem In Python - Earth's Orbit," provides further insights into simulating the three-body problem with a focus on Earth's orbital mechanics.