What's the most elegant way of rotating a 3-dimensional co-ordinate system? For two dimensional rotation of $x$ and $y$ axes anticlockwise by $\varphi$, the equation that transforms $P(x,y) \rightarrow P(x',y')$,  $x'=x \cos(\varphi)+y \sin(\varphi)$ and $y'=y \cos(\varphi)- x \sin(\varphi)$  are nice enough, but when I tried to do this for $\varphi$ and $\theta$ in three dimensions (that is, another iteration of the same thing) the result was horribly ugly.
So, starting and ending with a Cartesian point, what's the most elegant way of stating the ending co-ordinates as a function of the starting co-ordinates and other things, like the angle of rotation about various axes? 
 A: If you know the rotation axis and angle, I would consider Rodrigues' rotation formula.  
$$v' = v \cos \theta + (k \times v) \sin \theta + (1-\cos \theta) (k \cdot v) k$$
where $k$ is the unit vector along the axis of rotation and $\theta$ is the angle of rotation.  $v$ is the original vector, and $v'$ is the rotated vector.  You can then work out the components of $v'$ in terms of components of $v$ and $k$.
If you don't know the overall axis and angle of rotation (for example, because you're performing multiple rotations along several different axes), you may be forced to multiply rotation matrices and/or quaternions.  Quaternions are particularly easy to convert back and forth with axis-angle representation (and they have the benefit of being able to chain rotations through multiplication, like with rotation matrices), but they're somewhat more involved to work with.  Very useful, though, if you have the time and inclination to work with them.
A: Use a 3 by 3 matrix. 
Start by identifying the axis of rotation in spherical coordinates. It helps to make it unit length: $$[\vec{u}_3]_{\text{std}}=\begin{bmatrix}\cos\theta\cos\phi\\\sin\theta\cos\phi\\\sin\phi\end{bmatrix}$$ Extend to an orthonormal basis for $\Bbb{R}^3$. One option is $$\mathcal{B}=\left\{\begin{bmatrix}\sin\theta\\-\cos\theta\\0\end{bmatrix},\begin{bmatrix}\cos\theta\sin\phi\\\sin\theta\sin\phi\\-\cos\phi\end{bmatrix},\begin{bmatrix}\cos\theta\cos\phi\\\sin\theta\cos\phi\\\sin\phi\end{bmatrix}\right\}$$
So the matrix $$P=\begin{bmatrix}\sin\theta&\cos\theta\sin\phi&\cos\theta\cos\phi\\-\cos\theta&\sin\theta\sin\phi&\sin\theta\cos\phi\\0&-\cos\phi&\sin\phi\end{bmatrix}$$
is a change of basis matrix from $\mathcal{B}$ to the standard basis. Note that $P^{-1}=P^T$, since the columns are orthonormal.
Now if you want to rotate around this axis by angle $\psi$ (oriented by the right-hand rule) we can first construct $$R=\begin{bmatrix}\cos\psi&-\sin\psi&0\\\sin\psi&\cos\psi&0\\0&0&1\end{bmatrix}$$ which rotates by angle $\psi$ around the third vector in any given basis. The matrix that accomplishes your rotation is $M=PRP^T$. That is, for any $\vec{v}$ expressed as a column vector using standard coordinates, and letting $\vec{w}$ represent the rotated vector: $$\begin{align}
M[\vec{v}]_{\text{std}}&=PRP^T[\vec{v}]_{\text{std}}\\
&=PR[\vec{v}]_{\mathcal{B}}\\
&=P[\vec{w}]_{\mathcal{B}}\\
&=[\vec{w}]_{\text{std}}\\
\end{align}$$
I wouldn't actually multiply out $PRP^T$; doing so would only obscure the linear algebra theory that went in to finding $M$. (Of course if you are in the position of having to rotate many vectors by this same angle around this same axis, then for computational efficiency sake it makes sense to multiply them out.)
