What is gravity ?

Gravity is the force of attraction that governs the motion of all larger objects in the Universe, eg. moons, planets, stars and comets.
 

It was formulated by Newton as the law of general gravitation :

F = k * m * M / (r*r).

In words it states :

"Two bodies attract each other with a force that is proportional to the mass of each body and is inversely proportional to the square of their distance apart."
 

Calculating the acceleration

If more than two bodies are involved, then every body attracts every other body according to the law of general gravitation. Applying the rule of superposition the resulting force on one of them can be obtained by adding the force experienced from each of the other bodies.

Instead of calculating the force, the acceleration can be obtained directly using another law that Newton discovered :

a = F/m = k * M / (r*r).

where M is the mass of the other body, r is the distance between the bodies and k is the gravitational constant.

If the vector from m to M is (x1,x2,x3) then the distance is :

r = sqrt(x1*x1+x2*x2+x3*x3)

and the acceleration in each of the three directions are :

a1 = k * M * x1 / (r*r*r)
a2 = k * M * x2 / (r*r*r)
a3 = k * M * x3 / (r*r*r)

This calculation is repeated for each pair of bodies. If the universe contains N bodies this means N*(N-1)/2 pairs.
 

What is softening ?

When two bodies get very close the effect of gravity grows rapidly. The numerical methods used accumulate unacceptable large errors in this situation, and a 'trick' is therefore used to soften that effect.

Instead of (r*r*r) I use (r*r + epsilon)^(3/2) where the epsilon parameter is specified in the input file. Setting epsilon=0 will switch off softening.
 

Simulating movement

Now the combined effect from all the other bodies on a specific body has been calculated, giving us the acceleration of that body. Next step is to update the position and velocity of the body relative to a fixed coordinate system thus simulating the movement that has happened within one time slice. In mathematical terms the position, the velocity and the acceleration of a body are related by differential equations. A number of numerical methods are available for approximation of the solution to these equations. I chose the so called 'leap-frog' method, in each dimension the velocity and the position are updated by calculating :

velocity = velocity + timeslice*acceleration

position = position + timeslice*velocity

Although very simple this algorithm seems to conserve energy and momentum very well indeed.