rec.autos.simulators

  I downloaded and printed out a ton of information on rigid body dynamics and
inertia tensors for angular momentum calculations, and understand some of the
concepts.  For instance, angular velocity is usually not manipulated directly
or used at all in rigid body dynamic simulation, as by applying forces to the
body through the inertia tensor, angular momentum gets changed directly and
appropriately instead, something that is supposedly very difficult and
expensive to do correctly to angular velocities.  

  So, suppose I express my initial inertia tensor using the three seperate
polar moments of inertia of the sprung weight of the car (assuming it's the
same as a 3-D box) in matrix form like this:

 Inertia1    0       0
    0     Inertia2   0
    0         0    Inertia3

  I downloaded a paper written by David Baraff, from the Robotics Institute at
Carnegie Mellon University, available at Chris Hecker's web page:
http://www.racesimcentral.net/

  This explained to me everything I currently know on this subject up until now
(which isn't much, obviously!).  In it, he shows how a force vector and a
position vector for the location of the force on the body can be used to find
the total torque vector on the body (for 3 dimensions).  However, I'm confused
as to what to do with this torque and the inertia tensor.  I'm using a 3x3
matrix to define the rotational orientation of the sprung mass of the car.
What do I do with the torque vector?  It says that the inertia tensor itself
changes with orientation.  Should I be multiplying my orientation matrix by the
inertia tensor, then multiplying the result by the angular momentum matrix and
the time step in order to change the orientation matrix again?  Any ideas??  I
know GPL must do this for the body, as well as all four wheels independently,
because there's gyroscopic precession in the wheels, and this is the only way I
can imagine it might be done.  

  Help!!  
Thanks,
Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://www.racesimcentral.net/

First of all, for those not interested, I appologize for the longish
post.

You are just one small step away from doing it correctly.

There is much talk of the 6DOF models. These 6 degrees of freedom are
the position of the center of gravity (3 components) and then three
variables that describe the rotational orientation of the object
relative to its c.o.g. . As you say, and it's the easiest way to do it,
instead of using, say, three angles (say, pitch, yaw, roll) to describe
orientation, you are using a 3x3 matrix that is very well known to
anyone ever doing 3D graphics, however, usually in 3D graphics the
position vector is also included in the last row (or column) to form a
larger 4x4 matrix.

The 3x3 orientation matrix, let's name it O, is an orthogonal matrix,
meaning that if you exchange rows and columns and multiply it by the
original matrix )transpose it, O^T ), you should get the identity
matrix.

Now, the inertia tensor (let's name it J0 in the local coordinates of
the body) can also be represented as a 3x3 matrix. However, this tensor
is a constant only in the system that rotates with the vehicle. The
tensor that you created with components in pitch, yaw, and roll axes is
an example of such a tensor, although generally you would also get
components away from the diagonal.

If you've done 3D graphics, you know that the position of an arbitrary
vector v0 that is given in the local coordinates of the body needs to be
multiplied by the orientation matrix O, v = O v0 , to get the vector v
as the vector v0 oriented in the world coordinates. The tensor J0 in
world coordinates, J, however transforms a bit differently (as a
matrix). Without going into too much detail, to obtain the world
representation of the tensor J, you need to perform the following
operation:

J = O J0 O^T, (1)

where O^T is the transposed of O as given before.

I guess you are familiar with the concept of rotational momentum, gamma,
which is a vector. The relation

d(gamma)/dt = torque, (2)

is true, where torque is the torque vector. Now, gamma is given as

gamma = J omega, (3)

where omega is the rotational velocity vector, which points along the
axis of the rotation and whose magnitue gives the speed of rotation.
Actually,

(d O)/dt y0 = omega (x) (O y0) (4)

holds true, where (d O)/dt is the time derivative of the orientation
matrix, (x) denotes the vectorial product, and v is an arbitrary vector.
Similarly,

(d O^T)/dt y = - O^T (omega (x) y) . (5)

Now one needs to take into account that both the tensor of inertia in
the world (but not local) coordinates and omega change with time.
Inserting (3) into (2), and taking into account both (4), (5) and (1),
one gets

omega (x) O J0 O^T omega - O J0 O^T (omega (x) omega)
+ O J0 O^T (d omega)/dt = torque (6)

where the second term clearly vanishes, and from this one obtains

(d omega)/dt = J^(-1) (torque - omega (x) (J omega) ). (7)

While the derivation was done in an inertial frame of reference (world
coordinates for example), the equation (7) can be rewriten in any frame
of reference, as long as all the quantities are transformed correctly.
In the local coordinate system of the body, the J is just J0 (and J^(-1)
is the inverse of the J0 (J0^(-1)), while in the world system

J=O J0 O^T,  J^(-1)= O (J0^(-1)) O^T.

If you need help on how to update the matrix O from a given vector
omega, let me know, but basically you need to construct a matrix that
creates a small rotation around the vector omega and multiply the matrix
O by it from the left.

-Gregor

P.S. There might be some small errors above as I've last visited this
derivation a while ago. The result (7) that you need, however, is surely
correct.

I must say I would heartily recommend this approach over what I
suggested. Working with both the angular momentum and angular velocity
is easier to implement than just the angular velocity approach (which is
what my derivation focused on). Blame it on my narrowminded education
:).

-Gregor

  Geeze, Gregor!  The paper I looked at just skipped through that derivation
because it was so complicated!  You know your stuff!  

Thanks a ton for your help, I'll probably be needing more to get this working
right. :-)

Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com

  Ok, let me see if I follow you here.  I was going to have a symmetrical
inertia tensor:

inertia1    0         0
  0       inertia2    0
  0           0     inertia3

  But what I really want is the inverse of this.  Is this it?  Silly question,
I know.

Ibodyinv:
inertia3    0         0
  0       inertia2    0
  0           0     inertia1

  I'm a real amateur with matrices and vectors at this point, as I've just
started using these one or two months ago.  My 3-D graphics are a wireframe
Win32Api thing that hardly uses vectors or matrices at all, so I'm still
learning this stuff.

  The angular momentum vector is a 3x3 matrix, correct?  Or is it a three
dimensional vector?

  The angular velocity vector is also a 3x3 matrix, correct?

  Ok, to calculate the angular velocity vector, I multiply the current rotation
matrix by the inverse of the inertia tensor, which is calculated one time at
the beginning of the simulation.  The, I multiply the result of that by the
transpose of the current rotation matrix.  What is the transpose of the
rotation matrix?  
Rotation matrix is:
xx  yx  zx
xy  yy  zy
xz  yz  zz

   Ok, after the result of this is found, I get the new angular velocity vector
by multiplying the angular momentum vector by Ibodyinv (the inverse of the
inertia tensor).  

  Now, the angular velocity vector is essentially added to the old rotation
matrix to get the new rotation matrix?

  Geeze, I guess I have a lot of homework to do here!  I'll go through this
paper again with your post as a guide.  Hopefully I'll get this working.
Studying vectors and matrices more ought to help.  

  Thank you very much for your help Petri.

Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com

Not quite. :) This is what I meant:

1/inertia1   0               0
0               1/inertia2   0
0               0               1/inertia3

In general, inverting a matrix is more complex, but if every element off the
diagonal is zero you can take this handy shortcut.

I'm using OpenGL for my own simulation. I must say, if I had to pick one
thing about OpenGL that I like the most, it's the ease with which it can be
learnt and used. I've gone in a couple of months from spinning, unlit
single-colored cubes to multitextured, lit, environment-mapped 3d models. I
don't think I could've done that with D3D. Other people's mileage may vary,
of course. :)

It's a 3-d vector.

Nope, it's also a 3-d vector. :)

Transposing a matrix means simply mirroring it along the diagonal, in this
case the transpose is

xx  xy  xz
yx  yy  yz
zx  zy  zz

Not by Ibodyinv, but by the result of what we discussed above (which is
Ibodyinv transformed to world coordinates).

No, you first form a matrix from the elements of the angular velocity vector
like this:

 0   -z    y
 z    0   -x
-y    x    0

Then you multiply this by the old rotation matrix, and add the result of
this to the old rotation matrix to get the new one.

If you're wondering where that matrix in the middle comes from, I believe
the formulation is in the Baraff paper you downloaded.

Then there's also the fact that you'll need to normalize the rotation matrix
every now and then to keep it orthogonal. Normalizing quaternions is faster,
which is one good reason to use them.

Yep, that's the key. And once you have the rigid body dynamics working,
adding things like aerodynamic effects is a breeze, if you'll pardon the
pun. :)

Happy to help.

By the way, you mentioned something about precession of the wheels in GPL? I
haven't seen that game (shame on me :), I assume you mean the wheels precess
if they go flying through the air after a crash? Sounds cool.

Petri Blomqvist

  Yes, basically.  If you get some air off a steep hill, you're supposedly able
to rotate the car by steering the front wheels once airborne.  GPL seems to be
the most complete model out there, but I haven't tried N4, Nascar Heat, or
Viper Racing, so maybe it's just the hype.  I don't imagine many folks would go
through the trouble of modelling each wheel with an inertia tensor, but GPL
seems to have done this.  I know I don't, it's complicated enough as it is
right now!  Might as well do the flywheel and engine this way too.  That way,
high engine rpm's would cause the chassis to have a little more resistance to
turning.  Probably not enough to notice, of course.  Don't count on me doing
this!  lol  Real frame rate killer for an unnoticable effect.

  Thanks for the help again, Petri.  I'll see what I can do with this and come
back later if I get stuck.

Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com

Quick and dirty inverter of a 3x3 matric for you, plenty of room for
improvement, but as you only need to call once (or whenever you mass
distribution changes) then it will do :O)

Chris

BOOL InverseMassMatrix(MATRIX3X3 *source, MATRIX3X3 *dest)
{
    long col,row;
    float fDet;
    float fInvdet;

    dest->matrix[0][0] = source->matrix[1][1] * source->matrix[2][2] -
source->matrix[1][2] * source->matrix[2][1];
    dest->matrix[0][1] = source->matrix[0][2] * source->matrix[2][1] -
source->matrix[0][1] * source->matrix[2][2];
    dest->matrix[0][2] = source->matrix[0][1] * source->matrix[1][2] -
source->matrix[0][2] * source->matrix[1][1];
    dest->matrix[1][0] = source->matrix[1][2] * source->matrix[2][0] -
source->matrix[1][0] * source->matrix[2][2];
    dest->matrix[1][1] = source->matrix[0][0] * source->matrix[2][2] -
source->matrix[0][2] * source->matrix[2][0];
    dest->matrix[1][2] = source->matrix[0][2] * source->matrix[1][0] -
source->matrix[0][0] * source->matrix[1][2];
    dest->matrix[2][0] = source->matrix[1][0] * source->matrix[2][1] -
source->matrix[1][1] * source->matrix[2][0];
    dest->matrix[2][1] = source->matrix[0][1] * source->matrix[2][0] -
source->matrix[0][0] * source->matrix[2][1];
    dest->matrix[2][2] = source->matrix[0][0] * source->matrix[1][1] -
source->matrix[0][1] * source->matrix[1][0];

    fDet = source->matrix[0][0] * dest->matrix[0][0] + source->matrix[0][1]
* dest->matrix[1][0] + source->matrix[0][2] * dest->matrix[2][0];
    if ( fabs(fDet) <= 1e-06 )
     return FALSE;

    fInvdet = 1.0f / fDet;
    for ( row = 0; row < 3; row++)
    {
        for ( col = 0; col < 3; col++)
            dest->matrix[row][col] *= fInvdet;
    }

    return TRUE;

whens that demo ever gonna be released.......sorry just had to ask

neil canning

source->matrix[0][1]

Ruud van Gaal, GPL Rank +53.25
Pencil art    : http://www.marketgraph.nl/gallery/
Car simulation: http://www.marketgraph.nl/gallery/racer/

Thanks again,
Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com

  Thanks Ruud.  I've got it working already.  My compiler has some built in
Matrix functions, so even though this is all greek to me, I can continue with
my project anyway :-)  

Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com

  It works!  Thanks for everybody's help here.  It seems to work properly.
Download a little demo here at
http://performancesimulations.com/files/vector3b.exe

  Input an initial torque in world space about the x,y, and z axes.  This
torque only lasts for one frame, then the object spins on it's own after that.
It looks like it is working correctly to me.  What do you guys/gals think?  It
may not work on either Windows 2000 or ME, I forget which one, but it's fine on
Win95/98.  Thanks for any input, and thanks again for the help.

Todd Wasson
---
Performance Simulations
Drag Racing and Top Speed Prediction
Software
http://PerformanceSimulations.Com