Collision of 2 Spheres/Circles

Started by
7 comments, last by jckut10 16 years, 8 months ago
I have scoured the internet for help with these physics but to no avail. I seem to remember that I covered this in either maths or physics at school, but can't find anything. Any help with this would be great. I want to understand the physics and then implement it myself. I have a project I'm working on and I need to work out the resulting velocity (magnitude and direction) of 2 spheres/circles after collision. From this velocity, I need to convert back to x and y magnitudes. I am using C++ and OpenGL with Dev C++. I currently have a class that models a ball. It has the ball's x and y position (called x/y), the x and y movement (called movX/moxY) (the amount the x/y position is increased by each time the idle function is called. These values can be negative.) The radius (called rad) and the mass (called mass) are also stored. I will explain what I have so far, please correct my errors. First of all, I detect a collision when the distance between 2 balls is equal to or less than ( <= ) the sum of the radii. To find the magnitude of the velocity of a ball I use
sqrt(movX*movX + movY*movY)

To find the resulting velocity I use
float v1f = ((v1*(m1-m2))+(v2*(2*m2)))/(m1+m2) ;
float v2f = ((v1*(2*m1))+(v2*(m1-m2)))/(m1+m2) ;

To find the angle from y = 0 (in the positive x-direction) I use
float angle = atan2( movY, movX ) ;
if( movX < 0 && movY > 0 )
{ angle = PI - atan2( movY, movX ) ; }
if( movX < 0 && movY < 0 )
{ angle += PI ; }
if( movX > 0 && movY < 0 )
{ angle = (2*PI)/3 + ( PI/2 - atan2( movY, movX ) )  ; } //should be 2*PI - atan2(movY, movX) - but this doesn't work
return( angle ) ;

Having found the angles of each ball, I add the angle of atan( vy1, vx2 ) to the angle of the first ball, and atan( vy2, vx1 ) to the angle of the second ball - as suggested by the end of this page. I then try to find the x and y velocity magnitudes with: - new movX = final velocity * cos( angle ) new movY = final velocity * sin( angle ) Despite my best efforts, the resulting effect appears very unrealistic. In the initial collision (with the way I have it set up at the moment), the balls stick together for a moment before going their seperate ways. The full source code can be found here. More on my website (navigate Projects | GLUT-Circles).
GLUT Setup: -Dev C++ 4.9.9.2GPU: GeForce 7600GT (256mb)CPU: Core 2 Due 1.86RAM: 1Gig
Advertisement
You may need to change these two lines:

float v1f = ((v1*(m1-m2))+(v2*(2*m2)))/(m1+m2) ;float v2f = ((v1*(2*m1))+(v2*(m1-m2)))/(m1+m2) ;


to

float v1f = v1 + ((v1*(m1-m2))+(v2*(2*m2)))/(m1+m2) ;float v2f = v2 + ((v1*(2*m1))+(v2*(m1-m2)))/(m1+m2) ;


*

You seem to be computing the direction of your velocity separately from the magnitude. It's easier to do this all at once, using vectors. Here's a link:

http://www.wheatchex.com/projects/collisions/

This is a pseudocode version of the page above, with e = 1 for elastic collisions.

' aCoords, bCoords : locations of objects a and b' vAi, vBi : vectors representing initial velocities of A and B' vAf, vAf : final velocities of a and b' AtoB : unit vector from the center of a to b.    AtoB = Normalize(bCoords - aCoords)        c = DotProduct(AtoB, (vAi-vBi))    vAf = vAi + AtoB*(-2 * (bMass * c) / (aMass + bMass))    vBf = vBi + AtoB*(2 * (aMass * c) / (aMass + bMass))
Original post by jckut10
You may need to change these two lines:
float v1f = ((v1*(m1-m2))+(v2*(2*m2)))/(m1+m2) ;float v2f = ((v1*(2*m1))+(v2*(m1-m2)))/(m1+m2) ;

to
float v1f = v1 + ((v1*(m1-m2))+(v2*(2*m2)))/(m1+m2) ;float v2f = v2 + ((v1*(2*m1))+(v2*(m1-m2)))/(m1+m2) ;


Doing this sends the balls flying off at a very high speed...

I am not too sure how to impelent the Normalize and DotProduct functions mentioned.

Is AtoB simply the distance between the center of the balls? If so, I already have this calculated.
GLUT Setup: -Dev C++ 4.9.9.2GPU: GeForce 7600GT (256mb)CPU: Core 2 Due 1.86RAM: 1Gig
AtoB is a unit vector from the center of A to the center of B.

So if A was at (0, 0) and B was at (30, 0), AtoB would be (1, 0).
If A was at (10, 6) and B was at (16, 6), AtoB would still be (1, 0).
If A was at (0, 4) and B was at (0, 10) then AtoB would be (0, 1).
et cetera.

*

Do you have an understanding of vectors? They're essential for physics simulations, and especially collisions. Gamedev has a section on them:

http://www.gamedev.net/reference/list.asp?categoryid=28#260

Especially look over the sections on vectors and physics.

Vectors let you express a value as a magnitude and a direction. A 2D vector has two components (x, y) ; a 3D vector has 3 components (x, y, z).

Most programmers define a vector class:

Type VECTOR  ' 3D Vector datatypex as single ' Each component is a single precision floaty as singlez as singleEnd Type


Then you could reference the components declare a vector (v) and reference the components as v.x, v.y, v.z.
*

To take the magnitude of v:
Magnitude = SQR(v.x*v.x + v.y*v.y + v.z*v.z)

*

Adding vectors is simple, just add the components.
' c = a + bc.x = a.x + b.xc.y = a.y + b.y

Subtraction works exactly the same way.
*

The dot product of two vectors is a scalar value:
DotProduct = a.x*b.x + a.y*b.y + a.z*b.z


( The cross product of two vectors is a bit more complex, and returns a vector that is perpendicular to the original two. You won't need it for collisions. But it is detailed in the tutorials I linked to above. )
*

You can also scale a vector by a value, simply by multiplying each component by that value:
 ' b = v*kb.x = v.x * kb.y = v.y * kb.z = v.z * k

*

To normalize a vector is to set its magnitude to exactly one. A vector with magnitude = 1 is a unit vector. Divide each component by the vector's magnitude:

' n = unit vector of vMagnitude = SQR(v.x*v.x + v.y*v.y + v.z*v.z)n.x = v.x / Magnituden.y = v.y / Magnituden.z = v.z / Magnitude

*

Velocity is a vector, speed is a scalar value -- so the values you have defined as velocities in your code are actually speeds.

Here's the useful vector functions:

Type VECTOR  ' 3D Vector datatypex as single ' Each component is a single precision floaty as singlez as singleEnd TypeFUNCTION vecNormalize (v as VECTOR) AS VECTOR	' n = vector of magnitude 1 pointing in the same direction as v	DIM AS SINGLE Magnitude = SQR(v.x*v.x + v.y*v.y + v.z*v.z)	DIM AS VECTOR n	n.x = v.x / Magnitude	n.y = v.y / Magnitude	n.z = v.z / Magnitude	RETURN nEND FUNCTION FUNCTION vecAdd (a as VECTOR, B AS VECTOR) AS VECTOR	' c = a + b	DIM AS VECTOR c	n.x = a.x + b.x	n.y = a.y + b.y	n.z = a.z + b.z	RETURN cEND FUNCTION FUNCTION vecSubtract (a as VECTOR, B AS VECTOR) AS VECTOR	' c = a - b	DIM AS VECTOR c	n.x = a.x - b.x	n.y = a.y - b.y	n.z = a.z - b.z	RETURN cEND FUNCTION FUNCTION vecScale (v as VECTOR, k as single) AS VECTOR	' n = v*k	DIM AS VECTOR n	n.x = v.x *k	n.y = v.y *k	n.z = v.z *k	RETURN nEND FUNCTIONFUNCTION vecDotProduct (a as Vector, b as vector) AS SINGLE	RETURN a.x*b.x + a.y*b.y + a.z*b.zEND FUNCTION


ANd here's the same code from the first post, using those functions. Last time I used overloaded operators, so it was hard to tell what was a vector and what was a scalar. Using the new functions, it should be easier to distinguish between the two.

' This is the code that should be executed when a collision is detected.
' aCoords, bCoords : locations of objects a and b, stored as position vectors
' aMass, bMass: masses of a and b
' aVelocity, bVelocity: vectors representing the velocities of A and B.
' AtoB: unit vector from a to b.

DIM AS VECTOR AtoB = vecNormalize(vecSubtract(bCoords, aCoords))

DIM AS SINGLE c = vecDotProduct(AtoB, vecSubtract(aVelocity,bVelocity))
' Calculate the new velocities
aVelocity = vecAdd(aVelocity, vecScale(AtoB, (-2 * (bMass*c) / (aMass + bMass)))
bVelocity = vecAdd(bVelocity, vecScale(AtoB, (2 * (aMass*c) / (aMass + bMass)))



Again, I'd highly recommend reviewing vectors, and how they relate to physics.
Yeah, I understand the concept of vectors, but couldn't remember how to maniplate them.

I have made my own Vector class and have implemented the code you suggest, and the result is very good so far, everything looks like it's real...

Thanks very much.

I need to test it a bit more now, but I will post the code up in a little while...
GLUT Setup: -Dev C++ 4.9.9.2GPU: GeForce 7600GT (256mb)CPU: Core 2 Due 1.86RAM: 1Gig
After further testing, it seem that when the mass of balls that collide are different, every time they collide, they speed up.

Note: I am assuming a fictionless enviroment.

This speeding up accelerates very quickly to the point of not being able to see the images as they are moving so fast...

Also watching the program progress while the masses are the same, it seems that the speed at which they travel increases, as before but at a much slow rate. Am I correct in thinking that, even in a frictionless environment, with perfect elasticy (coefficient of restitution = 1 ) objects can not gain energy? I realise it can be conserved, but this is an overall gain. Both objects are gaining energy...
GLUT Setup: -Dev C++ 4.9.9.2GPU: GeForce 7600GT (256mb)CPU: Core 2 Due 1.86RAM: 1Gig
Quote:Am I correct in thinking that, even in a frictionless environment, with perfect elasticy (coefficient of restitution = 1 ) objects can not gain energy?


You would be correct.

Hm, that is a very odd bug.

Maybe the program is continuing to modify the velocity after the collision has taken place. Any time the program calls the routine that does the post-collision math* the velocities of your objects will be altered as if a collision has occurred, even if a collision hasn't occurred. This might create the effect you describe.

I usually don't work with C, but maybe if you pot a bit of your code things will make more sense.

---

*"the routine that does the post-collision math" That would be this bit, shown below. In the last post I didn't collapse it into a subroutine, but this is what it might look like:

SUB DoCollisionMath (	BYVAL aMass as SINGLE, BYVAL aCoords as VECTOR, BYREF aVelocity as VECTOR , _			BYVAL bMass as SINGLE, BYVAL bCoords as VECTOR, BYREF bVelocity as VECTOR )' This is the routine that should be executed when collision has occured.' aCoords, bCoords : locations of objects a and b, stored as position vectors' aMass, bMass: masses of a and b' aVelocity, bVelocity: vectors representing the velocities of A and B.'	They are passed *by reference*, since they are changed in this sub.' AtoB: unit vector from a to b.	DIM AS VECTOR AtoB = vecNormalize(vecSubtract(bCoords, aCoords))    	DIM AS SINGLE c = vecDotProduct(AtoB, vecSubtract(aVelocity,bVelocity))		' Calculate the new velocities	aVelocity = vecAdd(aVelocity, vecScale(AtoB, (-2 * (bMass*c) / (aMass + bMass)))	bVelocity = vecAdd(bVelocity, vecScale(AtoB,  (2 * (aMass*c) / (aMass + bMass)))END SUB
Ok, well first of all, I have my Vector class. Here is the header file
class Vector{   private:      float xVal,             yVal,            zVal ;   public:      Vector() ;      Vector( float a, float b, float c ) ;      ~Vector() ;            float x() ;      float y() ;      float z() ;      void vecNormalize() ;      void vecAdd( Vector v ) ;      void vecSubtract( Vector v ) ;      void vecScale( float k ) ;      float vecDotProduct( Vector c ) ;            float vecDotProduct( Vector v1, Vector v2 ) ;      Vector vecAdd( Vector v1, Vector v2 ) ;      Vector vecSubtract( Vector v1, Vector v2 ) ;      Vector vecScale( Vector v, float k ) ;      Vector vecNormalize( Vector v ) ;} ;


And just now, having gone over the code, I have found the error...
This is the collision code -
//COLLOSION DETECTION HEREif( hypo <= sumR )           {                              Vector aCoords = Vector( balls.getX(), balls.getY(), 0 ) ;               Vector bCoords = Vector( balls[j].getX(), balls[j].getY(), 0 ) ;               Vector aVelocity = Vector( balls.getMovX(), balls.getMovY(), 0 ) ;               Vector bVelocity = Vector( balls[j].getMovX(), balls[j].getMovY(), 0 ) ;               Vector temp, temp2 ;                              //THIS WAS NOT COMMENTED OUT, CAUSING THE ERROR               //aVelocity.vecSubtract( bVelocity ) ;                              Vector AtoB = temp2.vecNormalize( temp.vecSubtract( bCoords, aCoords ) ) ;                   float c = temp2.vecDotProduct(AtoB, temp.vecSubtract(aVelocity,bVelocity)) ;                aVelocity = temp2.vecAdd( aVelocity, temp.vecScale( AtoB, (-2 * (m2*c) / (m1 + m2))) ) ;               bVelocity = temp2.vecAdd( bVelocity, temp.vecScale(AtoB,  (2 * (m1*c) / (m1 + m2))) ) ;                              balls.setMovX( aVelocity.x() ) ;               balls.setMovY( aVelocity.y() ) ;                              balls[j].setMovX( bVelocity.x() ) ;               balls[j].setMovY( bVelocity.y() ) ;}

I originally has the Vector class without the functions that return Vectors, so to add a vector to another, you had to make an instance of it, then add other vector to it. This would have been the reason why my first try at implementing the code failed. This all works now, and perfectly. Now I have to make an empty instance that returns a new vector. Hence temp and temp2. Their member variables are always 0, and never called on.

I guess I could use this for 3D collisions as well. Would it work of non-spherical objects? I'm getting ahead of myself...

Thanks very much for your help jckut10.
GLUT Setup: -Dev C++ 4.9.9.2GPU: GeForce 7600GT (256mb)CPU: Core 2 Due 1.86RAM: 1Gig
Glad I could help. :)

Yes, this will work in 3D as well, as long your vector functions (vecAdd, vecSubtract, vecNormalize, et cetera) include the z component.

*

Obviously, in "The Real World" non-spherical objects sometimes spin as a result of a collision. The changes in velocity due to a collision are actually the due to the force one object applies on another . The applied force might generate torque, which affects the angular acceleration of the object.

The applied force on a given object is easy to get:

F = m*a = m * (Vf - Vi)

Vf is your final velocity vector (after the collision); Vi is the velocity vector before the collision. m is mass of the object.

Then you can find the torque vector using

t = CrossProduct(r, F)

Where t is the torque vector, r is the vector from the center of gravity to the point on the object where the force is applied (which can be difficult to determine), and F is the applied force vector.

Then,

t = I*alpha
(... therefore ...) alpha = t / I

Where t is torque, I is the rotational inertia (determined by the shape, size, and density of the object) and alpha is the angular acceleration vector, which then can be added to the current angular velocity to get the new angular velocity.

As you can see, the concept is very straightforward, but the math to determine 'r' and 'I' can be difficult, so I don't bother calculating it.

This topic is closed to new replies.

Advertisement