Given unit vectors $a$ and $b$, I want to find a rotation $R$ such that:
$$ b = Ra $$
and
$$ R = R_y(\phi)R_x(\theta) $$
where $R_x(\theta)$ is a right-handed rotation around the X axis by an angle of $\theta$ and $R_y(\phi)$ is the same for Y, for some $\phi$ and $\theta$.
I don't particularly want the actual angles, I just want to make sure that $R$ is a product of these two rotations. I approached this problem by looking at the matrix representation of $R$ and solving for the angles, but it seems inelegant. Is there a more direct way, perhaps using quaternions?
Here's the approach I have now:
If I expand R, I get:
$$ b = \begin{bmatrix} cos(\phi) & sin(\phi) sin(\theta) & sin(\phi) cos(\theta)\\ 0 & cos(\theta) & -sin(\theta)\\ -sin(\phi) & sin(\theta) cos(\phi) & cos(\phi) cos(\theta)\\ \end{bmatrix} a $$
Using the harmonic addition theorem,
$$ x \cdot cos(\phi)+y \cdot sin(\phi) = c \cdot cos(\phi + \lambda) = z\\ c = sgn(x) \sqrt{x^2+y^2}\\ \lambda = \operatorname{arctan} \left(- \frac{y}{x} \right)\\ \phi = \pm \operatorname{arccos} \left(\frac{z}{c} \right) - \lambda $$
I can solve the second row for $\theta$. This has two solutions, and I want the one where $| \theta | < \pi / 2$. Knowing this, I can again use the harmonic addition theorem to solve the first row for $\phi$. This again produces two solutions, and I pick the one that produces the correct $b_3$. This works, but the biggest problem is that solving the harmonic additions is undefined in several places ($a_2 = 0$, $b_2 > c$, $b_1 > c$). Is there a better way?