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?
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.