# 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 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();
BABYLON.Quaternion.RotationAxisToRef(axis, W.length()*lib.dt, dR);
dR.multiplyToRef(obj.rotationQuaternion,obj.rotationQuaternion);

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);
dR.copyFromFloats(D.x,D.y,D.z,0.0);
dR.multiplyToRef(R,dR);
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...

 // Angular velocity change // (e.g., rcs system) var rot=new vec3(0,0,0); if (lib.key['a']) rot.y-=1.0; if (lib.key['d']) rot.y+=1.0; if (lib.key['w']) rot.x+=1.0; if (lib.key['s']) rot.x-=1.0; if (lib.key['q']) rot.z+=1.0; if (lib.key['e']) rot.z-=1.0; // Angular acceleration, radians/s^2 var radsPerSec2=1.0; // Add new rotation rate to angular momentum W.pe(rot.t(radsPerSec2*lib.dt)); trace("W=" + W); // Euler angles by default: trace("Euler angles: "+o.box.rotation); // Press 't' for quaTernion mode if (lib.key_toggle['t']) { // quaternions if (!o.box.rotationQuaternion) o.box.rotationQuaternion=new BABYLON.Quaternion(); var R=o.box.rotationQuaternion; // This is the differential rotation to apply: var dR=new BABYLON.Quaternion(); if (!lib.key_toggle['r']) { // Quaternion differentiation formula: // R += 0.5 * dt * W * R trace("Quaternion derivative mode"); var D=W.t(0.5*lib.dt); dR.copyFromFloats(D.x,D.y,D.z,0.0); dR.multiplyToRef(R,dR); // Add quaternions R.copyFromFloats(dR.x+R.x, dR.y+R.y, dR.z+R.z, dR.w+R.w); } else { trace("Quaternion axis-angle mode"); // Axis-angle version: BABYLON.Quaternion.RotationAxisToRef(W.clone(),length(W)*lib.dt,dR); // W is in global coords, so left-multiply: dR.multiplyToRef(R,R); } R.normalize(); trace("Quaternion mode:"); trace("dR: "+dR); trace("R: "+R); } else { // We are in matrix mode o.box.rotationQuaternion=null; // Project W into local coordinates var Wl=new vec3(dot(m.X,W),dot(m.Y,W),dot(m.Z,W)); // Push the Z axis around by the X and Y rotations m.Z.pe(m.X.t(Wl.y*lib.dt)); m.Z.me(m.Y.t(Wl.x*lib.dt)); // Push the X and Y around by the Z rotation m.X.pe(m.Y.t(Wl.z*lib.dt)); m.Y.me(m.X.t(Wl.z*lib.dt)); // Orthonormalize m.X=normalize(cross(m.Y,m.Z)); m.Y=normalize(cross(m.Z,m.X)); m.Z=normalize(m.Z); // Output as euler angles (for BABYLON) o.box.rotation = BABYLON.Vector3.RotationFromAxis(m.X, m.Y, m.Z); trace("Matrix mode:"); trace(m.X); trace(m.Y); trace(m.Z); } // Rocket physics: o.P.pe(o.V.t(lib.dt)); if (o.P.z<0.0) { // rockets bounce, right? o.P.z=0.0; o.V.z=0.0; o.V.x*=1.0-lib.dt; o.V.y*=1.0-lib.dt; } var Fg=new vec3(0,0,-9.8*o.mass); // gravity force (in -z direction) // Engine physics: var Fe=new vec3(0.0,0.0,0.0); // newtons of engine force var drymass=1000.0; if (lib.key_toggle['z'] && o.mass>drymass) { trace("ROCKET ENGINE RUNNING"); var burn_rate=100.0; // kg/sec burn rate var dm=burn_rate*lib.dt; // kg propellant burned this step // Don't burn more propellant than exists... if (o.mass-dm