Simulating a real solar system with 70 lines of Python code
Animated Science
Computational astrophysics is a really fun subject where we use computers to simulate astronomical objects and phenomena. To demonstrate the power and beauty of computer simulations in the study of astronomy, I show here a program I wrote with 70 lines of Python code that simulates a (somewhat) realistic solar system.
We know that planets follow orbits determined mainly by the gravity from the Sun. To calculate the orbits of a planet, we need to calculate the gravitational force from the Sun and integrate it to get the velocity and position. We learned in general physics that the force following Newton’s law of gravity combined with Newton’s second law of motion gives the acceleration:
We use this formula to evolve the position and velocity of a planet:
This is a modified version of the Euler’s method of integration, where the first equation is forward and the second equation is backward. Unlike the normal Euler’s method, this modified version is stable.
The benefit of programming in Python is that Python handles vectors neatly. The three vector equations, each with three components, can be easily converted into the following code:
These three lines of code determine how most of the objects in the solar system move.
Now we have the integrator, we need the initial conditions, i.e. the position and velocity of a planet at some date. The precise locations and velocities of a planet are obtained from the NASA JPL Horizons online solar system database (https://ssd.jpl.nasa.gov/?horizons). This is done by calling a Python module that sends queries to the database and obtains the data. This part of the code looks like this:
where id=1 corresponds to Mercury. Then, obj['x'] gives the x position of Mercury on Jan 1st, 2017 (defined by "2017-01-01" in the second line).
Putting all together, including animation via matplotlib animation tool, the complete code that generates the movie at the beginning of this article is attached to the end of this article. It computes the orbits of the 4 inner planets (Mercury, Venus, Earth, and Mars) plus the Moon in years 2019 and 2020. This period can be easily adjusted in the code.
Q&A:
What are real and what are unreal in the animation?
The animation has realistic and accurate positions and speeds corresponding to the date shown at the top-left corner. The sizes of the objects are not to scale, although their relative sizes are correct. The Moon is on top of Earth because of this reason.
Does the code calculate the orbits in real-time, or it is just animating the data queried from the internet?
It is a real-time calculation. The initial conditions are from the internet. The program needs the positions and speeds of the planets at some date in order to predict their positions in future time.
Why is the orbit of Mercury so ugly?
Well, it is. Unlike any other planets in the solar system, Mercury has high orbital eccentricity, meaning its orbit is not very round.
Can you make this animation 3D?
Yes, but 2D seems to be good enough for this demonstration. I suspect it is possible to make it 3D while still keep the code in 100 lines. If a lot of readers are interested, I can make an update in the future.
Why only 4 planets, not all 8 of them?
Because starting from Jupiter, the 5th planet, the orbital radius becomes very large (>5 times that of Earth). Including them in the animation will make the inner planets nearly invisible.
