Rigid Body Dynamics

Rigid Body Rotation

 Intent

Although my main efforts were focused over about a week, for quite awhile I've wanted to develop a way to realistically move and rotate objects based on the physical properties of those objects, and be able to apply any number forces on an object in any direction at any point, and have the object respond in some semblance of reality. Also, it had to be realtime (not very hard with simple scenes) and be as general as possible (accomodating any shape of object, no matter how convoluted).

Probably many people have had some sort of first year physics course with a little rotational dynamics. None of the methods used in those courses can be applied to a general simulation, as they are confined to a few objects (cylinders, spheres, etc.), usually only 2-D, and have one or two forces acting on them. Even the moment of inertia, a critical quantity that is usually precomputed in textbooks (though they may mention an integral over volume) is in fact a 3x3 matrix of six unique values (called a second rank tensor). Those courses are a good intro, but the web proves to be a better resource than a general education in this (relatively arcane) case...

This is the result of my efforts. It's not that visually impressive, though it's nice to see the effects of not renormalizing.

These are two documents that I found invaluable in developing a rigid body dynamics simulation, and have nice equations and diagrams where this page is mostly text:

SIGGRAPH notes from 1997 in pdf. This file was most helpful (and I've not yet read the others). The intent of these is to present only the required material to an audience of non-mathematicians (familiarity with matrix operations is assumed, of course), and that makes it far superior to the large amounts of dense and unhelpful papers to be found on the web. It's main drawback, is the example code very annoyingly insists on converting parameters into and out of a single array (it is pretty old, so maybe there's some slight optimization from this or accomodation of old C standards)- but just ignore that, and ignore the collision detection/response material until you've got the non-contact dynamics down. It also fails to mention how to completely implement the ode solver function (which is actually quite simple if you're not too hung-up on accuracy).

Stephen Ehmann sums up much of those SIGGRAPH notes with this rigid body tutorial, with some psuedo-code superior to the c-code in the other notes.

I'm not trying to duplicate his tutorial, but rather to offer a couple of things I discovered while implementing it myself that is under-emphasized or ignored in those references.

 First Steps

I've already mentioned the moment of inertia tensor, and since it's a precomputed constant one of the first things to do is make a way to generate it (unless fancier 3d modelling programs can do it, but I'm just working with Wings3D) from a 3D object's surface data.

Here's a presently not very helpful page on doing just that, though I'll get around to making it better later. Actually it just tells how to find the approximate volume of an object, and the moment of inertia matrix and it's inverse is then computed from that and stored.

 ODE Solving

The SIGGRAPH notes leave it out entirely, but the other page manages to mention that the Rdot matrix can simply be added to the existing Rotation matrix just like the velocity can be added to the position in order to find the next position of an object. Both of those quantities actually need to be multiplied by the size of the time increment, but my simulation overlooks that for now (amounting to a time step of 1 unit).

 Quaternions and Matrix Drift
Both references suggest using quaternions, and there probably right, but I've not yet implemented them yet (and matrices are a bit more intuitive at the moment). The rotation matrix, if you recall, is somewhat redundant- they're unit length and they're orthogonal and therefore one row is simply the cross product of the other two. Quaternions distill the information down to four parameters instead of nine and therefore build up less error as the simulation evolves. Rotating the matrix in the conventional way with cos/sin pair matrices doesn't accumulate error very quickly, but you can see for yourself in my demo that strange things happen very rapidly with no correction, especially at higher torques. This can be corrected without quaternions at some additional expense (though quaternions themselves have to be renormalized and have some overhead be converted back and forth from matrices), by doing something like:
```	```
void matrix16f::Normalize(void)
{
vector3f right(	matrix[0],matrix[1],matrix[2]);
vector3f up(	matrix[4],matrix[5],matrix[6]);
vector3f out(	matrix[8],matrix[9],matrix[10]);

// make all axes orthogonal
out   = Cross(up,right);
right = Cross(out,up);

// normalize
right = right/right.Length();
up    = up/up.Length();
out   = out/out.Length();

matrix[0] = right.vertex[0]; matrix[4] = up.vertex[0]; matrix[8] = out.vertex[0];
matrix[1] = right.vertex[1]; matrix[5] = up.vertex[1]; matrix[9] = out.vertex[1];
matrix[2] = right.vertex[2]; matrix[6] = up.vertex[2]; matrix[10]= out.vertex[2];
}
```
```

For high torques, this may need to be done every time step.

One side effect of this correction is wobble, which should be evident in the demo, where even though a single axis is being rotated about, the object starts to drift along other axes. This is far superior to the skewing if no correction is used at all, in any case.