The mathematics of my simulation...

The mathematics involved are not super difficult, but do require some thought. Remember that I am only dealing with collisions that occur between two spheres. The approach that I took is below:

For the collisions to work both kinetic energy and momentum have to be conserved. To do this, the momentum equation was solved for the final velocity of the second object and that value was then plugged into the kinetic energy equation:

M1V1i + M2V2 i= M2V2f + M1V1f

(M1V1i + M2V2 i - M1V1f) / M2 = V2f

Letting W = M1iV1i + M2iV2 I, the new equation became:

M1V1i2 + M2V2i2 = M1V1f2 + (W - M1V1f) 2/M2

Foiling the added expression, the new equation became:

M1V1i2 + M2V2i2 = M1V1f2 + (W2 – 2WM1V1f + M12V1f2)/M2

Letting Z = M1V1i2 + M2V2i2 and splitting the numerator of the added expression the new equation became:

Z = M1V1f2 + W2/M2 – (2WM1V1f)/M2 + (M12V1f2)/M2

Next, the fact that this equation is a quadratic equation was used to find the possible roots for V2f. This was done by letting A = (M12/M2 +M1), B = (-2WM1)/M2, and letting C = (W2/M2 –Z):

V2f = (-B ± Ö(B2 –4AC))/(2A)

The positive root was assumed and since it is obvious that the final velocity could not be larger than the initial velocity of a sphere that has a mass larger than the second sphere that it is colliding with, if this was the case, the negative root was used. The value of V1f was calculated using the calculated V2f  value:

V1f =  (M1V1i + M2V2 i – M2V2f) / M1

Using this information and the angle of the collision between the two spheres, which came out to look like this:

q = arctan((Y1 – Y2) / (X1 – X2)) at the time of collision.

The program assigned new values to the Vx and Vy values of the two spheres:

Vx1= cos(q)V1f

Vx2= cos(q)V2f

Vy1= sin(q)V1f

Vy2= sin(q)V2f

To check if any two spheres had collided, the distance between them was compared to their total radii. Thus:

Collision occurred if Ö((X1-X2)2+ (Y1-Y2)2) >= R1 + R2