Simulation of the Earth-Moon orbit

Papadopol Lucian Ioan
l.i.papadopol@gmail.com

1. Objective of the simulation

In this exercise we want to simulate, in a simplified way, the orbital motion of the Moon around the Earth.

This motion:

1. is governed by the gravitational force;
2. is a dynamical system in time;
3. in the ideal case, the orbit should remain stable;
4. mechanical energy should be conserved.

2. Simplified physical model

We consider the Earth-Moon problem in the simplest possible form.

We make the following assumptions:

1. the Earth is fixed at the origin of the reference frame;
2. the Moon is treated as a point particle;
3. the motion takes place in a two-dimensional plane;
4. the only force considered is the Earth’s gravity;
5. we neglect the Sun;
6. we neglect the other planets;
7. we neglect everything else;
8. we assume that the Earth’s mass is much greater than the Moon’s mass.

In this way we obtain a simplified Kepler problem: a small body orbiting around a much more massive central body.

3. Reference frame

We place the Earth fixed at the origin of the Cartesian reference frame:

The position of the Moon is described by the vector:

where:

• x(t) is the horizontal coordinate of the Moon;
• y(t) is the vertical coordinate of the Moon;
• t is time.

The velocity of the Moon is:

4. Earth-Moon distance

The distance of the Moon from the Earth is the magnitude of the position vector:

This distance changes over time if the orbit is not perfectly circular.

In the ideal case of a circular orbit, instead, the distance remains approximately constant.

5. Gravitational force

The gravitational force exerted by the Earth on the Moon is described by Newton’s law of universal gravitation.

In magnitude:

where:

• G is the gravitational constant;
• M_T is the mass of the Earth;
• m_L is the mass of the Moon;
• r is the distance between the Earth and the Moon.

The force is always directed toward the Earth, therefore toward the origin of the reference frame.

6. Gravitational acceleration

To simulate the motion of the Moon, it is not necessary to compute the force directly. We can work directly with the acceleration.

From Newton’s second law:

and from the gravitational force we obtain:

We introduce the Earth’s gravitational parameter:

Then the acceleration becomes:

In components:

Given the Moon’s position vector, we can compute the acceleration due to the Earth’s gravity.

7. Differential system

The complete state of the system is formed by position and velocity:

The associated differential system is:

with:

This system describes how the position and velocity of the Moon change over time.

8. Position-velocity form

For the Leap-frog method it is better to update position and velocity separately.

where:

and:

Leap-frog updates position and velocity in a staggered way.

The idea is:

1. the velocity is updated at a half step;
2. the position is updated using that intermediate velocity;
3. the velocity is then completed using the acceleration at the new position.

Therefore, at each time step the method computes the gravitational acceleration from the current position of the Moon.

9. Initial conditions

To obtain an almost circular orbit, we can initialize the Moon at its average distance from the Earth:

Initially, we place the Moon on the x-axis:

The initial velocity is chosen perpendicular to the radius, therefore along the y-axis:

For an ideal circular orbit, the orbital speed is:

This choice produces an approximately circular orbit in the simplified model.

10. Orbital mechanical energy

In the ideal gravitational problem, without friction and without external perturbations, the total mechanical energy should be conserved.

The specific energy, namely the energy per unit mass, is:

where:

The first term is the specific kinetic energy:

The second term is the specific gravitational potential energy:

If the numerical method is good, the energy should remain approximately constant during the simulation.

11. Angular momentum

In the two-dimensional case, the specific angular momentum is:

The angular momentum should be conserved. This means that the orbit should not be artificially deformed over time.

First step – Python function representing the differential system

def model(t, u):
  x, y, vx, vy = u

  r = np.sqrt(x**2 + y**2)

  dxdt = vx
  dydt = vy
  dvxdt = -mu/r**3 * x
  dvydt = -mu/r**3 * y

  return [dxdt, dydt, dvxdt, dvydt]

Second step – Leap-frog method implemented in class

import numpy as np

def leapfrog_step(u_prev, u_curr, model, t_curr, dt):
  return np.asarray(u_prev, dtype=float) + 2.0 * dt * np.asarray(model(t_curr, u_curr), dtype=float)

def integrate_leapfrog(u0, model, t0, t_end, dt):
    """Time integration loop with the Leap-frog method"""

    n = int(np.ceil((t_end - t0) / dt))

    t = t0 + dt * np.arange(n + 1)

    u = np.zeros((n + 1,) + np.shape(u0), dtype=float)

    # First state
    u[0] = u0

    # Second state: needed to start Leap-frog
    # It is estimated with one Euler step
    u[1] = u[0] + dt * np.asarray(model(t[0], u[0]), dtype=float)

    # From here onward, the true Leap-frog scheme is used
    for k in range(1, n):
        u[k + 1] = leapfrog_step(u[k - 1], u[k], model, t[k], dt)

    return t, u

Third step – setting the initial conditions

# Physical constants

G = 6.67430e-11        # gravitational constant [m^3 kg^-1 s^-2]
M_earth = 5.972e24     # Earth mass [kg]

mu = G * M_earth       # Earth's gravitational parameter [m^3/s^2]

# Average Earth-Moon distance

r0 = 384400e3          # [m]

# Ideal circular orbital speed

v0 = np.sqrt(mu / r0)  # [m/s]

# Initial position of the Moon

x0 = r0
y0 = 0.0

# Initial velocity of the Moon
# For an ideal circular orbit it must be perpendicular to the radius

vx0 = 0.0
vy0 = v0

# Initial state: u = [x, y, vx, vy]

u0 = np.array([x0, y0, vx0, vy0], dtype=float)

# Time parameters

t0 = 0.0

day = 24 * 3600
t_end = 30 * day       # about 30 days of simulation

dt = 3600.0            # time step: 1 hour

Fourth step – running the simulation and visualizing the result

t, u = integrate_leapfrog(u0, model, t0, t_end, dt)

x = u[:, 0]
y = u[:, 1]
vx = u[:, 2]
vy = u[:, 3]

import matplotlib.pyplot as plt

plt.figure(figsize=(7, 7))

plt.plot(x / 1000, y / 1000, label="Moon orbit")
plt.scatter(0, 0, s=200, label="Earth")
plt.scatter(x[0] / 1000, y[0] / 1000, s=80, label="Initial Moon position")

plt.xlabel("x [km]")
plt.ylabel("y [km]")
plt.title("Earth-Moon orbit simulation with the Leap-frog method")
plt.axis("equal")
plt.grid(True)
plt.legend()
plt.show()

Fifth step – implementing Euler and RK2 to compare the result

def euler_step(u, model, t, dt):
    """Single explicit Euler step"""
    u = np.asarray(u, dtype=float)

    return u + dt * np.asarray(model(t, u), dtype=float)


def integrate_euler(u0, model, t0, t_end, dt):
    """Time integration loop with explicit Euler"""

    n = int(np.ceil((t_end - t0) / dt))

    t = t0 + dt * np.arange(n + 1)

    u = np.zeros((n + 1,) + np.shape(u0), dtype=float)
    u[0] = u0

    for k in range(n):
        u[k + 1] = euler_step(u[k], model, t[k], dt)

    return t, u

def rk2_midpoint_step(u, model, t, dt):
  """Second-order Runge-Kutta method with midpoint"""
  u = np.asarray(u, dtype=float)
  k1 = np.asarray(model(t,u), dtype=float)
  u_mid = u + dt/2*k1
  t_mid = t + dt/2
  k2 = np.asarray(model(t_mid, u_mid), dtype=float)
  return u + dt*k2

def integrate_rk2(u0, model, t0, t_end, dt):
  """Time integration loop"""
  n = int(np.ceil((t_end -t0)/dt))
  t = t0 + dt * np.arange(n+1)
  u = np.zeros((n+1,)+np.shape(u0), dtype=float)
  u[0] = u0
  for k in range(n):
    u[k+1]=rk2_midpoint_step(u[k], model, t[k], dt)
  return t,u

Sixth step – checking the difference between the three methods

t_euler, u_euler = integrate_euler(u0, model, t0, t_end, dt)

t_rk2, u_rk2 = integrate_rk2(u0, model, t0, t_end, dt)

t_lf, u_lf = integrate_leapfrog(u0, model, t0, t_end, dt)

x_euler = u_euler[:, 0]
y_euler = u_euler[:, 1]

x_rk2 = u_rk2[:, 0]
y_rk2 = u_rk2[:, 1]

x_lf = u_lf[:, 0]
y_lf = u_lf[:, 1]

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 8))

plt.plot(
    x_euler / 1000, y_euler / 1000,
    label="Explicit Euler",
    linewidth=2
)

plt.plot(
    x_rk2 / 1000, y_rk2 / 1000,
    label="RK2 midpoint",
    linewidth=2,
    linestyle="--"
)

plt.plot(
    x_lf / 1000, y_lf / 1000,
    label="Leap-frog",
    linewidth=2,
    linestyle=":"
)

plt.scatter(0, 0, s=220, label="Earth")
plt.scatter(u0[0] / 1000, u0[1] / 1000, s=80, label="Initial Moon position")

plt.xlabel("x [km]")
plt.ylabel("y [km]")
plt.title("Earth-Moon orbit comparison: Euler, RK2 and Leap-frog")

plt.axis("equal")
plt.grid(True)
plt.legend()
plt.show()

Checking energy conservation with the three methods

Energy function:

def orbital_energy(u, mu):
    """
    Specific orbital mechanical energy.

    State:
    u = [x, y, vx, vy]

    E = 1/2 * v^2 - mu/r
    """

    x = u[:, 0]
    y = u[:, 1]
    vx = u[:, 2]
    vy = u[:, 3]

    r = np.sqrt(x**2 + y**2)
    v2 = vx**2 + vy**2

    E = 0.5 * v2 - mu / r

    return E

Energy calculation with the three methods:

E_euler = orbital_energy(u_euler, mu)
E_rk2 = orbital_energy(u_rk2, mu)
E_lf = orbital_energy(u_lf, mu)

Relative error on the energy:

err_E_euler = (E_euler - E_euler[0]) / abs(E_euler[0])
err_E_rk2 = (E_rk2 - E_rk2[0]) / abs(E_rk2[0])
err_E_lf = (E_lf - E_lf[0]) / abs(E_lf[0])

plt.figure(figsize=(9, 5))

plt.plot(
    t_euler / day,
    E_euler / E_euler[0],
    color="darkred",
    linewidth=2,
    label="Explicit Euler"
)

plt.plot(
    t_rk2 / day,
    E_rk2 / E_rk2[0],
    color="orange",
    linewidth=2,
    linestyle="--",
    label="RK2 midpoint"
)

plt.plot(
    t_lf / day,
    E_lf / E_lf[0],
    color="blue",
    linewidth=2,
    linestyle="-.",
    label="Leap-frog"
)

plt.xlabel("Time [days]")
plt.ylabel("Normalized energy E(t)/E(0)")
plt.title("Normalized orbital energy")
plt.grid(True)
plt.legend()
plt.show()