# Astrophysics & AI with Python: Simulating Planetary Orbits with Kepler's Laws > Source: > Published: 2026-06-16 20:00:00+00:00 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. ``` python #!/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() ``` 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](http://tiny.cc/PythonAstrophysics). Check all the other 50 Programming & AI ebooks with python, typescript, swift, c#: [here](http://tiny.cc/ProgrammingBooks)