Simulating Rotation with Euler Angles, a Quaternion, or a Rotation Matrix

CS 482 Lecture, Dr. Lawlor

Rigid-body rotation is only a bit more complex than translation.

Linear Translation
Position P
   vec3, meters
Angular Orientation R (see last lecture)
   3 angles in radians, a quaternion, or a matrix
Velocity V=dP/dt
   vec3, m/s
Angular Velocity W = dR/dt
   3 rates in radians/sec, a quaternion derivative, a skew matrix
   a vec3, radians/sec 
Acceleration A=dV/dt
   vec3, m/s^2
Angular Acceleration α = dW/dt
   3 rates in radians/sec^2
Mass m
   float, Kg
Rotational inertia matrix I
   mat3, Kg*m^2
Linear Momentum p=mV
   vec3, N*s
Angular Momentum L=I W
   vec3, Kg*m^2/s
Force F=mA=dp/dt
   vec3, N
Torque T = R x F = dL/dt
   vec3, N*m

The standard trick here is to ignore angular acceleration.  Instead, sum up net torques in 3D, integrate (multiply by dt and sum) to get angular momentum, divide by rotational inertia (or multiply by the inverse inertia matrix) to get angular velocity, and use that to adjust the current angular orientation.  Note that everything's a vec3 until that last step.

If you're using Euler angles, you need to compute your roll rate, pitch rate, and yaw rate from your current roll, pitch, and yaw and intended angular velocity.  This takes a whole lot of trig (yuk!).

If you're using a quaternion, the simplest solution is an axis-angle approach, treating the angular velocity vector W as a rotation axis, and the length as the rotation rate:
     var axis=W.clone();  axis.normalize();
     BABYLON.Quaternion.RotationAxisToRef(axis, W.length()*lib.dt, dR);

There's also a handy
quaternion differentiation formula:
      dR/dt = 0.5 * W * R
Expanding the derivative into a finite difference, we have:
      R += dt * 0.5 * W * R
Since W is a vec3, we first scale it by dt*0.5.  We then need to multiply it by the quaternion orientation R, so we promote it to a vec4 by tacking on a w component of 0 and multiply.  Finally, we add the resulting product to R term-by-term.  In GLSL, this is simple other than the quaternion multiplication quaternion_multiply:
     vec4 R=...; // orientation quaternion
     vec3 W=...; // angular velocity vector (radians/sec)
     R+=quaternion_multiply(vec4(dt*0.5*W,0.0), R);

In BABYLON.js, there doesn't seem to be a way to convert from Vector3 to Quaternion, or add a quaternion in palce, so we need to work term by term:
        var D=W.t(0.5*lib.dt);
        R.copyFromFloats(dR.x+R.x, dR.y+R.y, dR.z+R.z, dR.w+R.w);

Axis-angle is simpler to understand, but the trig inside setFromAxisAngle makes it a bit slower than the quaternion differentiation approach.  In both cases, it's good practice to normalize the quaternion afterward to keep roundoff from changing the object scale.

If you're using a rotation matrix, you can add in a scaled copy of this "skew matrix", which is just the time derivative of the per-axis rotation matrix, and the exact analog of the "axis nudging" we derived by hand:
[ 0     -W.z   W.y ]
[ W.z 0 -W.x ]
[ -W.y W.x 0 ]
You do need to first rotate the skew matrix by the current rotation, since angular momentum applies in global coordinates.  You also need to orthonormalize, or the object will slowly grow in size. 

(See a more complete discussion of applying angular velocity here.)

Er, JavaScript or WebGL doesn't seem to be running, so basically all you're going to see is the bare code...