Simulating the Motion of Celestial Bodies in Python

Athos [Link] (s2549387)

March 2025

1 Introduction
The study of planetary motion often requires experiments that span years because the natural
timescales of these systems cannot be simply accelerated. The aim of this project is to code a
simulation of our solar system, perform experiments on said solar system, and compare to real data
to see if there is any validity in our simulation. Although the code implemented is relatively simple,
large-scale simulations have the potential to play a crucial role in exploring astronomical systems
that would otherwise remain beyond our reach due to the constraints of time.
The main aims of the project were to: produce a visual, animated simulation of the Solar System
(up to Jupiter), find the periods of each planet about the Sun, check whether energy is conserved
across all integration methods, and finally attempt to find a launch velocity and angle such that a
probe flies by Mars.

2 Methods

2.1 Classes

The code is governed by two different classes. The Planet class houses a bodies intrinsic values
such as it’s mass, position, velocity etc. The other class is the Space class. This is what holds
each Planet in a list called [Link] and stores an entire history of the positions, velocities
and accelerations of the Planets.

2.2 Initial Parameters

All parameters were saved to a JSON file which is then called into the code during the initialization
of the Space class. The methods of integration used in calculating a planet’s next position and
velocity required for a timestep to also be provided as an initial condition. In this code the
timestep used was 0.001 Earth years. A value for G was also passed in using revised units in terms
of Astronomical Units, Earth masses and years.
Within the JSON file, basic information is also saved for each Planet such as its name, mass, initial
position, etc.

2.3 Motion Functions

For simulating the motion itself, three different methods of integration were used. The Beeman ,
Euler-Cromer and Euler-Direct. The focus of the project was on the Beeman and all experiments
were conducted using it but it is still useful to discuss and compare different ways of calculating
the planets motions and how energy is conserved in each case.

1
For the Beeman to work it is necessary to have access to a planet’s current, previous and next
acceleration. The ’Planet’ class was structured such that each individual body will store its current
and previous accelerations. The ’next’ accelerations is always just calculated, used, and then set as
the new ’current’ acceleration at each timestep so there is no need to save it in the ’Planet’ class.
The project also uses two other integration methods called Euler-Cromer and Euler-Direct.
When using Euler-Cromer, we update the velocity and then use this updated velocity to calculate
the new position of each planet. Euler-Direct uses the same equations but reverses the order so
that we find the new position with the current velocity of the planet and then update the velocity.

2.4 Simulating Motion

When it came to simulating the code, the method used in the project was to just create a list of
lists, called [Link] positions, for all positions of each planet. The list of positions of a planet
correspond to the index assigned to the planet index in the [Link] list. The planets were
plotted as spheres using an imported Circle class and their motion was also animated using the
FuncAnimation package.

2.5 Conservation of Energy Experiment

Similar to our [Link] positions list, a [Link] velocities list was also generated and ap-
pended to at each timestep and was a history of all velocities for each planet so the Kinetic Energy
of each planet at each timestep could be easily calculated. For the calculation of the potential
energy of each planet, I had to ensure that each planet pair was only getting calculated once per
timestep or we would get double the total potential energy. To do this, I wrote a double for loop
going through the list of [Link]. The first went through a range of (i, len([Link])
and the second went through (i + 1, len([Link]). This structure ensured each planet pair
was only counted once.

2.6 Periods of Planets

All positions of planets over the entire simulation can be accessed but as they saved as position
vectors, we can also access each individual x and y component. I considered a planet to have
completed a full orbit once it crosses the positive x-axis of the Sun. For this section we just check
the y-component at each timestep and if the previous y-coordinate was negative, and the current
one is positive or 0, then we record the timestep the ’crossing’ occurred. We then append these
times and by subtracting these times from the one before them, I can get the interval between
crossings which is our period. Some planets will orbit the Sun multiple times and to account for
that the calculated periods are appended to another list and the period that is printed by the
computer is the average of the saved periods for each planet.

2.7 Satellite To Mars

The last experiment conducted was an attempt to get a satellite to do a fly-by of Mars. Note that
for this experiment I considered a fly-by to be any distance within 1 million km of Mars. I defined
the probe’s starting position to be 1/1000AU to the right of the starting position of Earth. It’s
velocity and launch angle was defined arbitrarily and was tuned to get the probe to Mars. When
the function corresponding to this experiment is called, it first creates the probe using the ’Planet’
class. It then runs the Beeman algorithm now including the satellite. It then checks the distance
between the satellite and Mars for each timestep and saves them to a list. The minimum value of
this list is then printed as the closest approach to Mars along with the time that this distance is
recorded. Now ideally we would want the probe to return to Earth so we also want to record the
distance between the Earth and Probe after it has reached Mars. To do this we can define a value

2
min index at which the closest approach to Mars is achieved and only check distances to Earth and
their time after this index.

3 Results

3.1 Simulation Results

When the motion of the bodies are calculated using the Beeman and the Euler-Cromer methods,
the motion is clear, smooth, circular and is stable over time. We can add some code to show a line
of the planet’s motion around the Sun for both methods and these are shown below.

Figure 1: Beeman Motion

3
 Figure 2: Cromer Motion

However when we use ”Euler-Direct” the system is clearly unstable. Planet’s slowly spiral out and
drift away from the Sun. Again this can be visualized and is shown below.

Figure 3: Direct Motion

4
3.2 Simulation Discussion

Overall the simulation behaves as it should for the Beeman and Euler-Cromer methods. The orbits
are stable over time, they all move in the same direction, etc. And despite it not being visible on
this scale, the Sun is also slowly drifting over time due to the collective ’tug’ of the Planets.
We see however that simulated motion for the Euler-Direct is not at all what we observe in real
life. This implies that the updated position always relies on the updated velocity and that changing
the order of operation does affect how the simulation will play out. This could point to some
underlying law that governs celestial mechanics and orbital motion.

3.3 Period Results

The simulation is run in time steps of 0.001 Earth years with the total number of iterations being
100000 such that the total run time is 100 years. The results obtained for the planets are as follows:

Mercury: 0.241 Earth years (87.965 days)

Venus: 0.615 Earth years (224.110 days)
Earth: 1.0 Earth years (364.635 days)
Mars: 1.881 Earth years (686.565 days)
Jupiter: 11.831 Earth years (4318.68 days)

Table 1: Table of Planets and their Periods

3.4 Period Discussion

Overall my periods agree with the real periods of the planets. There are small discrepancies here
and there but the biggest difference is in my period of Jupiter for which the real value is 4333 days.
I believe that these small difference arise from the fact that my simulation does not account for
the effects of Saturn, Neptune, and Uranus and maybe factoring in those planets would push my
results closer to the real values. Maybe a future improvement could be to also implement those
three planets and see how it affects the rest of the results.

3.5 Energy Results

For an easy comparison of the conservation of energy between the different integration methods we
can calculate the motion for each motion and plot the total energies all on the same graph. This is
shown below.
Here we can immediately see that when the motion is calculated using Euler-Direct, the total
energy is in fact not constant and confirms that the system would become unstable. However it not
immediately clear what the total energy looks like for the other two integration methods. However
if we use the built in zoom function for the graph we get a clearer picture.

5
Figure 4: Energy Comparison (Broad)

Figure 5: Energy Comparison (Zoom)

6
3.6 Energy Conservation Discussion

Thus we see that total energy is much better conserved when using the Beeman or the Euler-Cromer
methods. We also note that Beeman does not oscillate as much and is significantly more stable than
the other methods. The result derived from Euler-Direct is completely un-physical and the system
is gaining energy over time. One fact to note is that the Euler-Direct graph does not start at the
same point as the other two. This is quite interesting and not something I was expecting. I do not
know what is causing this issue so an improvement for the future could be to fix this.

3.7 Results for Satellite Mission

For the mission to Mars a satellite was generated with a mass of 700kg, an initial velocity of
≈ 10, 043m/s2 , a launch angle of θ = 3o (with respect to positive x-axis), and an initial position of
earth pos + (1/1000)AU. After the simulation is ran and the data is analyzed the following table
of results is printed for the trip to Mars.

Initial velocity taken as: 39853.87 m/s (Relative to Sun)

Closest distance to Mars achieved: 0.43177 million km
Total journey time being: 7.356 months

Table 2: Mars Trip Data

Now analyzing the data for the return to Earth, the following table is also printed.

Closest return to Earth is: 5.97 million km

Total return journey time being: 383.928 months

Table 3: Earth Return Data

3.8 Satellite Discussion

So we managed to get the satellite to get to Mars. The distance is within our 1 million km bound
and the journey time being 7 months is what we’d expect from this type of mission in the real
world. I did not manage to get the probe to return to Earth and so an improvement could be to
include Saturn, Neptune and Uranus and see how that affects the path of our satellite.

4 Conclusion
Overall I am pleased with the results. My orbits are stable and ’look’ like they are orbiting, the
energy of the system remains constant, the calculated periods are acceptable compared to the
observed values and the satellite mission to Mars was a success with a very close approach being
achieved. A few improvements could be made such as including the rest of the planets, giving
planets different starting y-coordinates, etc. One future endeavor with the project could be to try
and replicate more complicated systems such as binary star systems or even the 3-body-problem.