Resolution of a rotation Can the rotation of a point about some arbitrary axis $\vec{n}$ be described by some rotation of the point about z-axis followed by some rotation about y axis?
 A: Say we want to rotate point $(x_0, y_0, z_0)$ to $(x_1, y_1, z_1)$ where $\rho= \sqrt{x_0^2+ y_0^2+ z_0^2}= \sqrt{x_1^2+ y_1^2+ z_1^2}$.  First rotate about the z-axis so that $(x_0, y_0, z_0)$ is rotated to $(r, 0, z_0)$ where $r= \sqrt{x_0^2+ y_0^2}$.  That can be done with a matrix of the form $\begin{pmatrix}a & b & 0 \\ -b & a & 0 \\ 0 & 0 & 1\end{pmatrix}$.  We want $\begin{pmatrix}a & b & 0 \\ -b & a & 0 \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}x_0 \\y_0 \\z_0\end{pmatrix}= \begin{pmatrix}ax_0+ by_0 \\ -bx_0+ ay_0 \\ z_0\end{pmatrix}= \begin{pmatrix} r\\ 0 \\ z_0\end{pmatrix}$.  So we have $ax_0+ by_0= r$ and $-bx_0+ ay_0= 0$ and we want to solve for a and b.  Multiply the first equation by $x_0$ to get $ax_0^2+ bx_0y_0= rx_0$.  Multiply the second equation by $y_0$ to get $-bx_0y_0+ ay_0^2= 0$.  Adding those eliminates b: $a(x_0^2+ y_0^2)= ar^2= rx_0$.  $a= \frac{x_0}{r}$.  Then $-bx_0+ ay_0= 0$ becomes $-bx_0+ \frac{x_0y_0}{r}= 0$ so $bx_0= \frac{x_0y_0}{r}$ and $b= \frac{y_0}{r}$.
The first matrix, rotating about t the z-axis, is $\begin{pmatrix}\frac{x_0}{r} & \frac{y_0}{r} & 0 \\ -\frac{y_0}{r} & \frac{x_0}{r} & 0 \\ 0 & 0 & 1\end{pmatrix}$.
Now we want to rotate around the y-axis so that $\begin{pmatrix} r \\ 0 \\ z_0\end{pmatrix}$ rotates to $\begin{pmatrix} 0 \\ 0 \\ \rho\end{pmatrix}$. Rotating about the y-axis requires a matrix of the form $\begin{pmatrix} a & 0 & b \\ 0 & 1 & 0 \\ -b & 0 & a\end{pmatrix}$.  Then we must have $\begin{pmatrix} a & 0 & b \\ 0 & 1 & 0 \\ -b & 0 & a\end{pmatrix}\begin{pmatrix} r \\ 0 \\ z_0\end{pmatrix}= \begin{pmatrix}ar+ bz_0 \\ 0 \\ -br+ az_0 \end{pmatrix}= \begin{pmatrix} 0 \\ 0 \\ \rho \end{pmatrix}$ where $\rho= \sqrt{x_0^2+ y_0^2+ z_0^2}= \sqrt{r^2+ z_0^2}$.  Now we have the two equations $ar+ bz_0= 0$ and $-br+ az_0= \rho$ that we want to solve for a and b.
Multiply the first equation by r to get $ar^2+ brz_0= 0$ and multiply the second equation by $z_0$ to get $-brz_0+ az_0^2= \rho z_0$.  Adding the two equations eliminate b: $ar^2+ az_0^2= a(r^2+ z_0^2)= a\rho^2= \rho z_0$ so $a= \frac{z_0}{\rho}$.  Now, $ar+ bz_0= 0$ becomes $\frac{r z_0}{\rho}+ bz_0= 0$ so $bz_0= -\frac{r z_0}{\rho}$ and $b= -\frac{r}{\rho}$.
The second matrix, rotating around the y-axis, is $\begin{pmatrix}\frac{z_0}{\rho} & 0 & -\frac{r}{\rho} \\ 0 & 1 & 0 \\ \frac{r}{\rho} & 0 & \frac{z_0}{\rho} \end{pmatrix}$.
A: First you have to bring one of the axis (for instance $z$) onto $\bf n$.   
To do this, you can rotate around $z$ to bring e.g. $y$ normal to the plane $z, \bf n$, then rotate around the new $y$ to bring $z$ to $\bf n$.
After that, then you can rotate around the (new) $z$ for the required angle.
If the point to be rotated is given in advance, in the original base system, then before applying the last required rotation around $\bf n$, you shall express $P$ in the new coordinates.
Finally, the new point $P'$ will lay on a centered sphere same as $P$, so it will always possible to pass from one to the other with a latitude + longitude steps.
You could use spherical trigonometry, but I do not think that the process would be simpler.
