# Creating Three-Body Orbits Using Python: A Practical Guide

Written on

## Understanding the Three-Body Problem

In the realm of orbital mechanics, the three-body problem (3BP) examines the motion of a negligible mass influenced by two larger masses, such as a planet, moon, or star. This analysis is crucial for planning spacecraft trajectories where the gravitational effects of two significant bodies must be accounted for. For instance, the trajectory of the James Webb Space Telescope was designed with considerations of the 3BP in mind.

While many people associate orbits with a small mass following an elliptical path around a single celestial body, the introduction of an additional massive body complicates the derivation of motion equations for the smaller mass.

To simulate and visualize these orbits in the context of the 3BP, we need to derive the equations of motion for the negligible mass. Although I won’t delve into the derivation process here, I recommend checking out this article that provides a comprehensive explanation of the equations relevant to the three-body problem. Understanding the variables and equations you will be solving with Python is essential.

We will focus on the circular-restricted three-body problem (CR3BP). In this version, the two larger masses, referred to as the primaries, orbit their common center of mass in circular paths. This simplifies the equations governing the motion of a spacecraft or other negligible masses.

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

## Conclusion

Thank you for reading this overview of the circular-restricted three-body problem and the numerical integration of its motion equations. This topic is complex within the field of orbital mechanics, making it a challenge to master. Should you have any questions, feel free to reach out via LinkedIn or here. Don’t forget to follow for weekly insights on orbital mechanics, coding, and machine learning!

If you found this article intriguing, consider exploring additional resources that may deepen your understanding of these concepts.