# Modeling and simulation of planetary motion

One of my favourite problems in classical mechanics is the planetary motion, and one of my first attempts to solve it with computer algorithms dates from my early teens. My first approach was written in *BASIC*, and the results were completely wrong, as I didn’t know enough about the underlying maths. Subsequent attempts produced good results once I knew some basics of differential calculus during my high-school years.

In this post I’ll rescue one of those implementations, written in *C*, which uses the laws of motion to compute astronomical body paths.

## Physics background

Let’s suppose we have two bodies and , at positions and , and with masses and , respectively. Newton’s law of universal gravitation states that the force between the two bodies is proportional to the product of the two masses, and inversely proportional to the squared distance:

Where . For instance, the force with which is attracted towards , including magnitude and direction, and being the unitary vector from to , is

If we also use Newton’s second law of motion, which relates force and **acceleration**

and knowing that the acceleration is the second time derivative of **position**

we can express the problem with the following set of coupled non-linear **differential equations**:

If we provide the initial conditions, i.e. the initial positions and , and the initial velocities and , we can obtain, from the above equations, the paths and , regardless of the number of dimensions. Particularly, for a two-dimensional problem, where

the set of equations becomes

where we have to find , , and . The problem can be easily generalized to dimensions and bodies.

## Designing an numerical method to solve the equations

If we **discretize** the time domain

and we apply the following first order finite differences approximations of the differential operators (this is the easiest way, and it’s equivalent to Euler method):

We can compute both the velocity and position of a body , given the force exerted over it, by finding from the above equations

As simple as:

- We compute the forces for each body (the forces due to all the other bodies):
- Hence we have the acceleration over each body:
- We integrate it, once to obtain the velocity:
- And we integrate it again, for the position:

## The code

In the link below you have an implementation in *C* for the two-dimensional case. The code is ancient and platform-dependent, as it uses some basic graphic routines. Originally this code was written for MS-DOS, using the protected mode and the mode 13h. The code was originally compiled with the DJGPP programming platform under RHIDE. As this platform is from the past, you can either run it under an emulator like DOSBox, or compile it using the SDL wrapper that I’ve written for rescuing this and other old code (only tested under *GNU/Linux*). The wrapper is implemented in files `mode13h.h`

and `mode13h.c`

(the only platform-dependent code), and in order to switch between a native *MS-DOS* platform and *SDL* you must use the `DOS`

symbol (when present, it assumes a native *MS-DOS* platform, otherwise it will use *SDL*).

The main algorithm is implemented in file `newton.c`

, and, as I described before, it just computes the force exerted over each body, and then velocity and position.

A body is modelled with the struct `body_t`

, which considers the position and velocity. Before starting the algorithm, you must create the bodies with the `createBody()`

function, providing initial position, initial velocity and mass. Also, you can specify whether or not the body is fixed at its position, to simulate fictitious situations, suitable to being validated against Keppler’s laws (you can visually check how the two first laws hold when only two bodies interact). Another way to do this with real cases is by using the `bodyref`

variable: when this variable is not null, the camera will follow the selected body.

Other useful parameters are `scale_factor`

, to fit the required space in the screen, and `calculations_per_plot`

, that allows you to enable frame dropping, computing with enough accuracy without having to plot all the solutions.

The source code includes an example of an Earth–Moon system, where you can check that, given realistic data for the initial values and masses, some other parameters like the orbital period match the known values.

In this video you can watch an execution example of the three-body problem, under *DOSBox* emulator, where you can appreciate the chaotic behaviour.