Gravitron: A Rudimentary N-body Physics Simulator

While bored during a zoom meeting one day, I began to throw together some code in Python to simulate Newtonian gravitational interactions using the basic kinematic equations and the equation of universal gravitation. This may sound complicated, but it really only involves three or four different formulae. The first, and most obvious, is Newton’s law of universal gravitation:

\(\Large F_{1 on 2} = F_{2 on 1} = G \frac{m_1 m_2}{r^2}\)

A staple of high school physics, all this formula states is that the force of gravitational attraction between two given bodies is a function of their mass and distance scaled to the gravitational constant \(G\). It doesn’t account for things like relativity, but it is sufficient in most cases. So, in order to determine the sum of all forces on a body in a purely gravitational universe, we simply need to use this equation to calculate the force vector on every body. Since we know that \(F = ma\), and we know both the mass of every body and the sum of all forces exerted on it, we may find its position and velocity using the fundamental kinematic equations:

\(\Large v = v_0 + at\)

\(\Large \Delta x = v_0 t + \frac{1}{2}at^2\)

The code for this is implemented in only a few simple functions (helper functions for vector manipulation are omitted):

def fg(a, b): #grav formula (returns force vector)
    g = 6.674e-11 #grav constant
    return g*(a.getMass() * b.getMass())/(distance(a,b) ** 2)

def appendForce(self, obj):
    self.forces.append((fg(self, obj), angle(self, obj)))
    
def getNetForce(self):
    forceSum = (0,0)
    
    for force in self.forces:
        forceSum = addPolarVectors(forceSum, force)
        
    return forceSum

def calculateAccel(self):
    f = self.getNetForce()
    self.accel = (f[0]/self.mass, f[1])
    return self.accel
    
def calculateDisp(self, tick): #kinematic equation #3 
    v = (self.velocity[0] * tick, self.velocity[1])
   
    return addPolarVectors(v, ((0.5 * self.accel[0] * (tick ** 2)), self.accel[1]))
    
def calculateVelocity(self, tick): #calculate body velocity based on tick interval (kinematic equation #1)
    self.velocity = addPolarVectors(self.velocity, (self.accel[0] * tick, self.accel[1]))
    return self.velocity

After this is implemented, we can finally simulate bodies. As a test, I simply printed the positions and velocities of each body in the terminal, but I wanted a visual reference for the simulation. I first threw together a simple console display that put the bodies in a 2D array and printed it to the terminal, but I wanted something a little more polished. For this, I turned to Pygame. Pygame is a fairly popular Python library that provides SDL bindings, allowing you to draw 2D objects in a window. In order to visualize the locations of bodies, we can simply draw a circle at a place on the screen corresponding to the body’s location in the simulation. After doing this, we are finally able to see our objects moving. Shown below is a basic simulation of a binary star system and one orbiting planet:

As we can see, just these three or four simple equations and a day of programming is enough to get a realistic looking simulation of a star system. Now that I was able to simulate bodies, the next step was to add some more customization to the program. Originally I imported simulation data using JSON files that I manually coded. This worked well enough to test the program, but I wanted a way to add objects and save simulations completely in-app. To do this, I threw together a simple GUI that would pop up whenever the user hit the escape key. Since Pygame has no built-in support for GUIs, I had to code each element from scratch. The first two elements that I added were a textbox that a user could enter text and numbers into, and a button that could be used to trigger events. Adding these to the pause menu, I now had a way to create new objects in the simulation:

Here you can see all the parameters of a body. Mass, X/Y position, and velocity vector parameters are all part of the physics backend, while the object radius and RGB values affect how it is displayed in the window. The simulation can be zoomed in and out between different pixel/meter scale levels using the +/- keys on a keyboard, and the camera can be moved around using the arrow keys. The amount of space any object occupies on the screen is therefore scaled to this value. However, if the view becomes so far zoomed out that a given object would otherwise become invisible, its radius is defaulted to a minimum of 2 pixels to enable it to still be visible. At the solar system level, every object’s apparent size is below the minimum simply due to the vast amount of space between objects, however it is possible to see at the planetary level, like with this simulation of the fictional gas giant Jool from the video game Kerbal Space Program:

Another problem with simulating gravitational interactions at different scales is the difference in simulation time. For example, it takes the Moon only around 30 days to orbit the Earth, covering a distance of about 1.5 million miles over the course of its orbit. However, the Earth takes roughly 365.25 days to orbit the Sun, covering a distance of almost 600 million miles. Therefore, it is easy to see that simulating one Earth orbit is much more computationally taxing if we use the same time interval between each calculation. This is where having variable tick rates comes in handy. In this simulation, the position and velocity of each body is re-calculated after a “tick”, which just a constant arbitrary distance in time. In between ticks, objects in this universe are assumed to travel along their velocity vectors at a constant velocity (a straight line). Therefore, it is easy to see that a smaller tick interval will give us greater accuracy in our simulation since objects will re-calculate their gravitational interactions more often. Conversely, a larger tick interval will result in orbits that eventually become farther and farther apart in the simulation as the objects will move along their tangents longer before re-calculating their interactions, therefore resulting in a slightly weaker interaction than there should be. To overcome this, I coded in a slider class to the GUI so I could add a tick rate and frame rate slider to the simulation. The frame rate slider has no effect on the underlying simulation but it does determine how much simulation time passes before it delivers a new frame to the window. This can be useful for speeding up simulations without impacting their precision.

Finally, I added a means of loading and saving simulations to JSON files, as well as a little display in the pause menu where users can see all the bodies in the simulation and remove them individually if they wish. Here is what the final menu screen looks like:

Using this, I was able to create some very nice simulations, however I still had one problem. The camera’s frame of reference was fixed, meaning that it couldn’t track objects. This was ok if I wanted to simulate just a solar system or a moon orbiting a single planet, although it made things difficult if I wanted to look at moons orbiting planets that were in turn orbiting stars. My solution to this was to allow the camera to track objects simply by clicking on them. This works really well, as can be seen below in this simulation of a solar system: