matrix for a "canonical" rotation Let $p$ be an arbitrary point on the unit sphere $$ S = \{(x,y,z) \mid x^2+y^2+z^2 = 1 \},$$
other than the north or south poles $(0,0,\pm 1)$.
There is a one-dimensional family of rotations which take $n = (0,0,1)$ to $p$, but one rotation is canonical: the one that keeps $n$, $p$, and $0$ in a plane.  In other words the cross product $u = n \times p$ should be an axis of rotation.
Is there a nice way to write down the matrix $M \in SO(3)$ for this canonical rotation?
I observe that we have three equations:  (1) $M^T M = I$, (2) $M n = p$, and (3) $M u = u $.
Since I want $M \in SO(3)$ rather than simply $M \in O(3)$ we also have $\det M = 1$.
This should be enough to determine $M$, but I still don't see an obvious way to write down a formula.
 A: Let the point $p$ be given in cartesian coordinates by $(\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\phi)$, so that $\theta$ and $\phi$ are the usual spherical coordinates on $S^2$.  Then the matrix that rotates $n=(1,0,0)$ to $p$, and also preserves the plane passing through both of them and the origin, is given by  
$M=\left(\begin{array}{ccc}\sin^2\phi + \cos\theta\cos^2\phi & (-1+\cos\theta)\sin\phi\cos\phi &  \sin\theta\cos\phi\\
(-1+\cos\theta)\sin\phi\cos\phi & \cos^2\phi+\cos\theta\sin^2\phi & \sin\theta\sin\phi\\
-\sin\theta\cos\phi & -\sin\theta\sin\phi & \cos\theta
\end{array}\right).$
It's straightforward but tedious to check that this is indeed an orthogonal matrix with determinant $1$ that sends $n$ to $p$ and fixes $u=n\times p=(-\sin\theta\sin\phi,\sin\theta\cos\phi,0)$, and hence the desired transformation.
I produced this matrix by considering that, for each fixed $\phi$, what you desire is a one-parameter subgroup of $SO(3)$: the group of rotations that fix $u$.  It therefore is the family of exponentials of an infinitesimal transformation, i.e. an element of the Lie algebra $\frak{so}\mathrm{(3)}$.  For example, a tiny rotation tilting the $z$-axis toward the $x$-axis differs from the identity by 
$\left(\begin{array}{ccc}
 &  & \epsilon\\
 &  & \\
-\epsilon & & 
\end{array}\right)$, and
$\exp\left(\begin{array}{ccc}
 &  & \theta\\
 &  & \\
-\theta & & 
\end{array}\right) = \left(\begin{array}{ccc}
\cos\theta &  & \sin\theta\\
 & 1 & \\
-\sin\theta & & \cos\theta
\end{array}\right).$
Similarly, the Lie algebra element corresponding to an infinitesimal rotation of the $z$-axis in the direction of the $y$-axis is given by $\left(\begin{array}{ccc}
 &  & \\
 &  & \epsilon\\
 & -\epsilon& 
\end{array}\right)$, and
$\exp\left(\begin{array}{ccc}
 &  & \\
 &  & \theta\\
 & -\theta & 
\end{array}\right) = \left(\begin{array}{ccc}
1 &  & \\
 & \cos\theta & \sin\theta\\
 & -\sin\theta & \cos\theta
\end{array}\right).$
We want the Lie algebra element which has components $\cos\phi$ toward the $x$-axis and $\sin\phi$ toward the $y$-axis; thus we desire the exponential
$M=\exp\left(\begin{array}{ccc}
 &  & \theta\cos\phi\\
 &  & \theta\sin\phi\\
-\theta\cos\phi & -\theta\sin\phi & 
\end{array}\right)=\exp(A).$
This matrix $A$ is diagonalizable, so its exponential $M$ can be computed by hand, which is how I arrived at my answer above.
