I looking into the problem of finding solutions (there could be more than one) for the rotation matrix $R = [u_1, u_2, u_3] $ where $u_1, u_2, u_3$ are unit vectors and mutually orthogonal, and that in addition saisfies the following
$$ u_1 \cdot i = \alpha \\u_2 \cdot j = \beta \\u_3 \cdot k = \gamma $$
where $i,j,k$ are the unit vectors along the $x,y,z$ axes respectively.
In other words, I would like to find all the off-diagonal entries of orthogonal matrices $R$ whose diagonal entries are known.
I am looking for either a closed-form solution (this would be preferrable), or an algorithm for finding the solutions numerically.
The problem is inspired by this recent MSE post
The rotation matrix for a unit quaternion $\mathbf q = q_r + q_i\mathbf i + q_j\mathbf j + q_k\mathbf k$ is
$$\begin{bmatrix} 1 - 2(q_j^2 + q_k^2) & 2(q_iq_j - q_kq_r) & 2(q_iq_k + q_jq_r) \\ 2(q_iq_j + q_kq_r) & 1 - 2(q_i^2 + q_k^2) & 2(q_jq_k - q_iq_r) \\ 2(q_iq_k - q_jq_r) & 2(q_jq_k + q_iq_r) & 1 - 2(q_i^2 + q_j^2) \end{bmatrix}.$$
This gives us four linear equations in $q_r^2, q_i^2, q_j^2, q_k^2$:
\begin{gather*} q_r^2 + q_i^2 + q_j^2 + q_k^2 = 1, \\ 1 - 2(q_j^2 + q_k^2) = α, \quad 1 - 2(q_i^2 + q_k^2) = β, \quad 1 - 2(q_i^2 + q_j^2) = γ, \end{gather*}
which are easily solved:
\begin{gather*} q_r = ±\frac{\sqrt{1 + α + β + γ}}{2}, \quad q_i = ±\frac{\sqrt{1 + α - β - γ}}{2}, \\ q_j = ±\frac{\sqrt{1 - α + β - γ}}{2}, \quad q_k = ±\frac{\sqrt{1 - α - β + γ}}{2}. \end{gather*}
This yields up to 8 different rotation matrices (since $-\mathbf q$ yields the same matrix as $\mathbf q$).
There are also up to 8 orthogonal matrices that are rotary reflections rather than rotations, obtained as the negations of the rotation matrices for $-α, -β, -γ$.