Consider adding gravitational softening

[illuhad]

In nbody simulations, often gravitational softening is included. This is a small value added to the distance in the denominator of the force computation: F = G*m_j*m_2 normalize(r_j - r) * 1/(|r_j -r|^2+epsilon)

This rounds off the singularity at |r_j - r| = 0, and is in fact also physically meaningful since in practice bodies are not point-like particles anyway, so in the real world there is never a singularity.

Numerically, this has two advantages:

• When two particles have a close encounter, avoid forming bound pairs orbiting each other at extremely close distances and with velocities which may not be physically meaningful; this is particularly important when adaptive time stepping is involved (which we don't do) since this phenomenon can force time steps to be extremely small
• Avoid energy errors exploding in a close encounter in case of fixed time stepping as in our case