I noticed that the wrong mass is used in the bodyBodyInteraction in oclNbodyKernel.cl (and probably as wel in the CUDA version, haven’t checked).
The function is defined as: float4 bodyBodyInteraction(float4 ai, float4 bi, float4 bj, float softeningSquared)
r.x = bi.x - bj.x;
float s = bj.w * invDistCube;
where bj.w is intended to be the mass used for calculating the acceleration of bi caused by mass bj which is returned thru ai.
The function is called as follows from the gravitation() function:
accel = bodyBodyInteraction(accel, SX(i++), myPos, softeningSquared);
You will notice that bj corresponds with myPos. As a consequence, the mass of bi is used for calculating the acceleration of bi, which leads to rather surprising results if the masses used diverge: big masses are suddenly in a big hurry.
Therefore, the call should be:
accel = bodyBodyInteraction(accel, myPos, SX(i++), softeningSquared);
and the calculation of the distancevector r in bodyBodyInteraction() should be performed as:
r.x = bj.x - bi.x;
which corresponds to the direction of the attractive force.