I have an object being rotated by forces...

(vec3) torque = <some force>

(vec3) angMomentum += torque * dt;

(vec3) angVelocity = invInertiaTensor * angMomentum;

(quat) qrotation.Normalize();

(quat) qspin = quat(0.0, angVelocity) * 0.5 * qrotation;

qrotation += spin * dt; // current orientation

... I need to enforce angle limits about the x, y and z axis. So a torque force causing a rotation about the z axis for example wouldnt be able to rotate the object more then -45 or +45 degrees. What would be the best way to do this?

Usually hard limits on angular displacements are specified as constraints which are handed off to your dynamics engine which will (hopefully) use some computational magic to come up with a set of impulses that when applied to your rigid bodies will satisfy all the specified constraints simultaneously. That's the general approach anyway. It's a complex solution that's difficult to implement but it is very robust.

Perhaps you don't need that level of sophistication though. An easy solution is to detect when the constraint is violated and simply compute a corrective impulse in isolation such that the angular velocity is instantaneously changed and the constraint is not violated further. The angular displacement is likely to drift and violate the constraint even with the corrective impulse; you'll need to project it back into the valid space as well.

If rotation about x and y is prohibited and rotation about z is limited to +45 to -45 then a really simple way of enforcing this would be...

1) Convert the orientation quaterion to an axis and angle pair (v, theta)

2) Set v.x and v.y to zero and clamp theta to [-45, 45]

3) Convert the adjusted (v, theta) back into a quaternion and use this for the final orientation.

4) Set the x and y component of your angular velocity to zero. If theta was previously found to be outside of the [-45, 45] limit in step 2 then set the z component of the angular velocity to zero as well.

This is super naïve, especially steps 1, 2 and 3 but it should work.

First note that you integration is not correct:

You are multiplying the angular momentum from the *end* of the timestep with the inertia tensor from the *beginning* of the timestep. The inertia is a function of the orientation. In formulas:

omega( t + dt ) = InvInertia( t + dt ) * L( t + dt )

You have:

omega( t + dt ) = InvInertia( t ) * L( t + dt )

The is the essentially the same as dropping the gyroscopic term when integrating angular velocities.

To enforce your limits you need to decompose your quaternion into orthonormal factors (google for swing-twist decomposition) and then enforce the limit on your axis.

Given your quaternion q = ( x, y, z, w ). You decompose into q = q_xy * q_z where q_z = ( 0, 0, z, w ) / sqrt( z^2 + w^2 ) and q_xy = q * conjugate( q_z ). Then you apply the limit and multiply by q_xy.

HTH,

-Dirk

I also recommend normalizing at the end when using the quaternion derivative:

(quat) qrotation.Normalize();

(quat) qspin = quat(0.0, angVelocity) * 0.5 * qrotation;

qrotation += spin * dt; // current orientation

Should be:

(quat) qspin = quat(0.0, angVelocity) * 0.5 * qrotation;

qrotation += spin * dt; // current orientation

(quat) qrotation.Normalize();