My freshman year of college, I asked a physics professor if I could do some research with him, and he had me spend my free time programming a simulation of two bodies in motion using Newtonian mechanics… in Fortran. Yes, Fortran. So a few days ago I said to myself, hey, I can do this better in Python. So I did.
You can find my code on github. I use Python 2.7.6 and pylab (which is a part of matplotlib). The file labeled “nbodyproblem” has the classes and functions you need to make a simulation (I allowed arbitrary number of planets, even though I know that the nature of a chaotic system is that the model will probably quickly stop predicting accurately). The file “nbodysimulations” has some of the simulations that I had fun with.
For example, if you GraphSunEarth(), it runs a simulation of the sun and earth for 100 days (recalculating every minute by default), and you see a picture:
The yellow dot being the (very uninteresting) path of the sun, and the blue line being the path of the earth.
Or maybe this wasn’t you wanted, you want to see the path of the earth and sun for a longer period of time. So you run resetSunEarth() (or don’t — this resets them to the same starting positions, which you don’t need to do for any particular reason). If you run GraphSunEarth(24*300, 60), which will run the simulation for 300 days, and recalculate every hour rather than every minute. This gives a picture:
Okay, maybe we want more planets! We could run AMarsSimulation(), which runs the sun, earth, and mars for (by default) a year with 1 hour increments. This gives the picture:
Well that’s a bit boring. What about the sun, earth, and moon? This is a little problematic when graphing, because any picture that includes the sun makes the earth-moon system look like one point when the graph is large enough to include the sun. So I also graphed the close-up of just the earth and moon.
(Note that the places where the paths cross are not collisions, since the bodies are at that point at different points in time. As currently written, the code would give an error if there is a collision.)
The sun and a comet (i.e. a body with the wrong velocity) can be simulated with Comet():
Note the colors here were randomly generated.
Two suns (i.e. two bodies with large masses) can be simulated with TwoSuns():
But if the two suns don’t have nice starting velocities, weird things can happen. CrazySuns() simulations this:
The functions in nbodyproblem also make it easy to create new simulations. Just:
- Use the “planet” class to make your new bodies
- the “simulation” function inputs a list of planets and number of timesteps (and an optional input which allows you to recalculate less often), and outputs
- the “graphData” function inputs the output from simulation (and optionally, a list of colors to make the planets) and makes a graph.