Moving Planets - Newton's Laws Simulation in Python (2D)

In the 16th century, Tycho Brahe conducted extensive and unprecedentedly accurate astronomical observations. To achieve this, he had to invent new astronomical instruments.

In the first half of the 17th century, Johannes Kepler, using Tycho's observations, formulated three empirical laws, the most important of which states that the planets move around the Sun in elliptical orbits, with the Sun at one focus of the ellipse. This doesn't seem all that important and a big deal when you look at these ellipses on a schematic drawing of the Solar System. But when you consider that all the measurements of positions and velocities of the bodies moving through space were taken not from outside, but from one of the these bodies flying through space - Earth - you realize what a feat this was!

In the second half of the 17th century, Isaac Newton concluded that for planets to move in elliptical orbits, a force of attraction directly proportional to the masses of the planet and the Sun and inversely proportional to the square of the distance between them is sufficient. To achieve this generalization, he had to invent calculus (with another genius, Gottfried Wilhelm Leibniz) and to formulate the modern definition of force and the laws of mechanics.

In the first half of the 21st century I decided to explain this stuff to my grandchildren. This is a project I made with them a couple of years ago.

I wanted to show them in a visual form and to teach them some basic things about:

1. Newton's law of universal gravitation
2. Forces and Newton's second law of motion
3. Approximate calculation of derivatives and approximate integration
4. Object oriented programming - classes and inheritance

And I wanted to do it using the simplest tools - a programming language they already knew (Python ) and Turtle graphics to visualize the planets' motion.

I know that Python isn't the best language for teaching object-oriented programming and that it's not very efficient for tasks like simulation, but it's easy to write and run, and the boys already knew it and were using Turtle graphics at school (today one of them is already writing in C++).

And - last but not least - we had to do it in an hour or so, before two boys, aged 14 and 15, get bored and run away (okay, I didn't start the project from scratch).

Before turning it into a project for Instructables I had to clean it up and to write more comments.

Of course, this isn't a serious physics engine. It's a toy, but a very good one for those who want to see how mathematical formulas suddenly turn into realistic motion of physical bodies, play with the parameters of these bodies, feel how gravity attracts and accelerates bodies, and how the inertia of massive objects resists their acceleration... I understand that this is just how it should be, but it still feels like magic! I really enjoyed making this project, and I hope you enjoy it too.

If you want to skip the explanations you can proceed to step 10 to look at the examples of the simulation.

1. PyCharm IDE for developing and running Python;

1. Turtle Graphics:

1. abc module for defining Abstract Base Classes (see https://blog.teclado.com/python-abc-abstract-base-classes/):

All the definitions and functions we need for the simulation are concentrated in one module called NewtonCalc.py. We can launch simulation from small run modules with different configurations of physical body system.

NewtonCalc module is in fact a simple physical engine of the project. It consists of:

1. General definitions and functions such as setting screen size and zoom, functions converting radians to degrees and vice versa and so on;
2. Definitions of classes for representing physical values such as position, velocity, acceleration, forces and bodies;
3. The main loop of the simulation applying the forces to the bodies for the time quantum that has passed from the previous calculation and displaying the bodies on the screen.

I will provide link to the text of the module in the end of the Instructable in order to avoid too many technical details but I want to explain the classes I defined and the main features of the module.

Anyway, if you want to skip the explanation at all you can proceed to step 10 to look at the examples of the simulation.

Of course all the units and values in the code are purely fictional, I only wanted to show how the things work in general.

Object oriented programming in Python is a bit peculiar, you can read about it for example here: https://realpython.com/python3-object-oriented-programming/ .

Vector class is used for such things as position, velocity, acceleration and force.

For each vector I store both the value-direction (the angle from X axis) and width-height (X & Y components) pairs of variables - these is redundant but makes using the vectors simple.

Constructor: (__init__) in the terms of Python (it only resets the internal values):

"Set" methods:

A method for summing two vectors:

Force is an abstract class that contains abstract methods, i.e. only declarations of the methods that must be implemented in classes based on Force (inherited from Force). I could define forces in a more simple way but I wanted to teach the kids more concepts of object oriented programming.

I didn't start this project from scratch. I started it from another project we did with the kids, there we simulated tossing bodies in the field of gravity of the Earth. The forces we had in this previous project were the Earth gravity, drag force caused by the resistance of the air and drag force caused by wind. All of these forces describe interaction between a body and the environment. For example, though the gravity of the Earth is caused by the gravitational attraction between the Earth and a body, we don't consider the Earth as a physical object interacting with the body but use the free fall acceleration instead as an integral property of the environment.

In this new project we have another kind of force - the gravitational force between each pair of the bodies in the system. So I defined two types of forces in my code - "environmental" force and "interbody" force ('env' and 'inter'). Each specific force class inherited from Force class must implement a method of force calculation according to its type.

We don't need "environmental" forces for planet movement simulation but I decided not to erase them in order to make more general purpose (though primitive) "physical engine".

So Force abstract class has a constructor setting the type of the force:

And definitions of the methods of force calculations:

and

(I will describe class Body later though in Python it must be defined before it is used in Force).

Both methods must return a Vector describing the value and direction of the calculated force.

Gravitation is a class inherited from Force. Its calc_inter_force method implements the Newton's law of universal gravitation: the force of attraction between any two objects is directly proportional to the product of their masses and inversely proportional to the square of the distance between their centers. The law can be expressed by the formula F = G(m₁m₂/r²):

Other force classes are the environmental forces from the previous project: EarthGravity, Drag and Wind. They are not used for simulating planet movement, I only will write about them in a most general way:

1. EarthGravity is calculated as F = mg, where g is the Free fall acceleration (standard acceleration due to Earth's gravity): 9.81 m/s²
2. Drag is caused by the resistance of the air and is directly proportional to the surface of a body and to the square of its speed. The direction of the drag is opposite to the direction of the body's movement relative to the environment
3. Wind is similar to Drag, its direction is the same as the direction of the wind.

I'll show an example of simulation using these forces in Step 17 of this Instructable though it's not about space.

Body class contains the data of a physical body.

The parameters of the constructor are:

1. m - Mass
2. s - Surface (not used for planets)
3. x - X position
4. y - Y position
5. v - Velocity (vector)
6. shape - the shape of the "turtle" for drawing (I use "circle" for planets)
7. color - the color of the "turtle"

Besides the physical parameters Body contains the list of the forces applied to the body and their sum (the Net force):

In order to apply a force to body we must define it and add to the "Forces" list of the body, for example:

and append the body to the "Bodies" list of the system:

This is the heart of our "physical engine". The parameter of simulate function is a list of all the bodies of the system.

The function gets the current time and starts the main simulation loop:

Notice that the time in our simulation is the realtime. So the speed and the quality of the simulation depend on your computer.

In the end of each round I save the time when it began in prevTime variable and get the new time for the next round.

But what happens in between?

There is a loop for every Body of our system:

First of all we must calculate the Net force, i.e. the vector sum of all the forces applied to the body.

In a loop I calculate each force applied to the body according to the type of the force.

For an "environmental" force it's just one calculation:

For an "interbody" force we have to calculate the interaction of the current body with all the other bodies in the system:

Now we have calculated the Net forces applied to all the bodies of the system.

This is where all the Newton's magic happens.

We'll start a new loop for the bodies. Before this I calculate dt. In calculus, dt designates the differential of time, i.e. infinitely small change in time. Of course, in our case dt isn't infinitely small, it is just the smallest quantity of time we can measure - the time between the previous and the current calculation.

First of all we'll apply to each body the second Newton's law: we know the mass of the body and the Net force applied to it so we can calculate the acceleration: a = F/m.

All the calculations are performed separately for the X- and Y-components of all vector values:

Now we can start the integration.

Acceleration is the derivative of velocity with respect to time: a = dV/dt. So we can calculate the change in velocity dV as a * dt:

Velocity is the derivative of position with respect to time: V = ds/dt. So we can calculate the change in position ds as V * dt:

Now we have calculated the new values of velocity and position of the body. All we have to do is to store the new values in the body object and to move the body on the screen!

First of all, let's consider the simplest example: in our system there are two bodies with equal masses, and their initial velocities are zero.

In this case, the two bodies simply attract each other and eventually meet in the middle. We don't simulate collisions, so the bodies simply continue moving after they meet. Furthermore, as the distance between the bodies approaches zero, the simulation's result is unpredictable, so sometimes the bodies are "shot" and quickly "fly" away from the meeting point.

Here is the code to run:

In this example let's take two bodies with different masses (100000 and 50000 "mass units").

We can see that the less massive body (the blue one) is accelerated faster.

But what if the bodies have some initial velocity in the direction perpendicular to the line connecting them? They begin to move in elliptical orbits around their common center of mass! This already looks like the movement of planets or other space objects. Let's look at the case when two bodies have equal masses and equal initial velocities in opposite directions - one "up" and the second one "down":

We can see how the "planets" move in elliptical orbits, how they accelerate when they approach one another and decelerate when they move away.

Not every set of physical parameters (masses, initial positions, and velocities of the system's bodies) will yield a steady-state picture. And our engine isn't good enough to simulate transient processes that lead to a steady state — errors accumulate too quickly.

Furthermore, the overall momentum of the system must be balanced, otherwise it will drift across the screen - see the example: the velocities are like in step 12 (equal values in opposite directions) but the mass of the green planet is twice as big as that of the blue one and this makes the system drift up the screen.

Therefore, we have to experiment with the parameters to achieve a good result. But this is what allows you not only to understand but also to feel physics!

Now let's take two bodies with different masses and different velocities: the green planet is three times as massive as the blue one but the blue one has three times grater initial velocity:

Finally, let's create a small solar system... One body will be significantly more massive than the others - this will be the Sun. We won't force it to remain at the center and not to move - its mass will do it for us! We'll define four bodies and assign their initial velocities - the Sun and three planets, the greater the distance from the Sun, the slower the planet's velocity:

1. Sun: mass = 1000000, X=0, initial velocity = 0
2. Planet1: mass = 100, X=250, initial velocity = 180 # X is the distance from the Sun
3. Planet2: mass = 100, X=500, initial velocity = 110
4. Planet3: mass = 100, X=750, initial velocity = 90

But what happens if we let the simulation run for a long time? Yes, our physics engine isn't precise enough, and errors accumulate quite quickly. But watching the planets move is still mesmerizing... And look at the screenshot I took after an hour or two! Isn't it beautiful?

I only moved the Sun back to its original position in Photoshop because it had drifted an inch or two during the simulation.

This is an example showing how to use environmental forces: EarthGravity, Drag and Wind. This example simulates tossing physical objects at an angle from the Earth surface.

1. EarthGravity is calculated as F = mg, where g is the Free fall acceleration (standard acceleration due to Earth's gravity) - 9.81 m/s²
2. Drag is caused by the resistance of the air and is directly proportional to the surface of a body and to the square of its speed. The direction of the drag is opposite to the direction of the body's movement relative to the environment
3. Wind is similar to Drag, its direction is the same as the direction of the wind (it blows horizontally).

There are three bodies launched from the surface of the Earth at the same velocity and at the same angle:

1. The blue body is tossed without taking into account air resistance, as in a vacuum.
2. The green one takes into account air resistance, it rises to a lower altitude and falls down at a shorter distance from the start.
3. The violet one is the same as the green one but I add to it one more force: wind blowing horizontally. It rises to the same altitude as the green body but falls closer to the start because of the wind.

There is one more feature: when UP arrow key is pressed, the violet body "opens parachute". Technically this means that I increase the surface of this body and the Drag does its job - it slows down the speed. I also change the shape of the body (the "turtle") to a custom shape looking like a parachute.

Defining custom turtle shapes: https://www.geeksforgeeks.org/python/how-to-create-custom-turtle-shapes-in-python/

Using keyboard commands: https://www.geeksforgeeks.org/python/python-turtle-graphics-keyboard-commands/

This is a directory on Google drive with all the sources of the project:

https://drive.google.com/drive/folders/1kXQbh5s6-IjsSm-EPZdgXwHIly9ygnPA?usp=sharing

Thank you for reading and watching!