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 .

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 E_{tot} = 0.25. A consise description of these units can be found here.

read *t*=0 snapshot from input

calculate energy*dt*= 0.001

initialize*t*_{end} = 1

**WHILE** *t* < *t*_{end} **DO**

**END WHILE**

dump final snapshot

calculate energy

initialize

Calculate *a*_{i}(*r*_{i})

Advance positions*r*_{i} → *r*_{i+1}

Calculate*a*_{i+1}(*r*_{i+1})

Advance velocities*v*_{i} → *v*_{i+1}

*t* += *dt*

nsteps++

**IF** (nsteps%10 = 0)

**ENDIF**

Advance positions

Calculate

Advance velocities

nsteps++

calculate energy

print diagnostics (energy, energy error)

print diagnostics (energy, energy error)

dump final snapshot