Modeling and simulation of planetary motion

I believe that the maximum exponent of curiosity is trying to understand how our own world works; thus, as a curious minded person, I’ve always been very fond of physics and applied sciences.

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 mathematics. Subsequent attempts yielded 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 that we have two bodies $B_1$ and $B_2$ , at positions $\vec{p}_1$ and $\vec{p}_2$ , and with masses $m_1$ and $m_2$ , respectively. The 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

$F = G\,\dfrac{m_1\,m_2}{r^2}.$

Where $r = \left|\vec{p}_1-\vec{p}_2\right|$ . For instance, the force with which $B_1$ is attracted towards $B_2$ , including magnitude and direction, and being $\hat{r}$ the unitary vector from $\vec{p}_1$ to $\vec{p}_2$ , is

$\vec{F}_{12} = G\,\dfrac{m_1\,m_2}{\left|\vec{p}_1-\vec{p}_2\right|^2}\hat{r} = G\,\dfrac{m_1\,m_2}{\left|\vec{p}_1-\vec{p}_2\right|^3}(\vec{p}_2-\vec{p}_1)$

Solution example when two bodies interact.

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

$\vec{F} = m\,\vec{a},$

and knowing that the acceleration $\vec{a}(t)$ is the second time derivative of position $\vec{p}(t)$

$\vec{a}(t) = \dfrac{d^2}{dt^2}\vec{p}(t),$

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

$\begin{array}{c} \dfrac{d^2}{dt^2}\vec{p}_1(t) = G\,\dfrac{m_2}{\left|\vec{p}_1(t)-\vec{p}_2(t)\right|^3}(\vec{p}_2(t)-\vec{p}_1(t))\\ \dfrac{d^2}{dt^2}\vec{p}_2(t) = G\,\dfrac{m_1}{\left|\vec{p}_1(t)-\vec{p}_2(t)\right|^3}(\vec{p}_1(t)-\vec{p}_2(t)). \end{array}$

Same solution as before, but showing the relative motion. The result is an ellipse.

If we provide the initial conditions, i.e. the initial positions $\vec{p}_1(t_0)$ and $\vec{p}_2(t_0)$ , and the initial velocities $\vec{v}_1(t_0)$ and $\vec{v}_2(t_0)$ , we can obtain, from the above equations, the paths $\vec{p}_1(t)$ and $\vec{p}_2(t)$ , regardless of the number of dimensions. Particularly, for a two-dimensional problem, where

$\begin{array}{c} \vec{p}_1(t) = (x_1(t), y_1(t))\\ \vec{p}_2(t) = (x_2(t), y_2(t)), \end{array}$

the set of equations becomes

$\begin{array}{c} \dfrac{d^2}{dt^2}x_1(t) = G\,\dfrac{m_2}{((x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2)^\frac{3}{2}}(x_2(t)-x_1(t))\\ \dfrac{d^2}{dt^2}y_1(t) = G\,\dfrac{m_2}{((x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2)^\frac{3}{2}}(y_2(t)-y_1(t))\\ \dfrac{d^2}{dt^2}x_2(t) = G\,\dfrac{m_1}{((x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2)^\frac{3}{2}}(x_1(t)-x_2(t))\\ \dfrac{d^2}{dt^2}y_2(t) = G\,\dfrac{m_1}{((x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2)^\frac{3}{2}}(y_1(t)-y_2(t)), \end{array}$

where we have to find $x_1(t), y_1(t), x_2(t)$ and $y_2(t)$ . The problem can be easily spread to N dimensions and M bodies.

Designing an algorithm to solve the equations

If we discretize the time domain

$t_n = t_0 + n\,\Delta_t,\qquad n \in \mathbb{Z}^+$

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):

$\begin{array}{c} \vec{v}(t) = \dfrac{d}{dt}\vec{p}(t) \simeq \dfrac{\vec{p}(t+\Delta_t)-\vec{p}(t)}{\Delta_t}\\ \vec{a}(t) = \dfrac{d}{dt}\vec{v}(t) \simeq \dfrac{\vec{v}(t+\Delta_t)-\vec{v}(t)}{\Delta_t} \end{array}$

3D solution example

We can compute both the velocity and position of a body $B_i$ , given the force exerted over it, by finding $\vec{p}_i(t+\Delta_t)$ from the above equations

$\begin{array}{l} \vec{F}_i = \displaystyle\sum_{j \neq i}G\,\dfrac{m_i\,m_j}{\left|\vec{p}_i(t_n)-\vec{p}_j(t_n)\right|^3}(\vec{p}_j(t_n)-\vec{p}_i(t_n))\\ \vec{a}_i = \dfrac{\vec{F}_i}{m_i}\\ \vec{v}_i(t_n+\Delta_t) = \vec{v}_i(t_n) + \Delta_t\,\vec{a}_i\\ \vec{p}_i(t_n+\Delta_t) = \vec{p}_i(t_n) + \Delta_t\,\vec{v}_i(t_n) \end{array}$

As simple as that, we compute the forces for each body (the forces due to all the other bodies), hence we have the acceleration, and we integrate it twice to obtain both the velocity and the position.

The code

In the link below you have an implementation in C for the two-dimensional case. The code is 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, compiling 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 codes (only tested under GNU/Linux). The wrapper is coded in files mode13h.h and mode13h.c (the only platform-dependent code), and 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 coded in file newton.c, and, as I described before, it just computes the force exerted over each body, and then velocity and position.

Example of execution when we show the relative motion following one of the bodies. The result is an ellipse.

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 the video below you can watch an execution example of the three-body problem, under DOSBox emulator, where you can appreciate the chaotic behaviour.

Links

Source code

December 5th,2010 Uncategorizedby ibancg. | tags: C++, differential equations, finite difference method, gravitation, Keppler, mathematics, mechanics, modeling & simulation, Newton, physics | No Comments

A circle algorithm for ZX Spectrum

As like many people of my generation, my first steps in computer programming were with a ZX Spectrum. Following a natural order, Sinclair BASIC was the first language I dealt with, but quite often the hardware limitations forced me to move to the Z80 assembler language. A common practice was coding some time critical routines in assembler to avoid many of the bottlenecks, like some graphic routines.

In this post, I’ll explain, as an example, how to outperform the built-in circle drawing routine.

Algorithm explanation

The original built-in algorithm (command CIRCLE in BASIC) drew the circumference with an angular sweep, from 0 to 360 degrees, and although it was implemented in machine code, it looked quite slow.

It’s possible to make a more efficient implementation without using trigonometric functions nor square roots, and even without using multiplications. The following algorithm is equivalent to the now known mid-point algorithm (but with different error formulation), and it’s based on the implicit equation of the circumference (let’s assume the circumference is centered at (0,0))

$x^2 + y^2 - r^2 = 0$

To draw the contour, we establish $(x,y) = (r,0)$ as the starting point, and draw the first octant, with angles between 0 and 45 degrees. As below 45 degrees the y component goes faster than the x one, we’ll increment the y coordinate inside a loop, and correct the path decreasing the x component whenever we consider we are “far away” from the contour. To evaluate how far we are, let’s compute the squared error

$e(x,y) = x^2 + y^2 - r^2$

Let $\vec{p}_{k} = (x_{k},y_{k})$ be the coordinates of a pixel at a given step k; at each step we have to choose whether to move upwards or towards the upper-left pixel:

$\begin{array}{l} \vec{p}_{1,k+1} = \vec{p}_{1,k} + (0,1)\\ \vec{p}_{2,k+1} = \vec{p}_{2,k} + (-1,1) \end{array}$

Our decision will be based on a minimum squared error criterion; for the two options above we have the following expressions for the error (in our discrete space, we measure at the middle of the pixels, so the mathematical center is at the middle of the pixel (0,0))

$\begin{array}{ll} e_{1,k+1} = & e(\vec{p}_{1,k+1}) = x_{k}^2 + (y_{k}+1)^2 - r^2 = e_{k} + 1 + 2\,y_{k}\\ \\ e_{2,k+1} = & e(\vec{p}_{2,k+1}) = (x_{k}-1)^2 + (y_{k}+1)^2 - r^2\\ &= e_{1,k+1} + 1 - 2\,x_{k} = e_{k} + 1 + 2\,y_{k} + 1 - 2\,x_{k} \end{array}$

Notice that we can know the error at a given step if we know the error previous to it, with only a few additions. Whenever $|e_{1,k+1}| < |e_{2,k+1}|$ , we choose $\vec{p}_{1,k+1}$ as the next pixel, otherwise we take $\vec{p}_{2,k+1}$ . We can simplify the comparison in the following way, knowing that $e_{2,k+1}$ is always negative and $x_{k}$ and $y_{k}$ are integers

$\begin{array}{l} |e_{1,k+1}| < |e_{2,k+1}| \,\,\,\, \Rightarrow \,\,\,\, e_{1,k+1} < -e_{2,k+1} \,\,\,\, \Rightarrow \,\,\,\, e_{1,k+1} < -e_{1,k+1} - 1 + 2\,x_{k} \\ \,\,\,\, \Rightarrow \,\,\,\, 2\,e_{1,k+1} < 2\,x_{k} - 1 \,\,\,\, \Rightarrow \,\,\,\, e_{1,k+1} < x_{k} - \frac{1}{2} \,\,\,\, \Rightarrow \,\,\,\, e_{1,k+1} < x_{k} \end{array}$

The initial value for the error is $e_{0} = 0$ , as the first pixel $\vec{p}_{0} = (r,0)$ satisfies the circle equation. Observe the error evolution in the figure beside; the corrections produced by the x component decrements tend to keep the error around 0.

The algorithm can be coded in C as follows:

int x = r;
int y = 0;
int e = 0;

for (;;) {

drawPixel(xc + x, yc + y);

if (x <= y) {
break;
}

e += 2*y + 1;
y++;

if (e > x) {
e += 1 - 2*x;
x--;
}
}

Once we know how to draw the first octant, the rest of the circle can be drawn applying symmetry. Instead of drawing one single pixel, we draw 8 of them:

...
drawPixel(xc + x, yc + y);
drawPixel(xc + x, yc - y);
drawPixel(xc - x, yc - y);
drawPixel(xc - x, yc + y);
drawPixel(xc + y, yc + x);
drawPixel(xc + y, yc - x);
drawPixel(xc - y, yc - x);
drawPixel(xc - y, yc + x);
...