Force calculation

The force for each particle is calculated with a direct N-Body scheme, where all forces are calculated accurately from each object to each other object. The equation for the force is (boldface parameters are Cartesian vectors, see Wikipedia for more details.)

The acceleration follows from ${\bf a}_i = {\bf F}_i/m_i$.


The equations of motion are integrated via the Verlet-leapfrog
predictor corrector integrator, see Wikipedia for more details.

The equation to update the position:

The equation to update the velocity:


The units of the integration are so called Heggie-Units, in which case M=G=1. Here

As a consequence, the total energy of the system is Etot = 0.25. A consise description of these units can be found here.

The Algorithm

The gravitational N-body problem
read t=0 snapshot from input
calculate energy dt= 0.001
initialize tend = 1
WHILE t < tend DO
Calculate ai(ri)
Advance positions riri+1
Calculate ai+1(ri+1)
Advance velocities vivi+1
t += dt
IF (nsteps%10 = 0)
calculate energy
print diagnostics (energy, energy error)
dump final snapshot