Two-Body Problem in the Co-Moving Frame#
Recall the equations of motion for the two-body problem in a reference frame attached to \(m_1\), Eq. (40), repeated here for reference:
In this example, we will solve the two-body problem in relative coordinates. Although the problem is in relative coordinates, the solution procedure is the same as Two-Body Numerical Solution in an Inertial Frame. In Python, the first thing is to import the required libraries. This step is not necessary in Matlab.
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D # noqa: F401, needed for 3d plots
from scipy.integrate import solve_ivp
A 1000 kg satellite is in orbit around Earth. The initial position and velocity of the satellite are:
We want to determine the minimum and maximum altitude, and the satellite’s speed at those two locations. The altitude of the satellite is its height above the earth’s surface. The radius of the earth is 6378.12 km.
To avoid confusion with the \(y\) coordinate, we will use capital \(Y\) for the state vector. Remember that the state vector includes 6 components, 3 positions and 3 velocities. These all need to be stored in one array.
G = 6.67430e-20 # km**3/(kg * s**2)
m_1 = 5.97219e24 # kg
m_2 = 1000 # kg
mu = G * m_1
R_E = 6378.12 # km
r_0 = np.array((8000, 0, 6000)) # km
v_0 = np.array((0, 7, 0)) # km/s
Y_0 = np.hstack((r_0, v_0))
function two_body_relative_numerical_solution
G = 6.67430e-20; % km^3/(kg * s^2)
m1 = 5.97219e24; % kg
m2 = 1000; % kg
mu = G * m1; % km^3/s^2
RE = 6378.12; % km
r0 = [8000 0 6000]; % km
v0 = [0 7 0]; % km/s
Y0 = [r0 v0];
Then, we need to define the function that describes the motion of our system. This function needs to compute the derivative of the state vector and return it to the solver.
Inside the function, we use the values in the state vector to fill the Ydot vector. Then, we return Ydot back to the solver.
def relative_motion(t, Y):
"""Calculate the motion of a two-body system relative to $m_1$.
The state vector ``Y`` should be in the order:
1. Coordinates of $m_2$ relative to $m_1$
2. Velocity components of $m_2$ relative to $m_1$
"""
# Get the three position components
r_vec = Y[:3]
# Create the derivative vector and copy the velocities into it
Ydot = np.zeros_like(Y)
Ydot[:3] = Y[3:]
# Calculate the accelerations
r = np.sqrt(np.sum(np.square(r_vec)))
a_vec = -mu * r_vec / r**3
Ydot[3:] = a_vec
return Ydot
function Ydot = relative_motion(t,Y)
rvector = Y(1:3);
vvector = Y(4:6);
r = norm(rvector);
avector = -mu .* rvector ./ r.^3;
Ydot = [vvector; avector];
end
With the function defined, we can call the solver function. We need to tell it the function it should solve, the beginning and end times, the initial state vector, and then some information to help control the output.
We will guess the final time and check that increasing it does not change the answers for minimum and maximum altitude.
In Python, once the solver finishes, the solution is stored in sol.y. Each column of sol.y corresponds to a single timestep and each row corresponds to one of the state variables. It is more convenient to work with the transpose of this array, so we do that and define Y.
t_0 = 0 # seconds
t_f = 14_709 # seconds, period of one orbit
t_points = np.linspace(t_0, t_f, 1000)
sol = solve_ivp(relative_motion, [t_0, t_f], Y_0, t_eval=t_points)
Y = sol.y.T
r = Y[:, :3] # km
v = Y[:, 3:] # km/s
t0 = 0; % seconds
tf = 14709; % seconds, period of one orbit
[t,Y] = ode45(@relative_motion, [t0, tf], Y0);
rvec = Y(:, 1:3);
vvec = Y(:, 4:6);
Since we are looking for the minimum and maximum altitude, we need to calculate the distance of the satellite from the center of the earth and then subtract the radius of the earth. This will give the altitude above the surface:
where \(h\) is the altitude and \(R_E\) is the radius of the earth. Similarly, the speed is the magnitude of the velocity.
r_mag = np.sqrt(np.sum(np.square(r), axis=1))
# altitude is the distance above the surface of the Earth
altitude = r_mag - R_E
speed = np.sqrt(np.sum(np.square(v), axis=1))
rmag = vecnorm(rvec,2,2);
altitude = rmag - RE;
speed = vecnorm(vvec,2,2);
Now we can process the arrays to find the minimum and maximum altitude, and the speed at those locations. In Python, there are separate functions to find the minimum/maximum and its index; in Matlab, there is a single function to do both.
min_altitude = np.min(altitude)
i_min = np.argmin(altitude)
max_altitude = np.max(altitude)
i_max = np.argmax(altitude)
speed_at_min_alt = speed[i_min]
speed_at_max_alt = speed[i_max]
time_at_min_alt = sol.t[i_min]
time_at_max_alt = sol.t[i_max]
[min_altitude, imin] = min(altitude);
[max_altitude, imax] = max(altitude);
speed_at_min_alt = speed(imin);
speed_at_max_alt = speed(imax);
time_at_min_alt = t(imin);
time_at_max_alt = t(imax);
Finally, we can print the results to the screen.
print(
f"""\
The minimum altitude during the orbit is: {min_altitude:.2F} km
The speed at the minimum altitude is: {speed_at_min_alt:.2F} km/s
The time at minimum altitude is: {time_at_min_alt:.2F} s
The maximum altitude during the orbit is: {max_altitude:.2F} km
The velocity at the maximum altitude is: {speed_at_max_alt:.4F} km/s
The time at maximum altitude is: {time_at_max_alt:.2F} s
"""
)
fprintf('The minimum altitude during the orbit is: %8.2f km\n', min_altitude)
fprintf('The speed at the minimum altitude is: %8.2f km/s\n', speed_at_min_alt)
fprintf('The time at minimum altitude is: %8.2f s\n', time_at_min_alt)
fprintf('The maximum altitude during the orbit is: %8.2f km\n', max_altitude)
fprintf('The speed at the maximum altitude is: %8.2f km/s\n', speed_at_max_alt)
fprintf('The time at maximum altitude is: %8.2f s\n', time_at_max_alt)
The minimum altitude during the orbit is: 3621.88 km
The speed at the minimum altitude is: 7.00 km/s
The time at minimum altitude is: 0.00 s
The maximum altitude during the orbit is: 9431.85 km
The velocity at the maximum altitude is: 4.4353 km/s
The time at maximum altitude is: 7288.57 s
Interestingly, the results from Python and Matlab are slightly different. This is most likely because of different time resolution from the solvers. Now let’s plot the orbit. The central sphere is representative of the earth.
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(r[:, 0], r[:, 1], r[:, 2], label="Orbit")
ax.set_xlabel("km")
ax.set_ylabel("km")
ax.set_zlabel("km")
# This adds a sphere to the plot of the radius of the Earth
p = np.linspace(0, np.pi, 200)
t = np.linspace(0, 2 * np.pi, 200)
P, T = np.meshgrid(p, t)
X = R_E * np.cos(T) * np.sin(P)
Y = R_E * np.sin(T) * np.sin(P)
Z = R_E * np.cos(P)
ax.plot_surface(X, Y, Z)
ax.plot(r[i_min, 0], r[i_min, 1], r[i_min, 2], "ro", label="Min. Altitude")
ax.plot(r[i_max, 0], r[i_max, 1], r[i_max, 2], "go", label="Max. Altitude")
ax.legend()
figure()
[xx, yy, zz] = sphere(200);
surf(RE*xx, RE*yy, RE*zz)
colormap(light_blue)
caxis([-RE/100 RE/100])
shading interp
% Draw and label the X, Y and Z axes
line([0 2*RE], [0 0], [0 0]);
text(2*RE + RE/10, 0, 0, 'x')
line([0 0], [0 2*RE], [0 0]);
text(0, 2*RE + RE/10, 0, 'y')
line([0 0], [0 0], [0 2*RE]);
text(0, 0, 2*RE + RE/10, 'z')
hold on
plot3(rvec(:,1), rvec(:,2), rvec(:,3), 'k')
plot3(rvec(imin,1), rvec(imin,2), rvec(imin,3), '.r')
text(rvec(imin,1), rvec(imin,2), rvec(imin,3), 'Min. Altitude')
plot3(rvec(imax,1), rvec(imax,2), rvec(imax,3), '.g')
text(rvec(imax,1), rvec(imax,2), rvec(imax,3), 'Max. Altitude')
view([1,1,.4])
% Specify some properties of the graph
grid on
axis equal
xlabel('km')
ylabel('km')
zlabel('km')
function map = light_blue
r = 0.8; g = r; b = 1.0;
map = [r g b
0 0 0
r g b];
end
Fig. 13 The motion of a satellite in orbit around the earth.#