Image: Photo of Enceladus, Saturn's moon of ice.

My name is Michael Lozovsky. I'm an undergraduate student at The Department of Physics at Technion. You can contact me via e-mail: smichl(at)mail.technion.ac.il

The theme of my project is simulation of planetary system.

.

Some Quotation

Presently my soul grew stronger; hesitating then no longer,
But the fact is I was napping, and so gently you came rapping,
And so faintly you came tapping, tapping at my chamber door,
That I scarce was sure I heard you''-here I opened wide the door;-
Darkness there and nothing more.

Edgar Allan Poe, The Raven
Go To Top

Theory

Introduction

A planetary system can be simulated by the N-body problem. The interaction between every two bodies is gravitational, given by the formula:

Then are masses of the bodies, are locations and G is Newton's Constant.
The star is the center of the system and has large mass compared to the smaller bodies- planets. The result is the planets executes elliptic motion around a star. Kepler's laws of planetary motion are:

• The orbit of a planet is an ellipse with the Sun at one of the two foci.
• A line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time.
• The square of the orbital period of a planet is proportional to the cube of the semi-major axis of its orbit.

The problem is interaction between different planets in the same system. If the planets are far apart and small enough we can neglect the gravity and each planet will sweep out an ideal ellipse. But in reality this approximation doesn′t work perfectly and real tracks aren′t Keplerian.

In system with giant planets effects like ″migration″ and ″scattering″ can accrue. The giant planets can attract other ones and move them from their original track. The orbit of a moved planet can become unstable and the small planet can fall on the giant one or on a star or even scatter to infinity.

I ignore relativistic effects and assume the information form one planet comes to all the other bodies on the system immediately.

Interesting questions and expected result

If I take a long enough time scale, I expect most of the planets will fall on their stars and/or scatter. One of the things I can check out is this stability over time.

The question of stability may be checked by plotting a semi-major axis( a distance form a star) as a function of time.

Go To Top

Simulation

My setup

I want to create the simulation of N body problem and set the parameters of those bodies. The main body (star) will have huge mass and leaser bodies (planets) will be set near it.

The input of the algorithm shall be

• mass of a star
• number of planets (1- 15)
• mass, location(3D) of each planet

The algorithm will evaluate the motion of a star and planets and give me a location of all the N bodies in every moment. The motion of planets is guided by a gravity, and the gravity force as we know is a function of distance between every two bodies. And as we know force is proportional to acceleration, and acceleration is second derivative of the location on time. As you can understand this means the motion is guided by an ODE.
If we have two bodies, the problem can be solved analyticity and one can get full information about motion of the bodies. But real systems have multiple body motions, which means solving system of differential equations, and every equation linked to all the others. This cannot be solved analyticity.
But a problem which is impossible to solve analyticity can be solved by step-by-step numeric solution. I divide a continuous time into set of discrete moments then every moment is initial condition the following one. The motion in a given moment is described by a formula:

I evolve every body in time by time step dt. So I use:

In every moment I have to find a velocity and acceleration that will be used as initial conditions for the next moment. My idea is finding a force on each and every body from a formula:

This will give me N-1 vectors of force (supposing we have N massive bodies) which I can sum and get net force. Dividing by body mass will give me the acceleration. After finding acceleration I can find a new speed by:

And this procedure occurs to every body in the system in every moment. In this manner I get a vector of coordinates of every body at every moment. So I can visualize the motion.
The user sets initial conditions of masses and coordinates. The initial velocity is set to Keplerian velocity given by:

I create a .xyz file at every iteration. In Aviz I can open filelist.dat to run the output of a simulation.
Note: I used real sizes and masses of bodies, in MKS. So the time step has to be a huge one to see a real motion. This simulation does not check collusion and absorption, do not have stellar wind and friction. So it is not realistic, but close enough.

Go To Top

Instructions

The project is written to run on Phelafel, any other Linux server or personal computer with aviz installed. This instruction suppose you use Linux command line, if you use Windows PC, I assume you use Cygwin and WinSCP.

Installation and first run

If you using Linux computer:

• Place the PS.tar file you have downloaded to destination directory on a computer. For example your desktop.
• Open a terminal by right click on the desktop and choosing terminal.
• type in the terminal: tar -xvf PS.tar. This creates "Planet" directory on your desktop.
• type: cd Planet/
• type: g++ Prog.C. That will compile the program and create a.out file.
• type: ./a.out. That will run the program.

If you using remote Linux server via Windows PC

• Place the PS.tar file you have downloaded to destination directory on a server using WinSCP or any other program. For example root directory.
• Open Cygwin.
• Type in command line: startx
• If you are not connecting form the Technion network, you have to connect via t2 to tx server using ssh your_acount@t2.technion.ac.il -Y. If you are connected via wired network or TechSec ignore this line.
• In the window that pops up type: ssh your_acount@phelafel.technion.ac.il -Y(for example. you may use any server with Aviz).
• Extract the .tar file by typing: tar -xvf PS.tar. This creates "Planet" directory.
• type: cd Planet/
• type: g++ Prog.C. That will compile the program and create a.out file.
• type: ./a.out. That will run the program.

Using the program

The program allows you to chose between some planetary systems I've created (as Solar system or well knowen Kepler-11), or set your own system. Just follow the instructions.
running the programs creates .xyz output files which can be found in Planet/XYZout. Note that the time step has to be big (10000 etc.) to see real movement. Note that solar system is not good for visualization because of scale problems, so it is not recommend. To run the animation you have to open aviz:

• Type on previous window: aviz. This opens us a program to visualize the motion
• Go to: Style->Spheres.
• Go to: Quality->Final.
• Go to: File->Open->Open File List...
• Double click on XYZout and double click on filelist.dat
• In a window that pops click on Single step
• (optional)Go to: Settings->Hide Axis.
• Go to: Settings->Hide Contour.
• (optional)Go to: Elements->Atoms... and set you favourite colours and sizes. I recommend a small star and tiny planets
• Click: Cycle and watch the motion. (Note: If the motion too slow- you have chosen very small time step)
• You may click on Autosnap if you want to record the motion. Creating animation described here, at a bottom of the page.
• Note that the program outputs "table.txt" file which contains distances from a star as a function of time. This file can be imported to Exel or Matlab etc. to plot graphs. You can run the plotting by typing matlab and then open file "plotting.m" .

Note: My simulation is 2D in 3D space. I am ignoring any z input from user. The main reason of that is a fact what I have to determine a initial velocity automatically and z coordinate adds me a degree of freedom - the body can move theoretically at any direction on a sphere. In 2D a body can move only clockwise or counterclockwise. I decided it will move counterclockwise.
One can improve my code by removing the line `"R[3*i+2]=0; //REMOVE this if you what to go 3D "` from the source file. But then he has to think about ingenious way to determine a direction of the motion. You may think that I did is wrong, but as we know, the planets form from single protoplanetary disc, so their motion is mostly 2D an not 3D. The PS.tar file contains .cpp file if anyone want to run the program under Windows. And finally I have to say the motion is build from 500 frames, while the program creates 1000 files. Matlab graphs counts all the 1000.

Go To Top

Results

I to represent some of the simulation results. The yellow big body is a star and smaller ones are planets. The visualization does not distinguish between large and small planets.
First we are looking in Kepler-11 motion:

You can see in the phenomenon of "falling" on a star and flying away, this is called "scattering". As you can see, the most of the planets scatter and the inner ones scatter first. In more realistic simulation one should take in account the cross section of the planet and star, but in this one I did not study a collapse at all.
Now we look on next result:

Here we can see a motion of two different mass planes that almost antipodes to each other. Their motion looks like they do not interact each other. You may see the planets are slow- this is because I took large radii.
Now let us look at very interesting (fictional) system:

As you may see, all the planets start in the same radius but after a time, the motion become asymmetric. Notice the blue planet- I created it much heavier then all the others planets. It drags all the system after self and cause all the system to disperse.
Now let us look at the last animation:

This is simple system motion. You may see here that the outer planets are much more stable then inner one.

You may notice that in some cases the planet orbits "inside" a sun. This is one of a side effects of my simulation: all the bodies are "radius-less dots" and not balls. Effects of collision does not interest me in this simulation. In addition I have to say that while the Semi-major axes are up to scale, the radii of planets and star of course are not.

Now I can represent some of graphs of planetary dynamics. First we looking at solar system:

What we can see here is radius of a motion of every planet as a function of time. We see that first two planets (Venus and Mercury) "falls" on a sun until a moment then they suddenly fly away. This phenomenon called "scattering", and you had seen it before (in animation). Third planet(Earth) approaching a sun, but much less dramatic. All other planets seems to be stable.
Do not worry, scattering in this system is not real, because as we all know solar system is a stable one. This scattering is a result of numerical simulation. If I take the time step to be too small, we will not see a motion. But if we take a big one it gives us instability and almost chaotic results.

Now let″s look on Kepler-11 system:

As we see, the system is very dense in the beginning and all the planets scatter. But the interesting point here is a time of the scattering: the outer planets seem to be more stable than inner ones.

The next system that is interesting me is a fictional one - I set one big and some small planets in very close radii (1 AU). The result of this asymmetric setup is dispersion of planets:

The interesting thing in this case is the fact that although the planets had began in the same distance from a star, we can see they finally dispersed.

You may think All the system getting unstable in my simulation. It is not true. For example I crated a fictional system with 3 planets:
First on is 1 Earth mass in distance of 1.4AU, second one is 3 Earth mass in distance of 7.1AU and last one is 20 Earth mass in distance of 11.2AU.

You can see that a planet on 1.4AU falls on a star, while outer two are completely stable.
That means my simulation is unstable in cases of small axes (about 1AU) and stable on large one. Go To Top