Gravitational Simulation | Python | Astronomy | Physics | Computational Modeling

Introduction

In this project, we will develop a gravitational simulation using Python, showcasing the interaction between Earth and a hypothetical "Super Jupiter," a planet with 500 times the mass of Jupiter. This simulation aims to demonstrate the impact of such a massive body on Earth's motion, considering the immense gravitational forces at play. This project suits students and hobbyists passionate about physics, astronomy, and computational simulations. To achieve this, we will employ Python libraries like NumPy for numerical calculations and Matplotlib for visualizing the dynamic movements of the planets.

👀 Preview

🎯 Tasks

In this project, you will learn:

How to understand and apply Newton's Law of Universal Gravitation to model the interaction between celestial bodies.
How to use Python programming to create a computational model of a gravitational system.
How to employ the NumPy library for efficient numerical calculations in Python.
How to simulate the orbital mechanics of Earth in the presence of a "Super Jupiter" with 500 times the mass of Jupiter.
How to analyze and interpret the results of the simulation to understand the impact of massive celestial bodies on orbital dynamics.
How to implement Matplotlib to create visual representations of the simulation, showcasing the orbital paths and relative positions of the planets.
How to explore the concepts of force, mass, and acceleration in a cosmic context.
How to fine-tune simulation parameters like mass, distance, and time steps for different scenarios.
How to develop skills in debugging and optimizing Python code for scientific computations.

🏆 Achievements

After completing this project, you will be able to:

Apply fundamental principles of physics, specifically Newton's Law of Universal Gravitation, in a practical, computational context.
Create and run a physics-based simulation using Python.
Demonstrate proficiency in using NumPy for handling large-scale numerical computations efficiently.
Visualize complex data and simulations using Matplotlib, enhancing the interpretability of scientific results.
Understand the dynamics of planetary motion and the effects of gravitational forces from massive bodies.
Analyze and interpret the results of the simulation to draw meaningful conclusions about celestial mechanics.
Adjust and experiment with simulation parameters, leading to a deeper understanding of orbital mechanics.
Showcase enhanced problem-solving and debugging skills in a programming environment, particularly in the context of scientific computing.
Demonstrate a foundational knowledge of how gravitational forces shape the motion of celestial bodies, paving the way for further exploration in astrophysics and computational modeling.

Setting Up the Environment

First, we need to set up our project files and import necessary libraries.

Create a Python file named simulation.py:

touch ~/project/simulation.py

Then import the following libraries:

## Importing required libraries
from typing import Tuple
import numpy as np
import pylab as py
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.lines import Line2D
from tqdm import trange

In this initial step of the gravitational simulation project, several important Python libraries are imported to set up the environment for our simulation:

typing: Specifically, Tuple from the typing module is imported. This is used for type hinting in Python, which helps with code readability and debugging. It allows specifying what types of values can be grouped in a tuple, ensuring consistency and clarity in the function definitions.
numpy: The numpy library, imported as np, is a fundamental package for scientific computing in Python. It provides support for large, multi-dimensional arrays and matrices, along with a vast collection of high-level mathematical functions to operate on these arrays. It's crucial for handling numerical calculations in the simulation.
pylab: Imported as py, pylab is a module in Matplotlib that gets installed alongside Matplotlib. It provides a MATLAB-like interface, particularly useful for interactive calculations and plotting. However, its usage is generally discouraged in favor of explicitly importing matplotlib.pyplot and numpy.
matplotlib.pyplot: This module, imported as plt, is used for creating static, animated, and interactive visualizations in Python. It's a core aspect of this project for plotting the trajectories of celestial bodies.
matplotlib.animation: This submodule of Matplotlib provides the functionality for creating animations. It's essential for visualizing the dynamic movement of planets in the gravitational simulation.
matplotlib.lines.Line2D: This class from the Matplotlib library is used for creating line objects that can be drawn on plots. It's useful for more detailed or customized plotting in the simulation.
tqdm.trange: trange from the tqdm module is a variant of Python's range function with a built-in progress meter. It's useful for displaying the progress of loops in the simulation, especially when dealing with extensive calculations.

By importing these libraries and modules, we lay the groundwork for the complex numerical computations and visualizations required in the gravitational simulation project. This setup is critical for the subsequent steps where the actual physics and animation logic will be implemented.

Defining Constants

Define the constants that will be used in our simulation. This includes gravitational constant, astronomical units, and normalized masses of the Earth, Jupiter, and the Sun.

## Constants
G = 6.673e-11                 ## Gravitational Constant
AU = 1.496e11                 ## Astronomical Unit in km
YEAR = 365*24*60*60.0         ## Seconds in one year
MM = 6e24                     ## Normalizing mass
ME = 6e24/MM                  ## Normalized mass of Earth
MS = 2e30/MM                  ## Normalized mass of Sun
MJ = 500*1.9e27/MM            ## Normalized mass of Jupiter
GG = (MM*G*YEAR**2)/(AU**3)   ## Gravitational constant for simulation

In this step, we define several essential constants for the gravitational simulation. These constants are fundamental to accurately modeling the celestial mechanics and gravitational interactions in the simulation:

Gravitational Constant (G): Set to 6.673e-11 (in SI units), this constant is crucial in Newton's law of universal gravitation. It represents the strength of gravity and is used to calculate the gravitational force between two masses.
Astronomical Unit (AU): Defined as 1.496e11 kilometers, it represents the average distance from the Earth to the Sun. This unit is used to express distances within the solar system in a more manageable way.
Year (YEAR): Calculated as 365*24*60*60.0, it represents the number of seconds in a year. This is used to convert time-related calculations into years, a more intuitive time unit for astronomical simulations.
Normalizing Mass (MM): Set to 6e24, this value is used as a reference mass to normalize other masses in the simulation. It’s roughly equivalent to Earth's mass.
Normalized Mass of Earth (ME): Calculated as 6e24/MM, it represents Earth's mass in normalized units. It's essentially set to 1, as the normalization mass (MM) is based on Earth's mass.
Normalized Mass of Sun (MS): Set to 2e30/MM, it represents the Sun's mass in normalized units. The Sun's mass is a critical factor in simulating the gravitational pull on the planets.
Normalized Mass of Jupiter (MJ): Uniquely, in this simulation, Jupiter's mass is amplified to 500*1.9e27/MM to represent a "Super Jupiter" scenario. This significantly larger mass will demonstrate the impact of a much more massive planet on Earth's orbit.
Gravitational Constant for Simulation (GG): This is the gravitational constant adjusted for the simulation's scale and units, calculated as (MM*G*YEAR**2)/(AU**3). It adapts the real-world gravitational constant to the normalized units used in the simulation.

By defining these constants, we are equipped with the necessary parameters to accurately simulate the gravitational forces and movements of celestial bodies like Earth and this hypothetical "Super Jupiter." This step is vital for ensuring that the simulation's calculations reflect realistic astronomical dynamics. At the same time, you can modify these constants as much as you want to observe different results at the end.

Creating the Gravitational Force Function

Now, we will define a function to calculate the gravitational force between two bodies. This function will be crucial in determining the motion of the planets.

## Function to calculate gravitational force
def gravitational_force(m1: float, m2: float, r: np.ndarray) -> np.ndarray:
    """
    Calculate gravitational force between two bodies.
    """
    F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2
    theta = np.arctan2(np.abs(r[1]), np.abs(r[0]) + 1e-20)
    F = F_mag * np.array([np.cos(theta), np.sin(theta)])
    F *= -np.sign(r)
    return F

In this step, we define a function gravitational_force to calculate the gravitational force between two celestial bodies. This function is a critical component of the simulation, as it applies Newton's law of universal gravitation to determine the force exerted between any two masses in space. Here's a breakdown of the function:

Function Definition: The function gravitational_force takes three arguments:
- m1 and m2: The masses of the two bodies (float values).
- r: A NumPy array representing the displacement vector between the two bodies.
Calculating Force Magnitude (F_mag):
- The function first calculates the magnitude of the gravitational force using the formula F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2.
- GG is the gravitational constant for the simulation.
- np.linalg.norm(r) computes the Euclidean distance between the two bodies. A small number (1e-20) is added to prevent division by zero.
Determining Direction (theta and F):
- The angle theta is calculated using np.arctan2 to find the angle of the force vector in polar coordinates. Again, a small number (1e-20) is added to the denominator to avoid division by zero when the displacement in the x-direction (r[0]) is zero.
- The force vector F is then calculated using F_mag and the angle theta, with components along the x and y axes ([np.cos(theta), np.sin(theta)]).
Adjusting Force Direction:
- The force vector F is multiplied by -np.sign(r) to ensure that the force is always attractive (pointing towards the other body) as per gravitational laws.
Returning the Force Vector:
- Finally, the function returns the gravitational force vector F.

By implementing this function, we learn how to apply physical laws in a computational environment. The function is essential for calculating the gravitational interactions between Earth and the "Super Jupiter" in the simulation, which will influence their respective trajectories.

Implementing the RK4 Solver

Implement the Runge-Kutta 4th order solver to solve the differential equations of motion. This method is widely used for solving ordinary differential equations numerically.

## RK4 Solver
def RK4Solver(t: float, r: np.ndarray, v: np.ndarray, h: float, planet: str, r_other: np.ndarray, v_other: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    Fourth order Runge-Kutta solver for planetary motion.
    """
    def dr_dt(v: np.ndarray) -> np.ndarray:
        return v

    def dv_dt(r: np.ndarray, planet: str) -> np.ndarray:
        if planet == 'earth':
            return (gravitational_force(ME, MS, r) + gravitational_force(ME, MJ, r - r_other)) / ME
        elif planet == 'jupiter':
            return (gravitational_force(MJ, MS, r) - gravitational_force(MJ, ME, r - r_other)) / MJ

    k11 = dr_dt(v)
    k21 = dv_dt(r, planet)

    k12 = dr_dt(v + 0.5 * h * k21)
    k22 = dv_dt(r + 0.5 * h * k11, planet)

    k13 = dr_dt(v + 0.5 * h * k22)
    k23 = dv_dt(r + 0.5 * h * k12, planet)

    k14 = dr_dt(v + h * k23)
    k24 = dv_dt(r + h * k13, planet)

    y0 = r + h * (k11 + 2 * k12 + 2 * k13 + k14) / 6
    y1 = v + h * (k21 + 2 * k22 + 2 * k23 + k24) / 6

    return y0, y1

In this step, we implement a function named RK4Solver, which uses the Fourth Order Runge-Kutta (RK4) method for solving differential equations. This method is crucial for accurately simulating the motion of planets in our gravitational model. Here's an overview of what the RK4Solver function does:

Function Definition:
- RK4Solver takes several parameters:
  - t: The current time.
  - r: The current position vector of the planet.
  - v: The current velocity vector of the planet.
  - h: The time step for the simulation.
  - planet: A string indicating which planet (Earth or Jupiter) is being simulated.
  - r_other and v_other: The position and velocity vectors of the other planet.
Helper Functions within RK4Solver:
- dr_dt(v: np.ndarray) -> np.ndarray: This function returns the rate of change of position, which is simply the velocity v.
- dv_dt(r: np.ndarray, planet: str) -> np.ndarray: This function calculates the acceleration of the planet due to gravitational forces. It uses the gravitational_force function defined earlier to compute the forces exerted by the Sun and the other planet (Earth or Jupiter) and returns the acceleration.
Runge-Kutta Calculations:
- The RK4 method involves calculating four "slopes" (k11, k21, k12, k22, k13, k23, k14, k24) at different points within the time step and then combining them to get an accurate estimate of the position and velocity at the next time step.
- These slopes are calculated using the current position and velocity, as well as the derivatives of position (dr_dt) and velocity (dv_dt).
Updating Position and Velocity:
- The function calculates the next position (y0) and velocity (y1) of the planet using the weighted average of these slopes, according to the RK4 formula.
Returning the Results:
- Finally, the function returns the updated position and velocity vectors (y0, y1) for the planet.

By implementing the RK4Solver function, we learn an efficient and accurate method for solving ordinary differential equations numerically. This is particularly important in simulations involving complex systems like planetary motion, where precision is key to obtaining realistic results. The RK4 method strikes a good balance between computational efficiency and accuracy, making it a popular choice in many scientific and engineering applications.

Setting Up the Animation

Before running the simulation, we need to set up the animation. This step involves creating a plot and initializing the lines and markers representing the planets.

## Setup animation
def setup_animation() -> Tuple[py.Figure, py.Axes, Line2D, Line2D, py.Text]:
    """
    Set up the animation plot.
    """
    ## Creating a Plot Figure and Axes
    fig, ax = py.subplots()

    ## Setting Axes Limits and Ticks
    ax.axis('square')
    ax.set_xlim((-7.2, 7.2))
    ax.set_ylim((-7.2, 7.2))
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])

    ## Plotting the Sun
    ax.plot(0, 0, 'o', markersize=9, markerfacecolor="#FDB813",
            markeredgecolor="#FD7813")

    ## Initializing Lines for Earth and Jupiter
    line_earth, = ax.plot([], [], 'o-', color='#d2eeff',
                          markevery=10000, markerfacecolor='#0077BE', lw=2)
    line_jupiter, = ax.plot([], [], 'o-', color='#e3dccb', markersize=8,
                            markerfacecolor='#f66338', lw=2, markevery=10000)
    ## Adding a Text Object
    ttl = ax.text(0.24, 1.05, '', transform=ax.transAxes, va='center')

    ## Returning Components
    return fig, ax, line_earth, line_jupiter, ttl

In this step, we define a function called setup_animation that prepares the plot for animating the gravitational simulation. This function sets up the visual components of the simulation, such as the plot area, planetary markers, and initial positioning.

This function is essential for visualizing the simulated orbits of Earth and Jupiter. It sets the stage for the dynamic part of the simulation where the positions of these planets will be updated frame by frame to create an animation.

Creating the Animation Function

Create a function that will update the positions of the planets for each frame of the animation.

## Animation function
def animate(i: int) -> Tuple[Line2D, Line2D, py.Text]:
    """
    Animation function for the planetary motion.
    """
    earth_trail, jupiter_trail = 40, 200
    tm_yr = 'Elapsed time = ' + str(round(t[i], 1)) + ' years'
    ttl.set_text(tm_yr)
    line_earth.set_data(r[i:max(1, i - earth_trail):-1, 0],
                        r[i:max(1, i - earth_trail):-1, 1])
    line_jupiter.set_data(r_jupiter[i:max(
        1, i - jupiter_trail):-1, 0], r_jupiter[i:max(1, i - jupiter_trail):-1, 1])
    return line_earth, line_jupiter, ttl

In this step, we define the animate function, which is the core of creating the animation for the planetary motion simulation. This function is called repeatedly to update the positions of Earth and Jupiter on the plot, creating the effect of motion over time. Here's what happens in the animate function:

Trail Lengths for Earth and Jupiter:
- The variables earth_trail and jupiter_trail are set, determining how many previous positions (trail) of each planet are shown. This creates a visual trail effect, showing the path each planet has traveled.
Updating Time Display:
- The elapsed time (tm_yr) is calculated and set as the text of the ttl text object. This displays the simulation time in years on the animation.
Updating Earth's Position:
- line_earth.set_data updates the position of Earth in the animation. It uses a slice of the position array r to create a trail effect, showing where Earth has been recently in its orbit.
Updating Jupiter's Position:
- Similarly, line_jupiter.set_data updates Jupiter's position using a slice of the r_jupiter array. The longer trail length for Jupiter reflects its larger orbit.
Returning Updated Objects:
- The function returns the updated Line2D objects (line_earth, line_jupiter) and the text object (ttl). These are the elements that change in each frame of the animation.

This animate function is crucial for visualizing the simulation, showing how Earth and Jupiter move over time under their mutual gravitational influences. It brings the simulation to life by dynamically updating the plot with each frame.

Initializing Simulation Parameters

Before starting the simulation, initialize the time parameters, position, and velocity of Earth and Jupiter.

## Initialization
ti, tf = 0, 120  ## Initial and final time in years
N = 100 * tf     ## 100 points per year
t = np.linspace(ti, tf, N)  ## Time array
h = t[1] - t[0]  ## Time step

## Position and Velocity Initialization
r = np.zeros([N, 2])         ## Position of Earth
v = np.zeros([N, 2])         ## Velocity of Earth
r_jupiter = np.zeros([N, 2])  ## Position of Jupiter
v_jupiter = np.zeros([N, 2])  ## Velocity of Jupiter

## Initial Conditions
r[0] = [1496e8 / AU, 0]
r_jupiter[0] = [5.2, 0]
v[0] = [0, np.sqrt(MS * GG / r[0, 0])]
v_jupiter[0] = [0, 13.06e3 * YEAR / AU]

In this step, we set up the initial conditions and parameters for the simulation, which is crucial for starting the animation of Earth and Jupiter's motion. Here's a summary of what is being done:

Time Initialization:
- ti and tf are set as the initial and final times of the simulation, measured in years (0 to 120 years).
- N is defined as the total number of points in the simulation, calculated as 100 points per year.
- t is an array created using np.linspace to represent the time steps from the initial to the final time.
- h is the time step, calculated as the difference between the first two elements of the t array.
Position and Velocity Initialization:
- Arrays r and v are initialized to store the position and velocity of Earth at each time step. They are zero-initialized and have a shape to accommodate two coordinates (x and y) for each time step.
- Similarly, r_jupiter and v_jupiter are initialized to store Jupiter's position and velocity.
Setting Initial Conditions:
- The initial position of Earth (r[0]) is set to its distance from the Sun, normalized by the Astronomical Unit (AU).
- The initial position of Jupiter (r_jupiter[0]) is set to 5.2 AU from the Sun, reflecting its actual position in the solar system.
- The initial velocity of Earth (v[0]) is calculated based on the gravitational force exerted by the Sun.
- The initial velocity of Jupiter (v_jupiter[0]) is set to a specific value that reflects its actual orbital velocity.

This step is fundamental as it establishes the starting point for the simulation. The initial positions and velocities are crucial for calculating the subsequent motion of both planets under gravitational forces.

Running the Simulation

Execute the simulation over the defined time frame using the RK4 solver and update the positions of the planets.

## Running the simulation
for i in trange(N - 1, desc="Generating Animation"):
    r[i + 1], v[i + 1] = RK4Solver(t[i], r[i],
                                   v[i], h, 'earth', r_jupiter[i], v_jupiter[i])
    r_jupiter[i + 1], v_jupiter[i +
                                1] = RK4Solver(t[i], r_jupiter[i], v_jupiter[i], h, 'jupiter', r[i], v[i])

In this step, we execute the main loop to run the gravitational simulation. This is where the motion of Earth and Jupiter is computed over the specified time period using the RK4Solver function. Here's a simple explanation of this process:

Running the Simulation Loop:
- A loop is set to run for N - 1 iterations, where N is the total number of time steps in the simulation. This loop is crucial for advancing the simulation through each time step.
- The trange function from the tqdm module is used to iterate over the time steps. This provides a progress bar (described as "Generating Animation") to show the simulation's progress.
Updating Earth and Jupiter's Positions and Velocities:
- Within each iteration of the loop, the RK4Solver function is called twice:
  - Once to update Earth's position (r[i + 1]) and velocity (v[i + 1]) at the next time step.
  - Once to update Jupiter's position (r_jupiter[i + 1]) and velocity (v_jupiter[i + 1]) at the next time step.
- The RK4Solver function takes the current position and velocity of the planet, along with the position and velocity of the other planet, to calculate the new state.
Simulating Interactions:
- This loop essentially simulates the gravitational interaction between Earth and Jupiter (considering their influence on each other and the influence of the Sun) at each small time step, leading to a realistic representation of their orbits over time.

By the end of this loop, we have arrays r, v, r_jupiter, and v_jupiter filled with the positions and velocities of Earth and Jupiter at each time step, ready to be animated in the next phase of the project. This step is where the core physics and numerical methods are applied to simulate celestial dynamics.

Displaying the Animation

Finally, run the animation and display the results.

## Setting Up the Animation
fig, ax, line_earth, line_jupiter, ttl = setup_animation()
## Adding Scale and Labels
ax.plot([-6,-5],[6.5,6.5],'r-')
ax.text(-4.5,6.3,r'1 AU = $1.496 \times 10^8$ km')

ax.plot(-6,-6.2,'o', color = '#d2eeff', markerfacecolor = '#0077BE')
ax.text(-5.5,-6.4,'Earth')

ax.plot(-3.3,-6.2,'o', color = '#e3dccb',markersize = 8, markerfacecolor = '#f66338')
ax.text(-2.9,-6.4,'Super Jupiter (500x mass)')

ax.plot(5,-6.2,'o', markersize = 9, markerfacecolor = "#FDB813",markeredgecolor ="#FD7813")
ax.text(5.5,-6.4,'Sun')

## Creating the Animation
anim = animation.FuncAnimation(
    fig, animate, frames=4000, interval=1, blit=False)

## Displaying the Animation
plt.show()

In the final step of the project, we display the animation of the gravitational simulation that shows the motion of Earth and the hypothetical "Super Jupiter". This step involves visualizing the calculated positions of the planets over the simulation period.

Now we've completed all the steps, we can run the code in the desktop environment by using the following command:

cd ~/project
python simulation.py

This step is where the we can visually appreciate the results of their simulation, observing how the massive "Super Jupiter" influences Earth's orbit.

Summary

In this project, we have effectively developed a gravitational simulation featuring Earth and Jupiter using Python. We employed the Runge-Kutta 4th order (RK4) method for precise numerical integration, and Matplotlib for clear and engaging visualization. This project not only offers a hands-on application of physics and coding principles but also encourages students to experiment by altering the code to simulate celestial bodies of various masses, thereby exploring a range of gravitational interactions and dynamics.

Gravitational Simulation of Earth and Super Jupiter