Please note: What I'm looking for is $\alpha_0$ and $\beta_0$ in a form similar to Eq.(3) that will eventually give Eq.(1). Note that the inertia tensor is not a rotation matrix so I'm not trying to convert a rotation matrix into Euler form since $\mathbb{I}$ is not a rotation matrix to start with. As a result the proposed duplicates are not really relevant to what I'm trying to do.
It is well-known that the inertia tensor can be brought to diagonal form by a rotation: $$ \mathbb{I}=\left(\begin{array}{ccc} I_{xx}&I_{xy}&I_{xz}\\ I_{xy}&I_{yy}&I_{yz}\\ I_{xz}&I_{yz}&I_{zz}\end{array}\right) \to \mathbb{I}'=R(\Omega) \mathbb{I}R^{-1}(\Omega) =\left(\begin{array}{ccc} I_1&0&0\\ 0&I_2&0\\ 0&0&I_3\end{array}\right)\, , \qquad R(\Omega)\in SO(3)\, . $$ The rotation $R(\Omega)$ can be factorized in terms of Euler angles - say $$ R(\Omega)=R_z(\gamma)R_y(\beta) R_z(\alpha)\, . $$ Is there a systematic way of finding the Euler angles $(\gamma,\beta,\alpha)$ that will diagonalize the inertia tensor in terms of the entries of the $\mathbb{I}$ tensor?
Since $$ R_z(\alpha)=\left( \begin{array}{ccc} \cos (\alpha ) & \sin (\alpha ) & 0 \\ -\sin (\alpha ) & \cos (\alpha ) & 0 \\ 0 & 0 & 1 \\ \end{array} \right)\qquad R_y(\beta)=\left( \begin{array}{ccc} \cos (\beta ) & 0 & \sin (\beta ) \\ 0 & 1 & 0 \\ -\sin (\beta ) & 0 & \cos (\beta ) \\ \end{array} \right)\, . $$ one strategy is to find $\beta,\alpha$ so that $$ R_y(\beta)R_z(\alpha) \mathbb{I} R_z(-\alpha)R_y(-\beta)= \left(\begin{array}{ccc} A&B&0\\ B&C&0\\ 0&0&E \end{array}\right)\, . \tag{1} $$ A final rotation by $R_z(\gamma)$ would diagonalize the remaining $2\times 2$ subspace.
I started with the idea of finding consecutive transformations that would do the trick, i.e. started by finding $\alpha_0$ so that $$ R_z(\alpha_0) \mathbb{I} R_z(-\alpha_0) =\left(\begin{array}{ccc} a&b&c\\ b&d&0\\ c&0&e \end{array}\right)\, . \tag{2} $$ Finding $\alpha_0$ follows from matrix multiplication by setting $$ \cos(\alpha_0)=\frac{\mathbb{I}_{xz}}{\sqrt{\mathbb{I}_{xz}^2+\mathbb{I}_{yz}^2}}\, ,\quad \sin(\alpha_0)=\frac{\mathbb{I}_{yz}}{\sqrt{\mathbb{I}_{xz}^2+\mathbb{I}_{yz}^2}}\, , \tag{3} $$ so $\alpha_0$ is given explicitly in terms of entries of $\mathbb{I}$. The entries $a, b, c, d,e$ that result from this matrix multiplication now depend on $\alpha_0$ and other entries of $\mathbb{I}$, but the inertia tensor now has the form of Eq.(2) so I use this as the starting point for the next step.
I then tried to find $R_y(\beta_0)$ to make $c=0$ but this simply moves the $0$ to another location: $$ R_y(\beta_0)R_z(\alpha_0) \mathbb{I} R_z(-\alpha_0)R_y(\beta_0) = \left(\begin{array}{ccc} a'&b'&0\\ b'&d'&c'\\ 0&c'&e \end{array}\right)\, . $$ One can try variations on the $R_zR_yR_z$ form, say $R_xR_yR_x$ but it appears the same basic result remains.
Thus it seems the only way to proceed is to start with $R_y(\beta)R_z(\alpha) \mathbb{I} R_z(-\alpha)R_y(\beta)$ and simulteaneously find $\alpha,\beta$ that will bring $\mathbb{I}$ to the block diagonal form of (1), and this leads to highly non-trivial equations to solve, with possibly no analytical solution in terms of $I_{ab}$.
If it is not possible to solve the general case, is it possible to solve some specific specialized cases?
