Code: Particle System

Simulating a snowglobe, 50 particles with collision detection and response against a sphere with radius 1.5
User ability to shake the sphere left right up and down with corrected detection and response of particles.
Explanation:
First the position is updated using the Heun method, followed by velocity. Giving us our current state
Y(t) = (x(t),v(t)) and derivative d/dtY(t) = (v(t), F(t))
Notice mass is 1, and therefore trivial.
Since particles are not rigid bodies we assume no rotation, angular momentum, torque, inertia tensor etc...
We are assuming these spheres have only position and no orientation (rotation)

Following that collision detection is done by computing the distance to the center of the sphere and comparing against the radius.
Collision response is done first by updating position then updating velocity.
Position is updated by first finding the intersection of the line (created by the old position and the new position of the particle) and the sphere. A more detailed explanation of this implementation can be found here: http://paulbourke.net/geometry/sphereline/
Note: if b^2 -4*a*c is less than 0, we need to compute the normal of the center to the particle, and multiply it by the radius, adding it to the center, to give us a less accurate but close enough point on the sphere.
After updating the position on the sphere, and translating it by the distance the sphere has moved (from the user) the new velocity needs to be computed.
Finding the normal at the position and breaking down the velocity into tangential and normal components, it is as simple as flipping the normal and adding the two back together.



Note: debugging further would iron out small inaccuracies, such as too high velocity and other small cases
Also note that this is one of many many ways of implementing this. I supply this for your general knowledge, not as the optimal solution to this problem. This is just the first way I came up with.

using namespace std; 
double mytime=0.005; 
double NUMP =50; //Number of particles 


void onTimer(void *){

double howFast=0.005; //Note: Each frame is 0.005 seconds apart


if( DataManager::mUI->mControl->isPlaying()){
recordFrames();
DataManager::mUI->mControl->advanceFrame();
}


if(DataManager::mUI->mSim_but->value()){ // --------- begin simulation code-------------------
vector<Vec3d> newpartpos = DataManager::mParticlePos[DataManager::mParticlePos.size()-1]; vector<Vec3d> newpartvel = DataManager::mParticleVel[DataManager::mParticleVel.size()-1];

             for (int i=0;i<NUMP;i++){

                //store old values
double oldx = newpartpos[i][0];
double oldy = newpartpos[i][1];
double oldz = newpartpos[i][2];
//update position
newpartpos[i] = newpartpos[i] +
Vec3d(howFast*newpartvel[i][0],
 howFast*newpartvel[i][1],
 howFast*newpartvel[i][2]);//Heun Method
//update velocity
newpartvel[i] = Vec3d((newpartpos[i][0]-oldx)/howFast,
 (newpartpos[i][1]-oldy)/howFast,
 (newpartpos[i][2]-oldz)/howFast);
//compute distance to center
double distance2center = sqrt((newpartpos[i][0]-centerx)*(newpartpos[i][0]-centerx) +
 (newpartpos[i][1]-centery)*(newpartpos[i][1]-centery) +
 (newpartpos[i][2]-centerz)*(newpartpos[i][2])-centerz);

///Collision Detected
if (distance2center > 1.5){

double magc = sqrt((newpartpos[i][0]-centerx)*(newpartpos[i][0]-centerx) +
 (newpartpos[i][1]-centery)*(newpartpos[i][1]-centery) +
 (newpartpos[i][2]-centerz)*(newpartpos[i][2])-centerz);

Vec3d Nc = (newpartpos[i]-Vec3d(centerx,centery,centerz)) / magc;
//Line and Sphere intersection
double x1 = oldx; double x2 = newpartpos[i][0]; double x3 = centerx;
double y1 = oldy; double y2 = newpartpos[i][1]; double y3 = centery;
double z1 = oldz; double z2 = newpartpos[i][2]; double z3 = centerz;
double a = ((x2-x1)*(x2-x1))+((y2-y1)*(y2-y1))+((z2-z1)*(z2-z1));
double b = 2*(((x2-x1)*(x1-x3))+((y2-y1)*(y1-y3))+((z2-z1)*(z1-z3)));
double c = (x3*x3)+(y3*y3)+(z3*z3)+(x1*x1)+(y1*y1)+(z1*z1)-(2*((x3*x1)+(y3*y1)+(z3*z1)))-(1.5*1.5);

double u;
double inside = (b*b)-(4*a*c);

if(inside<0)
{
newpartpos[i] = Vec3d(centerx,centery,centerz) + (1.5)*(Nc) +
           Vec3d(centerx-oldcx,centery-oldcy,centerz-oldcz);
}

else{
double u1 = (-b + sqrt(inside))/(2*a);
double u2 = (-b - sqrt(inside))/(2*a);
if (abs(u2)>abs(u1)) u = u1;
else u = u2;

//snap ball to point on sphere
newpartpos[i] = Vec3d(oldx,oldy,oldz) + u *(newpartpos[i] - Vec3d(oldx,oldy,oldz)) +
           Vec3d(centerx-oldcx,centery-oldcy,centerz-oldcz);
}

//Cap velocity
svelox = (centerx-oldcx)/howFast;
if (svelox>5) svelox=5;
sveloy = (centery-oldcy)/howFast;
if (sveloy>5) sveloy=5;
sveloz = (centerz-oldcz)/howFast;
if (sveloz>5) sveloz=5;
//Compute magnitude & normal
double mag = sqrt((newpartpos[i][0]-centerx)*(newpartpos[i][0]-centerx) +
 (newpartpos[i][1]-centery)*(newpartpos[i][1]-centery) +
 (newpartpos[i][2]-centerz)*(newpartpos[i][2])-centerz);

Vec3d N = (newpartpos[i]-Vec3d(centerx,centery,centerz)) / mag;

//Break down vectors into normal and tangential
Vec3d Vn = (N[0]*(newpartvel[i][0]-svelox) +
    N[1]*(newpartvel[i][1]-sveloy) +
N[2]*(newpartvel[i][2]-sveloz))*N;

Vec3d Vt = newpartvel[i] - Vn;

//Update velocity after collision response
newpartvel[i] = Vt - (.7)*Vn;

}
}



int main(int argc, char ** argv){

// inits
// ---------------------------------------------
DataManager::mUI = new ParticleSystemUI;
// ---------------------------------------------

// testing only: init some points
vector<Vec3d> init_particles(NUMP, vl_0);
vector<Vec3d> init_velocity(NUMP, vl_0);

double i;
srand ( time(0) );
for (i=0;i<NUMP;i++){
init_particles[i] = Vec3d((rand()% 10-5.0)/5.5,(rand() % 10-5.0)/5.5,(rand() % 10-5.0)/5.5);
}
for (i=0;i<NUMP;i++){
init_velocity[i] = Vec3d((rand() % 10-5.0),(rand() % 10-5.0),(rand() % 10-5.0));
}

DataManager::mParticlePos.push_back(init_particles); // first frame of data
DataManager::mParticleVel.push_back(init_velocity);
// show UI
Fl::visual(FL_DOUBLE | FL_INDEX);
DataManager::mUI->Show();
Fl::add_timeout( 0.07, onTimer );
return Fl::run();
}