# Simulate The Solar System with Processing

It’s New Year’s eve and it’s time for horoscopes! Or maybe not. If you are a believer, or if you’re not, it’s still a good time to learn something more about celestial mechanics. We’ll make use of Processing, an open source programming language and environment for people who want to create images, animations, and interactions. You can download it from their website. You also need the Traer Physics library for Processing. It’s free and easy to install (just follow the instruction and be happy).

### The Basics: Newton’s Law

One of the greatest discovery of mankind is due to Isaac Newton. I can’t say if he took inspiration from the New Gospel (“whatever you bind on earth shall be bound in heaven“) or from a falling apple. The fact is that with a force that behaves like this:

$F = G_{N}\frac{mM}{r^{2}}$

one can calculate the trajectory (orbit) of the planets in the Solar System. If you don’t know it already, these trajectories belong in general to a class of curves called conic section. In particular, the orbit of the planets appear to be ellipses. I think it’s very important to remember that not all the objects that feel gravity follow a close path such as an ellipse: the bodies in the Solar System that did not follow such orbit has been ejected from our neighborhood long time ago (billions of years!) or crushed against one another. Only the ones with stable paths are still there, such as planets, moons and thebodies int the asteroid belt between Mars and Jupiter.

### Not so Basic: Long Range Interactions and Orbit Stability

The main problem with orbits is that gravity acts between every planet and the Sun, but also among the planets themselves. Our place is actually very peculiar: the large part of the mass present in the Solar System is actually located at the center, i.e. the Sun, whose position is fixed for all practical purposes. This is not true for example for binary stars. The overall effect is that the orbits are not exactly described by an ellipse. Some small corrections must be applied, but as long are they are tiny the trajectory appear stable. However, these small effects build up and after billions of years the orbit can become unstable. These tiny deviations from the perfect orbit can be actually measured and have been really used to conjecture the existence of Uranus, Neptune and Pluto. Astronomers were able to measure even tinier deviations in the orbit of Mercury. In this case the effect could not be explained by a never-seen-before planed (Vulcan) but to the precession of perihelion according to Einstein’s General Relativity.

### Set up things

With Processing/TraerPhysics one can create new particles, assign mass/position/momentum to them, define the interaction and then let the system evolve. This is actually very easy! Let’s have a look:

ParticleSystem physics;
Particle sun;
Particle p1; //earth
physics = new ParticleSystem( 0., 0. );
sun = physics.makeParticle( m_sun/m_earth, width/2, height/2, 0 );
// 20 * AU = width/2 // radius of the mini solar system within a screen of 800x800
float r_0 = ( width/70 );
physics.makeParticle(m_earth/m_earth, width/2 + r_0*r_earth, height/2, 0 );
p1.velocity().set( 0, sqrt(mu_solar_system_reduced/r_0), 0. );
physics.makeAttraction( sun, p1, G_N, r_0/100 );

Due to some constraints in the handling of floating point variables, it is better not to use Kilograms and meters, but normalize masses to the Earth mass (5.9 10^24 kg) and lengths to the Astronomical Unit (1AU = 1.496 10^11 m = average distance Sun-Earth). As for the initial velocity, it is possible to estimate the speed at the perihelion, completely directed along the y axis:

$v^{2} = GM_{sun} ( \frac{2}{r} - \frac{1}{a} )$

With r = distance from the Sun and a = major semi axis of the orbit. So, at the perihelion, for circular orbits ( r=a is rather good approximation):

$v = \sqrt{ \frac{GM_{sun}}{r}}$

At this point we can also add more planets, and interaction among them:

Particle p3; //jupiter
Particle p4; //saturn
p3 = physics.makeParticle( m_jupiter/m_earth, width/2 + r_0*r_jupiter , height/2, 0 );
p3.velocity().set( 0, sqrt(mu_solar_system_reduced/r_jupiter/r_0), 0. );

p4 = physics.makeParticle( m_saturn/m_earth, width/2 + r_0*r_saturn , height/2, 0 );
p4.velocity().set( 0, sqrt(mu_solar_system_reduced/r_saturn/r_0), 0. );
physics.makeAttraction( sun, p3, G_N, r_0/100 );
physics.makeAttraction( sun, p4, G_N, r_0/100 );
physics.makeAttraction( p3, p4, G_N, r_0/100 );

Now, enjoy the video on youtube (watch it full screen to see Venus, Earth and Mars better):

If you want to try it out, this is the complete code:

import traer.physics.*;
float AU = 1.496e11; // m
// 20 * AU = width/2 // radius of the mini solar system
float G_N = 6.67e-11;
float m_sun = 1.989e30;
float m_venus = 4.8e24;
float m_earth = 5.9e24;
float m_mars = 6e23;
float m_jupiter = 1.9e27;
float m_saturn = 5e26;
float m_uranus = 8.26e25;
float m_neptune = 1e26;
// standard grav parameter mu=GM
float mu_solar_system = G_N * m_sun;
float mu_solar_system_reduced = mu_solar_system / m_earth;
float r_venus = 0.7;
float r_earth = 1.; //AU
float r_mars = 1.4;
float r_jupiter = 5.;
float r_saturn = 9.;
float r_uranus = 20.;
float r_neptune = 29;
ParticleSystem physics;
Particle sun;
Particle venus; //venus
Particle earth; //earth
Particle mars; //mars
Particle jupiter; //jupiter
Particle saturn; //saturn
Particle uranus; //uranus
Particle neptune; //neptune
void setup()
{
size( 800, 800 );
smooth();

float r_0 = ( width/70 );
fill( 0 );
ellipseMode( CENTER );

println( "mu of the solar system " + mu_solar_system );
println( "speed of earth at perihelion = " + sqrt(mu_solar_system/AU) + " m/s");
println(" width = " + width + " r_0 = " + r_0 );

physics = new ParticleSystem( 0., 0. );
physics.setIntegrator( ParticleSystem.RUNGE_KUTTA );

sun = physics.makeParticle( m_sun/m_earth, width/2, height/2, 0 );
sun.makeFixed();

venus = physics.makeParticle( m_venus/m_earth, width/2 + r_0*r_venus, height/2, 0 );
venus.velocity().set( 0, sqrt(mu_solar_system_reduced/r_venus/r_0), 0. );

earth = physics.makeParticle( 1, width/2 + r_0*r_earth, height/2, 0 );
earth.velocity().set( 0, sqrt(mu_solar_system_reduced/r_0), 0. );

mars = physics.makeParticle( m_mars/m_earth, width/2 + r_0*r_mars , height/2, 0 );
mars.velocity().set( 0, sqrt(mu_solar_system_reduced/r_mars/r_0), 0. );

jupiter = physics.makeParticle( m_jupiter/m_earth, width/2 + r_0*r_jupiter , height/2, 0 );
jupiter.velocity().set( 0, sqrt(mu_solar_system_reduced/r_jupiter/r_0), 0. );

saturn = physics.makeParticle( m_saturn/m_earth, width/2 + r_0*r_saturn , height/2, 0 );
saturn.velocity().set( 0, sqrt(mu_solar_system_reduced/r_saturn/r_0), 0. );

uranus = physics.makeParticle( m_uranus/m_earth, width/2 + r_0*r_uranus , height/2, 0 );
uranus.velocity().set( 0, sqrt(mu_solar_system_reduced/r_uranus/r_0), 0. );

neptune = physics.makeParticle( m_neptune/m_earth, width/2 + r_0*r_neptune , height/2, 0 );
neptune.velocity().set( 0, sqrt(mu_solar_system_reduced/r_neptune/r_0), 0. );

println( "Start:" );
//println( "vx = " + earth.velocity().x() + " vy = " + earth.velocity().y() );

physics.makeAttraction( sun, earth, G_N, r_0/100 );
physics.makeAttraction( sun, mars, G_N, r_0/100 );
physics.makeAttraction( sun, jupiter, G_N, r_0/100 );
physics.makeAttraction( sun, saturn, G_N, r_0/100 );
physics.makeAttraction( sun, venus, G_N, r_0/100 );
physics.makeAttraction( sun, uranus, G_N, r_0/100 );
physics.makeAttraction( sun, neptune, G_N, r_0/100 );
physics.makeAttraction( earth, mars, G_N, r_0/100 );
physics.makeAttraction( earth, jupiter, G_N, r_0/100 );
physics.makeAttraction( earth, venus, G_N, r_0/100 );
physics.makeAttraction( mars, venus, G_N, r_0/100 );
physics.makeAttraction( mars, jupiter, G_N, r_0/100 );
physics.makeAttraction( jupiter, saturn, G_N, r_0/100 );

physics.makeAttraction( uranus, jupiter, G_N, r_0/100 );
physics.makeAttraction( uranus, saturn, G_N, r_0/100 );

physics.makeAttraction( neptune, jupiter, G_N, r_0/100 );
physics.makeAttraction( neptune, saturn, G_N, r_0/100 );
}
void mousePressed()
{

}
void mouseDragged()
{

}
void mouseReleased()
{

}
void draw()
{
//println( "(x, y) = ( " + earth.position().x() + ", " + earth.position().y() + ")" + " (vx, vy) = ( " + earth.velocity().x() + ", " + earth.velocity().y() + ")");

physics.tick(1000);

background( 0 );

color y = color( 255, 204, 0);
fill(y);
noStroke();
ellipse( sun.position().x(), sun.position().y(), 12, 12 );
color g = color( 0, 200, 100);
fill(g);
ellipse( venus.position().x(), venus.position().y(), 5, 5 );

color b = color( 0, 0, 255);
fill(b);
ellipse( earth.position().x(), earth.position().y(), 5, 5 );

color r = color( 255, 0, 0);
fill(r);
ellipse( mars.position().x(), mars.position().y(), 3, 3 );

color j = color( 255, 150, 100);
fill(j);
ellipse( jupiter.position().x(), jupiter.position().y(), 10, 10 );

color s = color( 100, 110, 100);
fill(s);
ellipse( saturn.position().x(), saturn.position().y(), 15, 4 );
color w = color( 200, 150, 100 );
fill(w);
ellipse( saturn.position().x(), saturn.position().y(), 8, 8 );

color b2 = color( 0, 200, 200);
fill(b2);
ellipse( uranus.position().x(), uranus.position().y(), 8, 8 );

color b3 = color( 0, 0, 250);
fill(b3);
ellipse( neptune.position().x(), neptune.position().y(), 8, 8 );

}