Two Body Problem

Typically, a two-body problem entails finding a solution to the motion of two masses that are affected by one another's gravitational pull. In orbital mechanics, the smaller of the two masses is often thought to have insignificant mass, having little impact on the larger mass. It is important to understand and solve for the motion of the insignificant mass. The smaller mass does influence the larger mass, but the impact is so negligible that it is effectively nil.

The equation of motion for a negligible mass under the influence of gravity of a larger mass are as follows:

# Importing Packages
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Earth Model
def model_2BP(state, t):
    mu = 3.986004418E+05  # Earth's gravitational parameter  
                          # [km^3/s^2]
    x = state[0]
    y = state[1]
    z = state[2]
    x_dot = state[3]
    y_dot = state[4]
    z_dot = state[5]
    x_ddot = -mu * x / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
    y_ddot = -mu * y / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
    z_ddot = -mu * z / (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)
    dstate_dt = [x_dot, y_dot, z_dot, x_ddot, y_ddot, z_ddot]
    return dstate_dt

# Initial Conditions
X_0 = -2500  # [km]
Y_0 = -5500  # [km]
Z_0 = 3400  # [km]
VX_0 = 7.5  # [km/s]
VY_0 = 0.0  # [km/s]
VZ_0 = 4.0  # [km/s]
state_0 = [X_0, Y_0, Z_0, VX_0, VY_0, VZ_0]

# Time Array
t = np.linspace(0, 6*3600, 200)  # Simulates for a time period of 6
                                 # hours [s]

# Solving ODE
sol = odeint(model_2BP, state_0, t)
X_Sat = sol[:, 0]  # X-coord [km] of satellite over time interval 
Y_Sat = sol[:, 1]  # Y-coord [km] of satellite over time interval
Z_Sat = sol[:, 2]  # Z-coord [km] of satellite over time interval

# Setting up Spherical Earth to Plot
N = 50
phi = np.linspace(0, 2 * np.pi, N)
theta = np.linspace(0, np.pi, N)
theta, phi = np.meshgrid(theta, phi)

r_Earth = 6378.14  # Average radius of Earth [km]
X_Earth = r_Earth * np.cos(phi) * np.sin(theta)
Y_Earth = r_Earth * np.sin(phi) * np.sin(theta)
Z_Earth = r_Earth * np.cos(theta)

# Plotting Earth and Orbit
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_surface(X_Earth, Y_Earth, Z_Earth, color='blue', alpha=0.7)
ax.plot3D(X_Sat, Y_Sat, Z_Sat, 'black')
ax.view_init(30, 145)  # Changing viewing angle (adjust as needed)
plt.title('Two-Body Orbit')
ax.set_xlabel('X [km]')
ax.set_ylabel('Y [km]')
ax.set_zlabel('Z [km]')

The code provided simulates a two-body orbit of a satellite around the Earth using numerical integration. Here's a breakdown of the code:

Importing Packages: The necessary packages (NumPy, SciPy, and Matplotlib) are imported for numerical computations and plotting.
Earth Model: The function model_2BP represents the differential equations governing the motion of the satellite in a two-body problem. It calculates the acceleration components (x_ddot, y_ddot, z_ddot) based on the satellite's position (x, y, z) and the Earth's gravitational parameter (mu).
Initial Conditions: The initial conditions (X_0, Y_0, Z_0, VX_0, VY_0, VZ_0) for the satellite's position and velocity are defined.
Time Array: The time array (t) is created to define the time span for the simulation.
Solving ODE: The odeint function from SciPy is used to solve the system of differential equations numerically (model_2BP) using the initial conditions and time array. The resulting solution (sol) contains the satellite's position coordinates (X_Sat, Y_Sat, Z_Sat) over the specified time interval.
Setting up Spherical Earth: The coordinates (X_Earth, Y_Earth, Z_Earth) for the spherical shape of the Earth are calculated based on the average Earth radius (r_Earth) and the angles (phi and theta).
Plotting Earth and Orbit: Matplotlib's plot_surface function is used to plot the Earth as a blue surface, and the satellite's orbit is plotted as a black line using the plot3D function. The viewing angle of the plot is adjusted using ax.view_init, and labels for the axes are set.