Solar system simulator

This repo contains a python script to simulate solar systems. The script has only been tested with our own Solar System.

This script was written in the summer of 2020 to improve on a program I had written for pico-8. Because of the pico-8's limited float and integer sizes (and the fact that I had to use Euler method for integration), the orbit calculations in that program were way off. This script is by no means an accurate representation of a planetary system -- it only calculates forces between each planet and the star -- but at least it seems to be fairly close to other simulations and measurements.

Requirements

The script requires the following modules:

• datetime for date calculations.
• spiceypy for initial planetary orbit positions.
• scipy for odeint() function to solve differential equations.
• numpy because math. And ndarrays.
• matplotlib to plot and animate the orbits.

How the script works

Enter the simulation start and end dates in the start_date_utc variables. The sec_per_iteration variable controls the precision of the calculations and the speed of the animation. The default value of 86 400 is equal to one day (60 sec/min * 60 min/hour * 24 hours/day), which means we perform one calculation of the orbits per day and that every frame of the animation represents one day.

Initialize a planetary system by doing something like Sun = planetary_system(10, "Sun", sun_mass_in_kg).

Initialize planets as Planet instances by doing something like earth = Planet(399, "Earth", Sun), where Sun is the planetary_system instance you constructed earlier. 399 is the NASA NAIF ID for the planet.

How accurate is the simulation?

I have not tested the accuracy of the simulation yet, but the orbits seem to be quite accurate (and Earth takes 1 year to orbit the Sun). The size of the planets and the Sun is not to scale. The sizes of the planets are correct relative to each other, but they are magnified relative to the Sun. The size of the planets and the Sun is also magnified relative to the distances between them.

Theory on planetary orbits

Planetary orbits are elliptical and orbit the system's barycenter.

The gravitational force on a planet with mass is given by:

Where is the gravitational constant, is the mass of the central object (the Sun in our case), is the mass of the planet and is the square of the distance between the objects.

We can also express the gravitational force as a vector:

where is the position vector and is the norm of such that .

From Newtons second law we have and substituting for gives:

We can find the velocity and position of the planet by integration as the acceleration of the planet is the derivative of the velocity and the double derivative of the position.

In this script we find the positions by numerical integration. To accomplish that we construct a vector such that:

This state vector contains both the position and velocity of the planet. We also define the function such that and use numerical integration to find the solution to our differential equation.

The initial values for our initial value problem are sourced from NASAs SPICE Toolkit.

Suggested reading on python simulation of the two-body problem:

• Scientific computing toolbox's introduction to odeint() in python
• C.P. Dullemond of the University of Heidelberg's introduction to the N-body model
• University of Oslo AST1100 lecture notes on the general solution to the two-body problem