**Code:**'======================= RK4 ==========================================
'Use current positions to compute accelerations. Store them in arrays a1x(), a1y(), a1z().
'zero out arrays a1x(), a1y(), a1z()
For k = 1 To objects
a1x(k) = 0: a1y(k) = 0: a1z(k) = 0
Next k
For k = 1 To objects
Mu = ObjMass(k) * G
For j = k + 1 To objects
Mu2 = ObjMass(j) * G
dx = (objx(j) - objx(k))
dy = (objy(j) - objy(k))
dz = (objz(j) - objz(k))
D = Sqr(dx * dx + dy * dy + dz * dz)
If Mu > 0 Then 'ignore massless particles
a1 = -(1 / D) * (1 / D) * Mu
a1x(j) = a1x(j) + (dx / D) * a1
a1y(j) = a1y(j) + (dy / D) * a1
a1z(j) = a1z(j) + (dz / D) * a1
End If
If Mu2 > 0 Then 'ignore massless particles
a1 = -(1 / D) * (1 / D) * Mu2
a1x(k) = a1x(k) - (dx / D) * a1
a1y(k) = a1y(k) - (dy / D) * a1
a1z(k) = a1z(k) - (dz / D) * a1
End If
Next j
Next k
'Copy the current velocities into an arrays v1x(), v1y(), v1z().
For k = 1 To objects
v1x(k) = objvx(k)
v1y(k) = objvy(k)
v1z(k) = objvz(k)
Next k
'Use the velocities stored in v1x(), v1y(), v1z() to predict positions of objects one-half a step into the future.
'Store the positions in arrays p2x(), p2y(), p2z()."
halfstep = timescale / 2
For k = 1 To objects
p2x(k) = objx(k) + v1x(k) * halfstep
p2y(k) = objy(k) + v1y(k) * halfstep
p2z(k) = objz(k) + v1z(k) * halfstep
Next k
'Use the accelerations stored in a1x(), a1y(), a1z to predict future velocities at half a step into the future.
'Store the velocities in arrays v2x(), v2y(), v2z().
For k = 1 To objects
v2x(k) = objvx(k) + a1x(k) * halfstep
v2y(k) = objvy(k) + a1y(k) * halfstep
v2z(k) = objvz(k) + a1z(k) * halfstep
Next k
'Use the positions stored in arrays p2x(), p2y(), p2z() (these are the positions at half a step in future)
'to compute accelerations on the bodies a half timestep into the future.
'Store these accelerations in arrays a2x(), a2y(), a2z()
'zero out arrays a2x(), a2y(), a2z()
For k = 1 To objects
a2x(k) = 0: a2y(k) = 0: a2z(k) = 0
Next k
For k = 1 To objects
Mu = ObjMass(k) * G
For j = k + 1 To objects
Mu2 = ObjMass(j) * G
dx = (p2x(j) - p2x(k))
dy = (p2y(j) - p2y(k))
dz = (p2z(j) - p2z(k))
D = Sqr(dx * dx + dy * dy + dz * dz)
If Mu > 0 Then 'ignore massless particles
a2 = -(1 / D) * (1 / D) * Mu
a2x(j) = a2x(j) + (dx / D) * a2
a2y(j) = a2y(j) + (dy / D) * a2
a2z(j) = a2z(j) + (dz / D) * a2
End If
If Mu2 > 0 Then 'ignore massless particles
a2 = -(1 / D) * (1 / D) * Mu2
a2x(k) = a2x(k) - (dx / D) * a2
a2y(k) = a2y(k) - (dy / D) * a2
a2z(k) = a2z(k) - (dz / D) * a2
End If
Next j
Next k
'Use the accelerations stored in arrays a2x(), a2y(), a2z() to compute velocities at half a step into the future.
'Store these velocities in arrays v3x(), v3y(), v3z().
For k = 1 To objects
v3x(k) = objvx(k) + a2x(k) * halfstep
v3y(k) = objvy(k) + a2y(k) * halfstep
v3z(k) = objvz(k) + a2z(k) * halfstep
Next k
'Use the velocities stored in arrays v2x(), v2y(), v2z() to predict positions of objects one-half a step into the future.
'Store these positions in arrays p3x(), p3y(), p3z().
For k = 1 To objects
p3x(k) = objx(k) + v2x(k) * halfstep
p3y(k) = objy(k) + v2y(k) * halfstep
p3z(k) = objz(k) + v2z(k) * halfstep
Next k
'Use the positions in arrays p3x(), p3y(), p3z() (these positions are at half a step in future)
'to compute the accelerations on the bodies a half time step into the future.
'Store them in arrays a3x(), a3y(), a3z()
For k = 1 To objects
a3x(k) = 0: a3y(k) = 0: a3z(k) = 0
Next k
For k = 1 To objects
Mu = ObjMass(k) * G
For j = k + 1 To objects
Mu2 = ObjMass(j) * G
dx = (p3x(j) - p3x(k))
dy = (p3y(j) - p3y(k))
dz = (p3z(j) - p3z(k))
D = Sqr(dx * dx + dy * dy + dz * dz)
If Mu > 0 Then 'ignore massless particles
a3 = -(1 / D) * (1 / D) * Mu
a3x(j) = a3x(j) + (dx / D) * a3
a3y(j) = a3y(j) + (dy / D) * a3
a3z(j) = a3z(j) + (dz / D) * a3
End If
If Mu2 > 0 Then 'ignore massless particles
a3 = -(1 / D) * (1 / D) * Mu2
a3x(k) = a3x(k) - (dx / D) * a3
a3y(k) = a3y(k) - (dy / D) * a3
a3z(k) = a3z(k) - (dz / D) * a3
End If
Next j
Next k
'Use the accelerations in arrays a3x(), a3y(), a3z() to compute velocities at one FULL step into the future.
'Store them in arrays v4x(), v4y(), v4z().
For k = 1 To objects
v4x(k) = objvx(k) + a3x(k) * timescale
v4y(k) = objvy(k) + a3y(k) * timescale
v4z(k) = objvz(k) + a3z(k) * timescale
Next k
'Use the veolcities stored in arrays v3x(), v3y(), v3z() to compute the positions of objects one FULL step into the future.
'Store them in arrays p4x(), p4y(), p4z().
For k = 1 To objects
p4x(k) = objx(k) + v3x(k) * timescale
p4y(k) = objy(k) + v3y(k) * timescale
p4z(k) = objz(k) + v3z(k) * timescale
Next k
'Use the positions stored in arrays p4x(), p4y(), p4z()(these are the positions at one FULL step in future)
'to compute accelerations on the bodies in the future.
'Store these accelerations in arrays a4x(), a4y(), a4z()
For k = 1 To objects
a4x(k) = 0: a4y(k) = 0: a4z(k) = 0
Next k
For k = 1 To objects
Mu = ObjMass(k) * G
For j = k + 1 To objects
Mu2 = ObjMass(j) * G
dx = (p4x(j) - p4x(k))
dy = (p4y(j) - p4y(k))
dz = (p4z(j) - p4z(k))
D = Sqr(dx * dx + dy * dy + dz * dz)
If Mu > 0 Then 'ignore massless particles
a4 = -(1 / D) * (1 / D) * Mu
a4x(j) = a4x(j) + (dx / D) * a4
a4y(j) = a4y(j) + (dy / D) * a4
a4z(j) = a4z(j) + (dz / D) * a4
End If
If Mu2 > 0 Then 'ignore massless particles
a4 = -(1 / D) * (1 / D) * Mu2
a4x(k) = a4x(k) - (dx / D) * a4
a4y(k) = a4y(k) - (dy / D) * a4
a4z(k) = a4z(k) - (dz / D) * a4
End If
Next j
Next k
'Apply the Runge-Kutta 4th Order method to these arrays to determine the slope of the position and velocities
'Then add these slopes to the current positions and velocities to get the new positions and velocities.
'Velocity is the derivative of position.
'Acceleration is the derivative of velocity.
'The RK4 method states that the slope for velocity is
'(1/6) * timestep * (the acceleration at the beginning of the interval, stored in a1x(), a1y(), a1z()
' + 2 * the acceleration at the first estimate of the half time step, stored in a2x(), a2y(), a2z()
' + 2 * the acceleration at the second estimate of the half time step, stored in a3x(), a3y(), a3z()
' + the acceleration at the estimate of the end of a full time step, stored in a4x(), a4y(), a4z() ).
'The RK4 method states that the slope for position is
'(1/6) * timestep * (the velocity at the beginning of the interval, stored in v1x(), v1y(), v1z()
' + 2 * the velocity at the first estimate of the half time step, stored in v2x(), v2y(), v2z()
' + 2 * the velocity at the second estimate of the half time step, stored in v3x(), v3y(), v3z()
' + the velocity at the estimate of the end of a full time step, stored in v4x(), v4y(), v4z() ).
For k = 1 To objects
deltavx = (timescale / 6) * (a1x(k) + 2 * a2x(k) + 2 * a3x(k) + a4x(k))
deltavy = (timescale / 6) * (a1y(k) + 2 * a2y(k) + 2 * a3y(k) + a4y(k))
deltavz = (timescale / 6) * (a1z(k) + 2 * a2z(k) + 2 * a3z(k) + a4z(k))
objvx(k) = objvx(k) + deltavx
objvy(k) = objvy(k) + deltavy
objvz(k) = objvz(k) + deltavz
deltapx = (timescale / 6) * (v1x(k) + 2 * v2x(k) + 2 * v3x(k) + v4x(k))
deltapy = (timescale / 6) * (v1y(k) + 2 * v2y(k) + 2 * v3y(k) + v4y(k))
deltapz = (timescale / 6) * (v1z(k) + 2 * v2z(k) + 2 * v3z(k) + v4z(k))
objx(k) = objx(k) + deltapx
objy(k) = objy(k) + deltapy
objz(k) = objz(k) + deltapz
Next k