Member-only story
Python in Astrophysics: Creating an N-Body Simulation
How do you know what happens when two galaxies collide? There are a few options here. The first is to just look into the sky and find different galaxies that are at some point of the collision process. Another option would be to just wait and watch (but you will be waiting for a very long time). Finally, you could model a collision.
Modeling the interaction between two galaxies is no small feat. Each galaxy can contain more than a billion stars and all of these stars have gravitational interactions with other stars. But still, the idea isn’t too bad. I just need to use the gravitational forces between all the stars (in both galaxies) to find out how the stars move.
This is what we call the n-body problem. It’s not just used for galaxies — but you can use it for any system containing lots of elements that all interact with each other (including particles in a gas).
I’m going to build this in python (and so can you) but let me make an important point. This is an educational exercise. This is not a professional n-body simulator. I’m going to do this from the most basic principles and try to use very little programming tricks.
The Three Body Problem
Before getting too deep — let’s just think about the case of three masses interacting (it’s still an n-body problem but n = 3). Also, just to make things simpler, I am going to use simple units such that the gravitational constant is G = 1. The mass of a star will be M = 1 and the orbital distance will be R = 1. This will let me use normal time steps on the order of fractions of a second. If you want to use different units, that’s up to you.
Here are three stars:
There is a gravitational interaction between pairs of stars. That means that we will need to calculate three forces: the force between 1 and 2, 2 and 3, and 3 and 1. Each star will then have a net force equal to the vector sum of two forces (from the two other stars). If we calculate the relative position between stars (say 1 and 2) as r12, then the force on star 2 would be:
