Given a unit vector $\vec{n}$, find the matrix of rotation about $\vec{n}$. In $\mathbb{R}^3$ suppose I have an arbitrary vector $\vec{n}$.  I want to find the rotation matrix about $\vec{n}$ through an angle $\theta$.
I can develop the rotation matrix in two dimensions and that makes it clear to me that the rotation matrices about the basis vectors are 
$$R_x=\begin{pmatrix}1 & 0 & 0 \\0 & \cos\theta &-\sin\theta\\ 0 & \sin\theta & \cos\theta\end{pmatrix}$$
and similarly for $R_y$ and $R_z$.  I've also looked at this answer: Vector rotation but I'm trying to follow the steps to deriving it, not just looking for the answer.  Also, I know nothing about tensors.
I looked at this answer: Getting a transformation matrix from a normal vector  I believe I understand why we would want a rotation that would take the plane normal to the normal, and align it with the $x$-$y$ plane, since having done that we might be able to then follow with a rotation about the $z$ axis to get the resulting vectors aligned with $x$ and $y$ axes.  The composition of these rotations would map the new bases to the old, and the inverse would map the old to the new.
So the task is to now rotate $\vec{n}$ onto $\vec{k}$, and I can see (somewhat) why we would want a vector $\vec{u}$ in the $x$-$y$ plane which is normal to this rotation, and how to get it by solving $\vec{n}\cdot \vec{u} = 0$ with $\vec{u} = a\vec{i}+b\vec{j}$.  
But once having $\vec{u}$, I'm not sure how I would find the matrix of the rotation about $\vec{u}$ which takes $\vec{n}$ onto $\vec{k}$, and the answer given there points to a Wikipedia article to do this.  But again I don't just want the answer, I want to understand the derivation.  
Besides that I'm somewhat confused by the methodology in this sense:  Don't we have a situation that is not very much better than the situation we started with?  Right now I have a vector $\vec{u}$ and want to rotate $\vec{n}$ about it ... Well isn't that the same as what we're trying to do generally?  Is anything gained by the fact that $\vec{u}$ is in the $x$-$y$ plane which makes that easier than just doing the same task with $\vec{n}$ directly?
And I would really prefer a derivation that went more like the way that the 2-D rotation proceeded, since I understood that clearly.  It would be nice to consider vector $\vec{n}$ and the rotation of $\vec{i}$ about it, write those coordinates into the first column of the answer, and then do likewise for $\vec{j}$ and $\vec{k}$.  Of course, perhaps we don't do it that way because it's very hard, and hence why I can't see how to do it that way on my own.  
 A: I am going to derive the rotation around a unit vector very slowly.
Given a unit vector $N$ and a vector $X$, the task is to rotate $X$ around $N$.
First suppose we are in the simple case in which the vector $X$ lies in the plane perpendicular to $N$, so:
$X \cdot N = 0$
To rotate $X$ on the plane we first need to find a orthonormal basis on that plane, made of two vectors $U_1$ and $U_2$.
$U1 \cdot U2 = 0$
$U1 \cdot U1 = U2 \cdot U2 = 1$
Since $X$ is already on the rotation plane we can use it as one of our basis vectors.
$U_1 = X / \|X\|$
$U_2 = N \times U_1$
Then we can form a $3 \times 3$ change-of-basis matrix $M$ that has as columns the vectors $U_1$, $U_2$ and $N$.
$M = [U_1 U_2 N]$
The matrix $M$ has the property of sending the cartesian basis vectors $E_1 = [1,0,0]$, $E_2 = [0,1,0]$ and $E_3 = [0,0,1]$ to our rotation plane.
You can check that:
$M E_1 = U_1$
$M E_2 = U_2$
$M E_3 = N$
Since $M$ is an orthonormal matrix its transpose is equal to the inverse $M^T = M^{-1}$ so:
$M^T U_1 = E_1$
$M^T U_2 = E_2$
$M^T N = E_3$
Now, let $R_3$ be a $3\times3$ rotation matrix around the cartesian vector $E_3$ (i.e., a rotation in the plane $E_1 E_2$).
$R = M R_3 M^T$
$R$ is a rotation in the plane $U_1$ and $U_2$ around $N$. You can check that:
$R X = X'$
$(M R_3 M^T) X = X'$
Since $M^T X = \|X\| E_1$ we have:
$M R_3 (M^T X) = M R_3 (\|X\| E_1)$
Lets call $E_1'$ to the rotation of $E_1$ by matrix $R_3$.
$\|X\| M (R_3 E_1) = \|X\| M E_1' $
$M (\|X\| E_1')  = X'$
And we get back:
$R X = X'$
Now lets remove the assumption that $X$ is in the rotation plane. We still need to find two vectors $U_1$ and $U_2$ in the rotation plane, so we project $X$ to the rotation plane to find $U_1$ and then we get $U_2$ as the cross product of $N$ and $U_1$.
$U1 = (X - (X \cdot N) N) / \|X - (X \cdot N) N\|$
$U_2 = N \times U_1$
$X$ is a linear combination of vectors $U_1$ and $N$. So applying the matrix $M^T = [U_1 U_2 N]^T$ to $X$ gives:
$M^T X = a E_1 + b E_3$
Where:
$a = \|X - (X \cdot N) N\| $
$b = (X \cdot N)$
Since $E_3$ is not changed by the rotation $R_3$. The rotation $R = M R_3 M^T$ around $N$ is going to change only the component of $X$ projected on the rotation plane. That does the "trick" of rotation around an axis.
Up to now we have used the vector $X$ for creating the vector $U_1$ in the rotation plane. But what happen if we only have the vector $N$ and the rotation angle? In that case we can choose any cartesian vector $E1$, $E2$ or $E3$ provided that they are linearly independent with $N$.
$U_1 = N \times E_1$   or   $U_1 = N \times E_2$  or  $U_1 = N \times E_3$
Of course there are better ways to choose the basis vectors, see the Rodrigues Formula derivation in wikipedia and perhaps you can understand now how it works.
https://en.m.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Derivation
