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

CS 482 Lecture, Dr. Lawlor

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

Linear Translation
 Rotation
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
or
   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
OK!
 huh?!

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();
     dR.setFromAxisAngle(axis, W.length()*lib.dt);
     R.multiplyQuaternions(dR,R);

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 quat_mult:
     vec4 R=...; // orientation quaternion
     vec3 W=...; // angular velocity vector (radians/sec)
     R+=quat_mult(vec4(dt*0.5*W,0.0), R);

In THREE.js, there doesn't seem to be a way to convert from Vector3 to Quaternion, or add quaternions, so we need to work term by term:
     var v=W.t(0.5*lib.dt); // scale as vec3
     dR.x=v.x; dR.y=v.y; dR.z=v.z; dR.w=0.0; // copy vec3 to quaternion
     dR.multiplyQuaternions(dR,R); // W*R
     R.x+=dR.x; R.y+=dR.y; R.z+=dR.z; R.w+=dR.w; // add quaternions

Axis-angle is simpler to understand, but the trig inside setFromAxisAngle makes it slower than this quaternion derivative 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, or it will stop working after a quarter turn.  You then need to orthonormalize, or the object will slowly grow in size. 

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