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.

The Physics: From Kepler’s Curves to Newton’s Force

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

1. The Phenomenology: Kepler’s Laws

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:

The Law of Ellipses: Orbits aren't perfect circles; they are ellipses with the Sun at one focus. This dictates the shape of our simulation.
The Law of Equal Areas: A planet speeds up when it's close to the Sun and slows down when it's far away. This dictates the dynamics of the simulation.
The Harmonic Law: The relationship between the size of the orbit ( aa ) and the time it takes to complete it ( TT ) is fixed ( T2∝a3T^2 \propto a^3 ). This is our ultimate verification metric.

2. The Causality: Newton’s Gravity

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

F=Gm1m2r2F = G \frac{m_1 m_2}{r^2}

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

a=−GMr2r^\mathbf{a} = -G \frac{M}{r^2} \hat{\mathbf{r}}

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.

The Computational Challenge: Why Stability Matters

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:

Update Velocity: vnew=vold+aoldΔt\mathbf{v}{new} = \mathbf{v}{old} + \mathbf{a}_{old} \Delta t
Update Position: rnew=rold+voldΔt\mathbf{r}{new} = \mathbf{r}{old} + \mathbf{v}_{old} \Delta t

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:

Update Velocity: vnew=vold+aoldΔt\mathbf{v}{new} = \mathbf{v}{old} + \mathbf{a}_{old} \Delta t
Update Position: rnew=rold+vnewΔt\mathbf{r}{new} = \mathbf{r}{old} + \mathbf{v}_{new} \Delta t (Note the change here!)

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

The Python Simulation Code

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

# --- 1. PHYSICAL CONSTANTS AND INITIAL SETUP ---

# Gravitational Constant G (SI units: m^3 kg^-1 s^-2)
G = 6.67430e-11
# Mass of the Central Body (Sun) (SI units: kg)
M_sun = 1.989e30

# Initial Conditions (Approximating Earth's orbit at perihelion)
# We use NumPy arrays to represent 2D vectors (x, y)
# Position vector r0 (meters). Placed on the positive x-axis. (~1 AU)
R_AU = 1.496e11 # 1 Astronomical Unit in meters
r0 = np.array([R_AU, 0.0])

# Velocity vector v0 (meters/second). Initial velocity is purely in the y-direction.
V_EARTH = 2.978e4 # Earth's average orbital speed
v0 = np.array([0.0, V_EARTH])

# --- 2. SIMULATION PARAMETERS ---

# Total duration of the simulation (1 Earth year in seconds)
TOTAL_TIME = 365.25 * 24 * 3600
# Time step (dt) for integration (6 hours in seconds). Smaller dt means higher accuracy.
DT = 3600.0 * 6
# Calculate the total number of discrete steps
NUM_STEPS = int(TOTAL_TIME / DT)

# --- 3. DATA STORAGE AND INITIALIZATION ---

# Array to store the X and Y positions at every step.
# Dimensions: (Number of steps, 2 dimensions)
positions = np.zeros((NUM_STEPS, 2))
# Store the initial position in the first row
positions[0] = r0

# Initialize the current state variables for the loop
r_current = r0.copy()
v_current = v0.copy()

# --- 4. ACCELERATION FUNCTION ---

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

    # Check for division by zero (e.g., if the planet hits the Sun)
    if r_mag == 0:
        return np.zeros(2)

    # Calculate the unit vector pointing from the central body to the planet
    r_hat = r_vec / r_mag

    # Calculate the magnitude of the gravitational acceleration
    # The negative sign ensures the acceleration vector points toward the origin (attraction)
    a_mag = - (G * M_central) / (r_mag**2)

    # The acceleration vector is the magnitude times the unit vector
    return a_mag * r_hat

# --- 5. MAIN INTEGRATION LOOP (Euler-Cromer Method) ---

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

for i in range(1, NUM_STEPS):
    # 5a. Calculate acceleration (a) at the current position (r_current)
    a = calculate_acceleration(r_current, M_sun)

    # 5b. Update velocity (v_current) based on acceleration (a) and time step (DT)
    # This is the 'Euler' part of the update for velocity
    v_current = v_current + a * DT

    # 5c. Update position (r_current) using the *newly calculated* velocity (v_current)
    # This crucial step makes it the 'Cromer' or Semi-Implicit Euler method,
    # which ensures better energy conservation than standard Euler.
    r_current = r_current + v_current * DT

    # 5d. Store the new position vector for plotting
    positions[i] = r_current

print("Simulation complete.")

# --- 6. VISUALIZATION ---

# Scale the positions back to Astronomical Units (AU) for readable plotting
x_au = positions[:, 0] / R_AU
y_au = positions[:, 1] / R_AU

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

# Formatting and Display
plt.title(f'Orbital Simulation using Euler-Cromer Integration ({NUM_STEPS} steps)')
plt.xlabel('X Position (AU)')
plt.ylabel('Y Position (AU)')
# Crucially, ensure the aspect ratio is equal so the orbit isn't distorted
plt.gca().set_aspect('equal', adjustable='box')
plt.grid(True, linestyle=':', alpha=0.7)
plt.legend()
plt.show()

Code Breakdown: How It Works

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

Vector Math: We treat position and velocity as vectors (arrays with an x and y component). This allows us to calculate the direction of the gravitational pull easily. np.linalg.norm is the key function here—it calculates the distance rr between the planet and the Sun.
The Acceleration Function: This function implements Newton's law. It takes the position, calculates the distance, finds the unit vector (direction), and scales it by the magnitude of gravity. The negative sign is crucial—it ensures the acceleration points towards the Sun (attraction).
The Main Loop: This is where the magic happens.
- We calculate acceleration based on where the planet is.
- We update velocity.
- The Euler-Cromer Trick: We immediately use that new velocity to update the position. This small change prevents the simulation from "drifting" over time.
Visualization: We convert the raw meters back to Astronomical Units (AU) for a human-readable graph. The command 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.

Conclusion: The Digital Cosmos

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.

Let's Discuss

In the code, we treat the Sun as stationary (infinitely massive). How would you modify the Euler-Cromer loop to account for the Sun's motion if it were not stationary (the true Two-Body Problem)?
If we changed the time step 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

推荐订阅源

DEV Community