After a lot of searching, I have finally found a reference that has what I needed. Although it doesn't appear in a 'Quaternion to Euler' or 'Quaternion to Tait-Bryan' google search, when I gave up and started looking for 'Quaternion to Axis-Angle' with the intention of going through that representation as an intermediate step, I came across the following wonderful document:
Technical Concepts
Orientation, Rotation, Velocity and Acceleration,
and the SRM
Version 2.0, 20 June 2008
Author: Paul Berner, PhD
Contributors: Ralph Toms, PhD, Kevin Trott, Farid Mamaghani, David Shen,
Craig Rollins, Edward Powell, PhD
http://www.sedris.org/wg8home/Documents/WG80485.pdf
It covers a lot of the formalisms, but most importantly, shows derivations and solutions for 3-1-3 and 3-2-1 Euler angle representation. It also seems to cover inter-conversion between pretty much every other rotation representation I'm aware of, and so I would also recommend it as a good general reference.
Oh, and the actual solution for a 3-2-1 ($z-y-x$) Tait-Bryan rotation convention from that reference:
$$
\phi = \operatorname{arctan2}\left(q_2 q_3 + q_0 q_1,\frac{1}{2}-(q_1^2 + q_2^2)\right) \\
\theta = \arcsin(-2(q_1 q_3 - q_0 q_2)) \\
\psi = \operatorname{arctan2}\left(q_1 q_2 + q_0 q_3,\frac{1}{2}-(q_2^2 + q_3^2)\right)
$$
Note that the gimbal-lock situation occurs when $2(q_1 q_3 + q_0 q_2) = \pm1$ (which gives a $\theta$ of $\pm \frac{\pi}{2} $), so it can be clearly identified before you attempt to evaluate $\phi$ and $\psi$.
(Convention for arctan2 is $\operatorname{arctan2}(y, x)$, as hinted on page 3.)
I think I have solved my problem. The main reference that finally clicked with me was this page, but I will continue and lay out my work here.
Here is a badly-drawn version of the coordinate system I'm using here (it looks better in my notes, I swear...)
bad plane drawing*
- $\psi$ represents a rotation around the $z$-axis, yaw
- $\theta$ represents a rotation around the $x$-axis, pitch
- $\phi$ represents a rotation around the $y$-axis, roll
The names of the angles are the traditional ones ($\psi$ for yaw, etc.) but the axes they are paired with are not.
Tait-Bryan rotations here will be performed in the order yaw, pitch, roll in the body-fixed reference frame. To apply them in the body-fixed frame, they must be applied in the reverse order: roll, pitch, yaw.
Angles to Quaternion
To convert from these three angles to a quaternion, I will use the definition of quaternion for each axis and then multiply them to get an overall quaternion.
$\begin{align}
Q_p &= \left(\cos\frac{\theta}{2}, \sin\frac{\theta}{2}\left(\begin{array}{ccc}1\\0\\0\end{array}\right)\right)\\
Q_r &= \left(\cos\frac{\phi}{2}, \sin\frac{\phi}{2}\left(\begin{array}{ccc}0\\1\\0\end{array}\right)\right)\\
Q_y &= \left(\cos\frac{\psi}{2}, \sin\frac{\psi}{2}\left(\begin{array}{ccc}0\\0\\1\end{array}\right)\right)
\end{align}$
Thus the overall quaternion I seek is $Q_R = Q_y\times Q_p\times Q_r$. Since the action of quaternions is applied right-to-left, this series of multiplications corresponds to the desired roll, pitch, yaw order.
Substituting and using the associative property:
$Q_R = \left(\cos\frac{\psi}{2} + \sin\frac{\psi}{2}\hat{k}\right)\times \left(\left(\cos\frac{\theta}{2} + \sin\frac{\theta}{2}\hat{i}\right)\times\left(\cos\frac{\phi}{2} + \sin\frac{\phi}{2}\hat{j}\right)\right)$
where $\hat{i}, \hat{j}, \hat{k}$ represent the three basis elements of the quaternions corresponding to rotations of $\pi/2$ about the $x$, $y$, and $z$ axes, respectively.
Performing the multiplication:
$\begin{align}
Q_R =& \left(\cos\frac{\psi}{2} + \sin\frac{\psi}{2}\hat{k}\right)\times\left(\cos\frac{\theta}{2}\cos\frac{\phi}{2} + \cos\frac{\theta}{2}\sin\frac{\phi}{2}\hat{j} + \sin\frac{\theta}{2}\cos\frac{\phi}{2}\hat{i} + \sin\frac{\theta}{2}\sin\frac{\phi}{2}\hat{i}\hat{j}\right) \\
=& \cos\frac{\psi}{2}\cos\frac{\theta}{2}\cos\frac{\phi}{2} +
\cos\frac{\psi}{2}\cos\frac{\theta}{2}\sin\frac{\phi}{2}\hat{j} +
\cos\frac{\psi}{2}\sin\frac{\theta}{2}\cos\frac{\phi}{2}\hat{i} +
\cos\frac{\psi}{2}\sin\frac{\theta}{2}\sin\frac{\phi}{2}\hat{i}\hat{j} +\\
& \sin\frac{\psi}{2}\cos\frac{\theta}{2}\cos\frac{\phi}{2}\hat{k} +
\sin\frac{\psi}{2}\cos\frac{\theta}{2}\sin\frac{\phi}{2}\hat{k}\hat{j} +
\sin\frac{\psi}{2}\sin\frac{\theta}{2}\cos\frac{\phi}{2}\hat{k}\hat{i} +
\sin\frac{\psi}{2}\sin\frac{\theta}{2}\sin\frac{\phi}{2}\hat{k}\hat{i}\hat{j}
\end{align}$
Using the standard results for multiplication of quaternion elements,
$\begin{array}{ccccc}
\times & 1 & \hat{i} & \hat{j} & \hat{k} \\
1 & 1 & \hat{i} & \hat{j} & \hat{k} \\
\hat{i} & \hat{i} & -1 & \hat{k} & -\hat{j} \\
\hat{j} & \hat{j} & -\hat{k} & -1 & \hat{i} \\
\hat{k} & \hat{k} & \hat{j} & -\hat{i} & -1
\end{array}$
these terms can be collected and placed in a column vector
$Q_R = \left(\begin{array}{c}w\\x\\y\\z\end{array}\right)
=\left(\begin{array}{c}
\cos\frac{\psi}{2}\cos\frac{\theta}{2}\cos\frac{\phi}{2} - \sin\frac{\psi}{2}\sin\frac{\theta}{2}\sin\frac{\phi}{2} \\
\cos\frac{\psi}{2}\sin\frac{\theta}{2}\cos\frac{\phi}{2} - \sin\frac{\psi}{2}\cos\frac{\theta}{2}\sin\frac{\phi}{2} \\
\cos\frac{\psi}{2}\cos\frac{\theta}{2}\sin\frac{\phi}{2} - \sin\frac{\psi}{2}\sin\frac{\theta}{2}\cos\frac{\phi}{2} \\
\cos\frac{\psi}{2}\sin\frac{\theta}{2}\sin\frac{\phi}{2} - \sin\frac{\psi}{2}\cos\frac{\theta}{2}\cos\frac{\phi}{2}
\end{array}\right)$
and this is my first desired result!
Quaternion to Angles
The next step is to get back from a quaternion $(w,x,y,z)^T$ to the angles $\psi, \theta, \phi$. This eventually condenses down to three simple formulae, but the derivation takes a few steps.
First, realize that the rotation described by the quaternion can also be described by the product of three rotation matrices
$M_R = M_y\times M_p\times M_r =
\left[\begin{array}{ccc}
\cos\psi & -\sin\psi & 0 \\
\sin\psi & \cos\psi & 0 \\
0 & 0 & 1
\end{array}\right]\times
\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos\theta & -\sin\theta \\
0 & \sin\theta & \cos\theta
\end{array}\right]\times
\left[\begin{array}{ccc}
\cos\phi & 0 & \sin\phi \\
0 & 1 & 0 \\
-\sin\phi & 0 & \cos\phi
\end{array}\right]$
$M_R = \left[\begin{array}{ccc}
\cos\psi\cos\phi-\sin\psi\sin\theta\sin\phi & -\sin\psi\cos\theta & \cos\psi\sin\phi+\sin\psi\sin\theta\cos\phi \\
\cos\psi\sin\theta\sin\phi+\sin\psi\cos\phi & \cos\psi\cos\theta & \sin\psi\sin\phi-\cos\psi\sin\theta\cos\phi \\
-\cos\theta\sin\phi & \sin\theta & \cos\theta\cos\phi
\end{array}\right]$
The elements of this matrix can be used to extract expressions for the angles:
$\frac{M_{R21}}{M_{R22}} = \frac{-\sin\psi\cos\theta}{\cos\psi\cos\theta} = -\tan\psi, M_{R32} = \sin\theta, \frac{M_{R31}}{M_{R33}} = \frac{-\cos\theta\sin\phi}{\cos\theta\cos\phi} = -\tan\phi$
A quaternion can also be represented as a rotation matrix (some elements omitted since they are unncessary for this derivation--full matrix here: H-T-T-P://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm):
$M_Q = \left[\begin{array}{ccc}
- & 2xy-2wz & - \\
- & 1-2x^2-2z^2 & - \\
2xz-2wy & 2yz+2wx & 1-2x^2-2y^2
\end{array}\right]$
Since, $M_R=M_Q$, the above expressions involving the elements of $M_R$ can be written with the elements of $M_Q$:
$\begin{align}
-\tan\psi &= \frac{2xy-2wz}{1-2x^2-2z^2} &\Rightarrow& \psi = \arctan\left(\frac{-2xy+2wz}{1-2x^2-2z^2}\right) \\
\sin\theta &= 2yz+2wx &\Rightarrow& \theta = \arcsin\left(2yz+2wx\right)\\
-\tan\phi &= \frac{2xz-2wy}{1-2x^2-2y^2} &\Rightarrow& \phi = \arctan\left(\frac{-2xz+2wy}{1-2x^2-2y^2}\right)
\end{align}$
Implementation Notes
There are two difficulties with implementing the conversion from quaternion to Tait-Bryan angles in a program. The first is what happens when the denominator in an $\arctan$ expression is 0. Most programming languages have solved this problem with a function arctan2(y,x) which computes $\arctan(y/x)$ while avoiding the possible divide-by-zero error. Check the documentation of your particular language, because some use the convention $\arctan(y/x)$ and some $\arctan(x/y)$.
The second difficulty involves the singularities of Tait-Bryan/Euler angles (gimbal lock), which are well-documented elsewhere. In this formulation, the singularities occur when $\theta = \pm \pi/2$. This can be solved with the technique described here under "Derivation of Equations": H-T-T-P://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToEuler/index.htm
*Anyone with the proper reputation, feel free to place that image inline and fix my extra links...
Best Answer
If you know the sequence of roll, pitch and yaw that you want to do in succession, you can find their matrices and multiply them together in the right order, and the three operations are performed "simultaneously."
The order is important because the group of rotations is nonabelian. If you can produce a matrix for each of these operations, then you can produce a single matrix combining all of them.