Building an Interactive Celestial Body Simulation (Part 1)

I recently developed a program that simulates the attractive force of two objects in space time using Newtonian Mechanics. In this post, I will walk you through my process of building this program.

0. Preface

My original (and current) idea for this project was/is to create a video game-like simulation software for users to see how any number of bodies in space would interact with each other. The current iteration is not quite there yet, I still have to refine some bugs, and develop an interactive platform that makes it easy for even a young child to play with.

1. Defining the Problem

The goal of this project is to create a simulation that models the gravitational forces between two celestial bodies and tracks their motion over time. I knew it would be more of a challenge for me to code the project than it would be to define the mathematics, so I brushed up on some Python libraries that I planned to use, like matplotlib.

First, let's start with what I know best, the mechanics: I am using Newtonian Mechanics for this problem even though they are known to not work for very large masses, like black holes. But since I don't know General Relativity, this is the best it's gonna get (for now).

The majority of the heavy lifting is done by the Inverse Square Law:

\vec{\mathbf{F}} = G \frac{m_1 m_2}{\vec{\mathbf{R}}^2}

where

\vec{\mathbf{R}} = x\vec{\mathbf{i}} + y\vec{\mathbf{j}}

Next we define some basic kinematic equations, which are updated with small time steps given by the Inverse Square Law defined above. I used:

\vec{\mathbf{v}}(t) = \vec{\mathbf{a}}t

and

\vec{\mathbf{r}}(t) = \vec{\mathbf{v}}t + \frac{1}{2} \vec{\mathbf{a}}t^2

For sake of completeness, I will also define how I computed the x and y components of the force vector:

F_x = \vec{\mathbf{F}} \frac{x_2 - x_1}{|\vec{\mathbf{R}}|}

and

F_y = \vec{\mathbf{F}} \frac{y_2 - y_1}{|\vec{\mathbf{R}}|}

Where \( x_2 \) and \( y_2 \) correspond to the position of the second body, and \( x_1 \) and \( y_1 \) to the first body.

Now let's move onto the code...

2. Designing the Classes and Methods

The key part of the design is the CelestialBody class. Each celestial body has a mass, position, and velocity. These are initialized in the constructor, and the class includes methods for updating the body's position and velocity:


class CelestialBody:
    def __init__(self, mass, position, velocity):
        self.mass = mass
        self.position = np.array(position, dtype=float)
        self.velocity = np.array(velocity, dtype=float)
    
    def update_position(self, dt, acceleration):
        self.position += self.velocity * dt + 0.5 * acceleration * dt ** 2

    def update_velocity(self, dt, acceleration_new, acceleration_old):
        self.velocity += 0.5 * (acceleration_old + acceleration_new) * dt

The update_position method, as its name suggests, updates the position based on velocity and acceleration, while the update_velocity method updates the velocity using the average of old and new accelerations.

3. Gravitational Force Calculation

The next crucial part of the simulation was computing the gravitational force between two bodies.


def compute_gravitational_force(body1, body2, G):
    r_vec = body2.position - body1.position
    r_mag = np.linalg.norm(r_vec) + 1e-10  # Prevent division by zero
    F_mag = G * body1.mass * body2.mass / r_mag ** 2
    F_vec = F_mag * (r_vec / r_mag)
    return F_vec

This uses the Inverse Square Law which we defined earlier.

4. Integrating the Motion

Once the forces are calculated, I used the Verlet integration method to update the positions and velocities of the celestial bodies step by step. The step function performs these calculations and updates for each time step:


def step(dt, body1, body2, G):
    F12 = compute_gravitational_force(body1, body2, G)
    a1_old = F12 / body1.mass
    a2_old = -F12 / body2.mass

    body1.update_position(dt, a1_old)
    body2.update_position(dt, a2_old)

    F12_new = compute_gravitational_force(body1, body2, G)
    a1_new = F12_new / body1.mass
    a2_new = -F12_new / body2.mass

    body1.update_velocity(dt, a1_new, a1_old)
    body2.update_velocity(dt, a2_new, a2_old)

This approach models the interactions step-by-step in a loop, allowing us to animate the motion of the bodies.

5. Animating the Simulation

To visualize the motion, I used matplotlib.animation to create an animation. The positions of the celestial bodies are updated over time, and trails are drawn to show their paths:


anim = FuncAnimation(fig, animate, init_func=init, interval=20, blit=True)
plt.show()

The result is a dynamic animation showing the gravitational dance between the two bodies. The animation also includes trails to track the movement of the Moon relative to Earth.

In The End...

This project was a great exercise in understanding gravitational dynamics and implementing numerical integration techniques, but it is far from done. I still need to implement support for multiple bodies as well as an interactive front end for ease of use. Plan to see more blog posts about this project in the future. To see the full project, you can visit the repository on my GitHub.