>  collide(particles_now, forces_now, gravity);

>  step(particles_now, particles_stepped, forces_now, dt/2);
>  collide(particles_stepped, forces_A, gravity);

>  step(particles_now, particles_stepped, forces_A, dt/2);
>  collide(particles_stepped, forces_B, gravity);

>  step(particles_now, particles_stepped, forces_B, dt);
>  collide(particles_stepped, forces_C, gravity);

If I am guessing correctly, "particles_stepped" includes velocity 
information?  Well you need to keep 4 separate copies of the velocity 
information just like you did with the force information.

>  for (int i=0; i<size; i++) {
>   forces_final[i] =
>    (forces_now[i]+2*forces_A[i]+2*forces_B[i]+forces_C[i])/6.0f;

Then you need to calculate the weighted sum of the velocities in this manner 
too here, and use the weighted velocity in the final step call.

I find it easier if you keep two distinct data structures, one for the state 
information (position, velocity), and one for the derivatives 
(momentum,force or velocity,acceleration it doesn't matter for point 
masses).

Then it is clear that "collide" must return a "derivative of state" data 
structure which must be then used in the following step function.  At the 
end, you do the weighted average of the "derivative of state" data structure 
as a whole.

This allows you to add more stuff easily in future (eg orientation and 
angular momentum) without having to worry about the actual RK4 core code.

Post a reply to this message