Astrophysics & AI with Python: Simulating Planetary Orbits with Kepler's Laws

Ever looked up at the night sky and wondered how we know the precise dance of the planets? It’s not magic—it’s math. Specifically, it’s the elegant interplay between Johannes Kepler’s observations and Isaac Newton’s laws of gravity.

But how do we take those centuries-old equations and turn them into a living, breathing simulation on a computer screen?

In this chapter of our journey through computational physics, we are moving from theory to practice. We will bridge the gap between abstract physics and concrete code by building a robust orbital simulator. We'll explore the stability of numerical methods and write Python code that can predict the path of a planet around the Sun with frightening accuracy.

To simulate an orbit, we first need to understand the two pillars of celestial mechanics that make it possible.

Johannes Kepler spent years staring at Tycho Brahe’s astronomical data until three distinct patterns emerged. These are the "rules of the road" for the cosmos:

Kepler told us what happens; Newton told us why. His Law of Universal Gravitation provides the engine for our simulation:

By combining this with Newton's Second Law ( F=ma ), we get the acceleration vector that drives the motion:

This equation is the heart of our code. However, because physics happens in continuous time and computers operate in discrete steps, we have to solve a differential equation numerically. This brings us to the critical challenge of the Two-Body Problem.

Imagine driving in a fog. You know your current speed and direction, so you guess where you'll be in one minute. Then you check your new position and guess again. This is numerical integration.

The simplest way to do this is the Standard Euler Method:

The Problem: This method is unstable. Over time, it fails to conserve energy. The planet will slowly spiral out into space or crash into the Sun, violating the laws of physics.

The Solution: We will use the Euler-Cromer Method. It’s a tiny tweak with a massive impact. Instead of using the old velocity to update the position, we use the newly calculated velocity:

This "semi-implicit" approach is symplectic, meaning it keeps the energy error bounded. The orbit remains stable for millions of simulated years.

Here is the complete Python script to simulate Earth's orbit around the Sun. We use numpy

for efficient vector math and matplotlib

for visualization.

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt


G = 6.67430e-11
M_sun = 1.989e30

R_AU = 1.496e11 # 1 Astronomical Unit in meters
r0 = np.array([R_AU, 0.0])

V_EARTH = 2.978e4 # Earth's average orbital speed
v0 = np.array([0.0, V_EARTH])


TOTAL_TIME = 365.25 * 24 * 3600
DT = 3600.0 * 6
NUM_STEPS = int(TOTAL_TIME / DT)


positions = np.zeros((NUM_STEPS, 2))
positions[0] = r0

r_current = r0.copy()
v_current = v0.copy()


def calculate_acceleration(r_vec, M_central):
    """
    Calculates the acceleration vector (a) due to gravity.
    a = - (G * M / r^2) * r_hat
    """
    r_mag = np.linalg.norm(r_vec)

    if r_mag == 0:
        return np.zeros(2)

    r_hat = r_vec / r_mag

    a_mag = - (G * M_central) / (r_mag**2)

    return a_mag * r_hat


print(f"Starting simulation: {NUM_STEPS} steps over {TOTAL_TIME/3600/365.25:.2f} years.")

for i in range(1, NUM_STEPS):
    a = calculate_acceleration(r_current, M_sun)

    v_current = v_current + a * DT

    r_current = r_current + v_current * DT

    positions[i] = r_current

print("Simulation complete.")


x_au = positions[:, 0] / R_AU
y_au = positions[:, 1] / R_AU

plt.figure(figsize=(8, 8))
plt.plot(x_au, y_au, label='Simulated Orbit')
plt.plot(0, 0, 'o', color='gold', markersize=10, label='Sun (Origin)')

plt.title(f'Orbital Simulation using Euler-Cromer Integration ({NUM_STEPS} steps)')
plt.xlabel('X Position (AU)')
plt.ylabel('Y Position (AU)')
plt.gca().set_aspect('equal', adjustable='box')
plt.grid(True, linestyle=':', alpha=0.7)
plt.legend()
plt.show()

If you are new to numerical physics, the code might look intimidating. Let's break down the logic flow:

np.linalg.norm

is the key function here—it calculates the distance r between the planet and the Sun.plt.gca().set_aspect('equal')

is vital; without it, a perfect circular orbit would look like an ellipse because the x and y scales would be different.By implementing the Euler-Cromer method, we have successfully translated the laws of Kepler and Newton into a computational algorithm. We didn't just plot a circle; we simulated the dynamic interplay of forces that keeps the planets in orbit.

This code is more than just a visualization tool—it is a foundation. The principles used here (calculating acceleration, updating state variables, and verifying conservation laws) are the exact same steps used in advanced astrophysics to simulate galaxy collisions, satellite trajectories, and the formation of planetary systems.

DT

to be very large (e.g., 1 day instead of 6 hours), what physical phenomena would you expect to see in the resulting orbit plot?The concepts and code demonstrated here are drawn directly from the comprehensive roadmap laid out in the ebook

Astrophysics & AI: Building Research Agents for Astronomy, Cosmology, and SETI. You can find it here. Check all the other 50 Programming & AI ebooks with python, typescript, swift, c#: here