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;
}

 A: 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$.
A: The rotation you are looking for is in fact the combination of


*

*a rotation of $\pi/4$ around the $z$ axis, followed by a rotation of $-(\pi/2-\alpha)$ (right hand rule) around the $y'$ axis (note the $'$), or

*a rotation of $-(\pi/2-\alpha)$ around the $y$ axis (no $'$), followed by a rotation of $\pi/4$ around the $z$ axis.


Denoting the rotation matrices as (the column vectors are the coordinates of the rotated ref. relative to the base ref.)
$$
{\bf R}_{\,{\bf z}} (\varphi ) = \left( {\matrix{
   {\cos \varphi } & { - \sin \varphi } & 0  \cr 
   {\sin \varphi } & {\cos \varphi } & 0  \cr 
   0 & 0 & 1  \cr 
 } } \right)\quad {\bf R}_{\,{\bf y}} (\beta ) = \left( {\matrix{
   {\cos \beta } & 0 & {\sin \beta }  \cr 
   0 & 1 & 0  \cr 
   { - \sin \beta } & 0 & {\cos \beta }  \cr 
 } } \right)
$$
then the required rotation is
$$
\eqalign{
  & {\bf R} = {\bf R}_{\,{\bf y'}} ( - \pi /2 + \alpha )\;{\bf R}_{\,{\bf z}} (\pi /4)\; = {\bf R}_{\,{\bf z}} (\pi /4)\;{\bf R}_{\,{\bf y}} ( - \pi /2 + \alpha )\; =   \cr 
  &  = \left( {\matrix{
   {\sqrt 3 /3} & { - \sqrt 2 /2} & { - \sqrt 3 \sqrt 2 /6}  \cr 
   {\sqrt 3 /3} & {\sqrt 2 /2} & { - \sqrt 3 \sqrt 2 /6}  \cr 
   {\sqrt 3 /3} & 0 & {\sqrt 3 \sqrt 2 /3}  \cr 
 } } \right) \cr} 
$$
because the rotation around the new $y$ axis, resulting upon the rotation
around the $z$ axis, is
$$
{\bf R}_{\,{\bf y'}} ( - \pi /2 + \alpha )\; = {\bf R}_{\,{\bf z}} (\pi /4)\;{\bf R}_{\,{\bf y}} ( - \pi /2 + \alpha )\;{\bf R}_{\,{\bf z}} (\pi /4)^{\, - \,{\bf 1}} 
$$
To conclude that your $xy$ vector, parallel to $y'$, should be $(-1,1,0)$ and the rotation angle (if RH rule is respected) shall be $-\pi/2+\alpha$.
