Resolving the n-t Coordinate System

Now that we have the final collision response, we need to figure out how to get the initial values for (v_ai)n, (v_ai)t, (v_bi)n, (v_bi)t and then when the problem is solved we have to convert the values in the n-t axes back into values in the x,y axes. Let's begin by first finding the vectors n and t.

To find n we want a vector that is unit length (length equal to 1.0) and along the line from the center of ball A(x_a0,y_a0) to the center of ball B(x_b0,y_b0). Let's begin by finding a vector from ball A to B, calling it N, and then normalizing it:

Equation 8: Computation of n and t

N = B – A = <xb0 – xa0, yb0 – ya0>

Normalizing N to find n, we get

n = N/|N| = <nx,ny>

Now we need the tangential axis t which is perpendicular to n. We could find it again using vector geometry, but there's a trick we can use; if we rotate n 90 degrees clockwise, that's what we want. However, when a 2D vector <x,y> is rotated in a plane 90 degrees clockwise, the rotated vector is just t=<-y, x>=<t_x, t_y>. Take a look at Figure 13.35.

Now that we have both n and t and they are both unit vectors (since n was unit, t is unit also) we're ready to go. We want to resolve the initial velocities of each ball into terms of n and t for ball A and B, respectively:

vai=<xvai, yvai>
vbi=<xvbi, yvbi>

This is nothing more than a dot product. To find (v_ai)n, the component of the initial velocity of ball A along axis n, here's what you do:

(vai)n = vai . n = <xvai, yvai> . <nx,ny>
       = (xvai*nx + yvai* ny).

Note the result is a scalar, as it should be. Computing the other initial velocities is the same as shown in Equation 9.

Equation 9: Components of v_ai along n and t:

(vai)n = vai . n = <xvai, yvai> . <nx,ny>
 (xvai*nx + yvai* ny)
=
(vai)t = vai . t = <xvai, yvai> . <tx,ty>
       = (xvai*tx + yvai* ty)

Components of v_bi along n and t:

(vbi)n = vbi . n = <xvbi, yvbi> . <nx,ny>
       = (xvbi*nx + yvbi* ny)

(vbi)t = vbi . t = <xvbi, yvbi> . <tx,ty>
       = (xvbi*tx + yvbi* ty)

Now we're ready to solve the problem completely. Here are the steps:

I'll leave step 4 up to you. Now let's talk about tensors. Just kidding, just kidding. Let's finish this bad boy off. At this point we have the final velocities:

Final velocity for Ball A in terms of n,t:

(vaf)nt=<(vaf)n,(vaf)t>

Final velocity for Ball B in terms of n,t:

(vbf)nt=<(vbf)n,(vbf)t>

Now, let's forget about collisions and think about vector geometry. Take a look at Figure 13.36; it illustrates the problem we have.

Stated in plain Vulcan, we have a vector in one coordinate system n-t that we want to resolve into x,y. But how? Again, we are going to use dot products. Take a look at the vector (v_af)_nt in Figure 13.36, forgetting about the n,t axes. We just want (v_af)_nt in terms of the axis x,y. To compute this we're going to need the contribution of (v_af)_nt along both the x- and y-axes. This can be found with dot products. All we need to do are the following dot products:

For Ball A, v_a in terms of n,t is

va = (vaf)n * n + (vaf)t * t
     (vaf)n * <nx,ny> + (vaf)t * <tx,ty>

Therefore writing with dot products,

xaf = <nx,0> . (vaf)n + <tx,0> . (vaf)t
    = nx*(vaf)n + tx*(vaf)t

yaf = <0,ny> . (vaf)n + <0,ty> . (vaf)t
    = ny*(vaf)n + ty*(vaf)t

For Ball B:

vb = (vbf)n * n + (vbf)t * t
     (vbf)n *<nx,ny> + (vbf)t * <tx,ty>

Therefore,

xbf = <nx,0> . (vbf)n + <tx,0> . (vbf)t
    = nx*(vbf)n + tx*(vbf)t

ybf = <0,ny> . (vbf)n + <0,ty> . (vbf)t
    = ny*(vbf)n + ty*(vbf)t

Send the balls off with the above velocities and you're done! Wow, that was a whopper, huh? Now, since I think that the code is a lot easier to understand than the math, I have listed the collision algorithm here from an upcoming demo:

void Collision_Response(void)
{
// this function does all the "real" physics to determine if there has
// been a collision between any ball and any other ball; if there is a
// collision, the function uses the mass of each ball along with the
// initial velocities to compute the resulting velocities
// from the book we know that in general
// va2 = (e+1)*mb*vb1+va1(ma - e*mb)/(ma+mb)
// vb2 = (e+1)*ma*va1+vb1(ma - e*mb)/(ma+mb)

// and the objects will have direction vectors co-linear to the normal
// of the point of collision, but since we are using spheres here as the
// objects, we know that the normal to the point of collision is just
// the vector from the centers of each object, thus the resulting
// velocity vector of each ball will be along this normal vector direction

// step 1: test each object against each other object and test for a
// collision; there are better ways to do this other than a double nested
// loop, but since there are a small number of objects this is fine;
// also we want to somewhat model if two or more balls hit simultaneously

for (int ball_a = 0; ball_a < NUM_BALLS; ball_a++)
     {
     for (int ball_b = ball_a+1; ball_b < NUM_BALLS; ball_b++)
         {
         if (ball_a == ball_b)
            continue;

         // compute the normal vector from a->b
         float nabx = (balls[ball_b].varsF[INDEX_X] –
                      balls[ball_a].varsF[INDEX_X] );
         float naby = (balls[ball_b].varsF[INDEX_Y] –
                      balls[ball_a].varsF[INDEX_Y] );
         float length = sqrt(nabx*nabx + naby*naby);

         // is there a collision?
         if (length <= 2.0*(BALL_RADIUS*.75))
            {
            // the balls have made contact, compute response

            // compute the response coordinate system axes
            // normalize normal vector
            nabx/=length;
            naby/=length;

            // compute the tangential vector perpendicular to normal,
            // simply rotate vector 90
            float tabx =  -naby;
            float taby =  nabx;

            // draw collision
            DDraw_Lock_Primary_Surface();

            // blue is normal
            Draw_Clip_Line(balls[ball_a].varsF[INDEX_X]+0.5,
               balls[ball_a].varsF[INDEX_Y]+0.5,
               balls[ball_a].varsF[INDEX_X]+20*nabx+0.5,
               balls[ball_a].varsF[INDEX_Y]+20*naby+0.5,
               252, primary_buffer, primary_lpitch);
            // yellow is tangential
            Draw_Clip_Line(balls[ball_a].varsF[INDEX_X]+0.5,
               balls[ball_a].varsF[INDEX_Y]+0.5,
               balls[ball_a].varsF[INDEX_X]+20*tabx+0.5,
               balls[ball_a].varsF[INDEX_Y]+20*taby+0.5,
               251, primary_buffer, primary_lpitch);

             DDraw_Unlock_Primary_Surface();

            // tangential is also normalized since
            // it's just a rotated normal vector

            // step 2: compute all the initial velocities
            // notation ball: (a,b) initial: i, final: f,
            // n: normal direction, t: tangential direction

            float vait = DOT_PRODUCT(balls[ball_a].varsF[INDEX_XV],
                                     balls[ball_a].varsF[INDEX_YV],
                                     tabx, taby);

            float vain = DOT_PRODUCT(balls[ball_a].varsF[INDEX_XV],
                                     balls[ball_a].varsF[INDEX_YV],
                                     nabx, naby);

            float vbit = DOT_PRODUCT(balls[ball_b].varsF[INDEX_XV],
                                     balls[ball_b].varsF[INDEX_YV],
                                     tabx, taby);

            float vbin = DOT_PRODUCT(balls[ball_b].varsF[INDEX_XV],
                                     balls[ball_b].varsF[INDEX_YV],
                                     nabx, naby);


            // now we have all the initial velocities
            // in terms of the n and t axes
            // step 3: compute final velocities after
            // collision, from book we have
            // note: all this code can be optimized, but I want you
// to see what's happening :)

            float ma = balls[ball_a].varsF[INDEX_MASS];
            float mb = balls[ball_b].varsF[INDEX_MASS];

            float vafn = (mb*vbin*(cof_E+1) + vain*(ma - cof_E*mb))
                          / (ma + mb);
            float vbfn = (ma*vain*(cof_E+1) - vbin*(ma - cof_E*mb))
                         / (ma + mb);

            // now luckily the tangential components
            // are the same before and after, so
            float vaft = vait;
            float vbft = vbit;
            // and that's that baby!
            // the velocity vectors are:
            // object a (vafn, vaft)
            // object b (vbfn, vbft)

            // the only problem is that we are in the wrong coordinate
            // system! we need to
           // translate back to the original x,y
            // coordinate system; basically we need to
            // compute the sum of the x components relative to
            // the n,t axes and the sum of
            // the y components relative to the n,t axis,
            // since n,t may both have x,y
            // components in the original x,y coordinate system

            float xfa = vafn*nabx + vaft*tabx;
            float yfa = vafn*naby + vaft*taby;

            float xfb = vbfn*nabx + vbft*tabx;
            float yfb = vbfn*naby + vbft*taby;

            // store results
            balls[ball_a].varsF[INDEX_XV] = xfa;
            balls[ball_a].varsF[INDEX_YV] = yfa;

            balls[ball_b].varsF[INDEX_XV] = xfb;
            balls[ball_b].varsF[INDEX_YV] = yfb;

            // update position
           balls[ball_a].varsF[INDEX_X]+=
                    balls[ball_a].varsF[INDEX_XV];
            balls[ball_a].varsF[INDEX_Y]+=
                    balls[ball_a].varsF[INDEX_YV];

            balls[ball_b].varsF[INDEX_X]+=
                    balls[ball_b].varsF[INDEX_XV];
            balls[ball_b].varsF[INDEX_Y]+=
                    balls[ball_b].varsF[INDEX_YV];

            }  // end if

         }  // end for ball2

     }  // end for ball1

}  // end Collision_Response

The code follows the algorithm almost identically. However, the code is from a demo that simulates a pool table system, so added loops test all collision pairs. Once inside the loop the code follows the math. To see the algorithm in action check out DEMO13_8.CPP|EXE (16-bit version, DEMO13_8_16B.CPP|EXE). A screen shot from it is shown in Figure 13.37. The demo starts up with a number of randomly moving balls and then the physics takes over. At the bottom of the screen the total kinetic energy is displayed. Try changing the coefficient of restitution with the right and left arrows keys and watch what happens. For values less than 1, the system loses energy, for values equal to 1, the system maintains energy, and for values greater than 1, the system gains energy—I wish my bank account did that!