Moon and Earth epicycloid

In this article, I will show you how to calculate the moon's orbit around the Sun and why this orbit is an epicycloid. The video with this story is published on YouTube.

Epicycloid

Despite it being discovered long ago, the trajectory of the Moon, not many students can quickly answer how long it takes for the Moon to make one rotation around the Earth, and what the Moon's trajectory is in the solar system. It is a big surprise that the Moon is rotating around the Sun, and the Earth only makes some tiny changes in its orbit, changing it from a pure ellipse to an epicycloidal ellipse.

But let's plot these trajectories with Python and visually check how they looks like

Ephemirides

Of course, it is possible to calculate the trajectories based on physics. By knowing the mass of the Sun, Earth and Moon, it is possible to model the trajectories, but we can use ephemerides with coordinates of all these spatial objects at any moment and just plot their coordinates or trajectories. DE421.bsp ephemeride will be good enough for these calculations.

Animated code

To easily calculate the positions of the Sun, Earth and Moon from de421 ephemeride, I create the following Python-based code:


# Import necessary libraries
import numpy as np                      # For numerical operations and arrays
import matplotlib.pyplot as plt         # For plotting and visualization
import matplotlib.animation as animation  # For making animations
from skyfield.api import load           # For precise astronomical calculations

# Load NASA JPL ephemerides (planetary positions)
# 'de421.bsp' contains solar system body data (Sun, Earth, Moon, planets...)
eph = load('de421.bsp')
sun, earth, moon = eph['sun'], eph['earth'], eph['moon']

# Create a time scale
ts = load.timescale()

# Generate times from Jan 1, 2025 onward:
#   np.arange(0, 365*24*60, 6*60) → every 6 hours for 1 year
#   divided by 60 because Skyfield expects hours as floats
t = ts.utc(2025, 1, 1, np.arange(0, 365*24*60, 6*60)/60)

# Get positions of Moon and Earth relative to the Sun
#   .at(t)        → position at each moment in array "t"
#   .observe(sun) → vector from body to Sun
#   .apparent()   → includes light-time correction
#   .ecliptic_position().km → position in ecliptic coordinates, in kilometers
moon_pos = moon.at(t).observe(sun).apparent().ecliptic_position().km
earth_pos = earth.at(t).observe(sun).apparent().ecliptic_position().km

# Sun is always at the origin in heliocentric coordinates
sun_pos = np.zeros((3, len(t)))

# Create a figure for animation
fig, ax = plt.subplots()
ax.set_aspect('equal')          # Keep x and y axis scale equal
ax.set_xlim(-2e8, 2e8)          # Show ±200 million km (enough for Earth's orbit)
ax.set_ylim(-2e8, 2e8)

# Create empty markers for Sun, Earth, and Moon
# 'yo' = yellow circle, 'bo' = blue circle, 'ro' = red circle
# markersize sets the visual size on screen (not real scale!)
line_sun, = ax.plot([], [], 'yo', markersize=120)
line_earth, = ax.plot([], [], 'bo', markersize=30)
line_moon, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation (start with nothing drawn)
def init():
    line_sun.set_data([], [])
    line_earth.set_data([], [])
    line_moon.set_data([], [])
    return line_sun, line_earth, line_moon
    
# Animation function: updates the positions frame by frame
def animate(i):
    # Place Earth at its current x,y
    line_earth.set_data([earth_pos[0,i]], [earth_pos[1,i]])
    # Place Moon at its current x,y
    line_moon.set_data([moon_pos[0,i]], [moon_pos[1,i]])
    return line_sun, line_earth, line_moon

# Create animation: run through all frames of "t"
ani = animation.FuncAnimation(fig, animate, frames=len(t), init_func=init, blit=True)

# Save the animation to an MP4 video
ani.save('moon.mp4', fps=60, dpi=150)

# Close the figure (prevents it from showing interactively)
plt.close(fig)

# --- Static plot of a selected time interval ---

# Define time window (June 1 – July 12, 2025)
start = ts.utc(2025, 6, 1)
end   = ts.utc(2025, 7, 12)
mask = (t.tt >= start.tt) & (t.tt <= end.tt)  # Boolean mask to slice arrays

# Slice Moon and Earth positions for this time interval
moon_cut = moon_pos[:, mask]
earth_cut = earth_pos[:, mask]

# Plot Earth and Moon paths (trajectories) in that window
plt.figure(figsize=(8,8))
plt.plot(earth_cut[0], earth_cut[1], 'b-', label="Earth")  # Earth's orbit
plt.plot(moon_cut[0], moon_cut[1], 'r-', label="Moon")     # Moon's orbit

plt.legend()
plt.show()  # Display the static plot

The program starts by importing the tools it needs: numerical calculations with NumPy, plotting and animation with Matplotlib, and astronomical data from Skyfield. It then loads NASA's DE421 ephemerides, a precise dataset describing the positions of the Sun, Earth, Moon, and other solar system bodies.

Next, a time sequence is defined: every six hours throughout the entire year of 2025. At each of these times, the code asks Skyfield for the position of the Earth and the Moon as seen from the Sun, expressed in kilometres in ecliptic coordinates. The Sun itself is kept fixed at the origin, since we are working in a heliocentric frame.

With this data prepared, the program builds a Matplotlib animation. A figure is created where the Earth, Moon, and Sun are represented as markers of different colours and sizes. An initialisation function clears the plot, and an animation function updates the Earth and Moon's position at each step in time. The Sun stays fixed at the centre. This animation is then saved as a video file, showing the Earth circling the Sun and the Moon tracing its path around the Earth, producing the characteristic epicyclic motion.

After the animation is finished, the code produces a second visualisation. This time, instead of moving markers, it shows the trails left by the Earth and Moon over a shorter time span - from June 1st to July 12th, 2025. By applying a time filter, the program extracts just that portion of the data and plots the Earth's and Moon's orbits as continuous blue and red lines. This static figure makes it easier to see the shape of the Moon's path compared to Earth's orbit around the Sun.

In summary, the script first generates a dynamic animation of the Earth–Moon–Sun system across a full year, and then a static diagram of their orbital trajectories across a chosen two-month period.

Published: 2025-09-07 09:09:11
Updated: 2025-09-07 09:10:30