3D rotation of a square Suppose we have a square in 3D space with its normal pointing in the direction of the positive $z$ axis, i.e. $(0,0,1)$. The four vertices of the square are $(1,1,0)$, $(-1,1,0)$, $(-1,-1,0)$, $(1,-1,0)$. 
Say I now rotate this normal using an azimuth and zenith angle $$0 <\phi <2\pi,\quad 0 <\theta <\pi/2$$
 so that $n'=(\cos \phi \sin\theta, \sin\phi\sin\theta,\cos\theta)$.
I want to obtain the new positions of the four points. Do we just apply rotations about the $z$ and $x$ axis using azimuth and zenith angles to the four points to get their new positions? In what order?
 A: $\newcommand{\Reals}{\mathbf{R}}\newcommand{\Basis}{\mathbf{e}}\newcommand{\vec}{\mathbf{x}}$Because three-dimensional rotations about the origin don't commute (unlike rotations of the plane about the origin), there's not enough information in your question to give a unique answer.
My crystal ball says you probably have in mind a rotation about a single axis sending the old normal vector $n = (0, 0, 1)$ to the new normal $n' = (\cos\phi \sin\theta, \sin\phi \sin\theta, \cos\theta)$. The axis of your rotation is therefore
$$
\Basis_{3} = \frac{n \times n'}{\|n \times n'\|} = (-\sin\phi, \cos\phi, 0).
$$
Extend this to an orthonormal basis of $\Reals^{3}$ by taking
$$
\Basis_{1} = (0, 0, 1),\qquad
\Basis_{2} = \Basis_{3} \times \Basis_{1} = (\cos\phi, \sin\phi, 0).
$$
Since you're rotating about $\Basis_{3}$ by $\theta$, with $\Basis_{1}$ moving toward $\Basis_{2}$, your rotation sends the ordered triple $(\Basis_{i})_{i=1}^{3}$ to the ordered triple $(\Basis_{i}')_{i=1}^{3}$ with
$$
\Basis_{1}' = (\cos\theta) \Basis_{1} + (\sin\theta) \Basis_{2},\qquad
\Basis_{2}' = -(\sin\theta) \Basis_{1} + (\cos\theta) \Basis_{2},\qquad
\Basis_{3}' = \Basis_{3}.
$$
Express each vertex as a linear combination of the $(\Basis_{i})_{i=1}^{3}$ and use the resulting coefficients with the $(\Basis_{i}')_{i=1}^{3}$ to get the transformed vertices. The end result is that a vertex $\vec = (x, y, z)$ maps to
$$
(\vec \cdot \Basis_{1}) \Basis_{1}' + (\vec \cdot \Basis_{2}) \Basis_{2}' + (\vec \cdot \Basis_{3}) \Basis_{3}'.
$$
The result with $\theta = \pi/4$ and $\phi = \pi/3$ looks like this:

