How to rotate a vector to get the magic angle? I have vector x = (1,0,0). I need to find a set of rotations which transform it in such a way that it forms the magic angle with each axis. The magic angle is equal to $\arccos \left(\frac{1}{\sqrt{3}}\right)$.
It seems intuitive to rotate 45 degrees over z axis and then 90-magic over xy axis, but it doesn't work. Or perhaps there is something wrong with my code? I past it below.
Eigen::Vector3d x(1, 0, 0);
Eigen::Vector3d y(0, 1, 0);
Eigen::Vector3d z(0, 0, 1);
Eigen::Vector3d xy(1, 1, 0);

Eigen::Matrix3d rot = turboRotation(90-magic, xy);
rot = rot * turboRotation(45, z);

cout << acos((rot*x).dot(x))*180/M_PI << endl; //checking the angles
cout << acos((rot*x).dot(y))*180/M_PI << endl;
cout << acos((rot*x).dot(z))*180/M_PI << endl;

turboRotation function is:
Eigen::Matrix3d turboRotation(double deg, Eigen::Vector3d rotation_line)
{
    rotation_line /= rotation_line.norm();
    double a = deg*(M_PI / 180);
    Eigen::Matrix3d rotationMatrix;
    double x = rotation_line(0);
    double y = rotation_line(1);
    double z = rotation_line(2);
    rotationMatrix << cos(a)+x*x*(1-cos(a)), x*y*(1 - cos(a))-z*sin(a), x*z*(1 - cos(a))+y*sin(a),
                      y*x*(1 - cos(a))+z*sin(a), cos(a)+y*y*(1 - cos(a)), y*z*(1 - cos(a)) + x*sin(a),
                      z*x*(1 - cos(a)) + y*sin(a), z*y*(1 - cos(a)) + x*sin(a), cos(a) + z*z*(1 - cos(a));
    return rotationMatrix;
}


Let $a=\tfrac{1}{\sqrt{3}}$, $b=\tfrac{1}{\sqrt{2}}$ and $c=ab.$
Consider matrix:
$$R=\left(\begin{array}{rrr}a & b & c\\a & -b & c\\a & 0 & -2c \end{array}\right).$$
It is an orthogonal matrix (its columns $C_1,C_2,C_3$ are orthogonal with norm 1). As $\det R = 6abc = 1$, it is a rotation matrix. 
Moreover, $R$ sends vector $U$ with coordinates $1,0,0$ onto $C_1$ with a dot product equal to : $\cos(\theta)=1 \times a+0 \times a + 0 \times a = a = \tfrac{1}{\sqrt{3}}$, thus with the correct "magic" :) angle.
Moreover, setting $d=\sqrt{2}a,$ we have the following decomposition:
$$R=\left(\begin{array}{rr|r}b & -b & 0\\b & b & 0\\ \hline 0 & 0 & 1 \end{array}\right)\left(\begin{array}{r|rr}1 & 0 & 0\\ \hline 0 & -d & -a\\0 & a & -d \end{array}\right)\left(\begin{array}{rr|r}0 & -1 & 0\\1 & 0 & 0\\ \hline 0 & 0 & 1 \end{array}\right)$$
with respect to $x$, $z$ and $x$ axes in succession ($y$ axis is not used, as usual with Euler angles, here $\pi/2$ for the first one, $arcsin(\Phi)$ for the 2nd angle ($\Phi$ denotes the magic angle), and $\pi/4$ for the last rotation to be applied.
Remark: although I had no choice for $C_1$, I have taken $C_2$ as an arbitrary vector orthogonal to $C_1$, subsequently normalized ; $C_3$ has been taken as the cross-product of $C_1$ and $C_2$.
