Given are two unnormalized, non-parallel vectors, the rotation axis n and the vector r to be rotated.
Decompose r into two orthogonal components:

Clearly, n and x are orthogonal. Define further y as a cross product, a vector orthogonal to the plane containing n, r, and x,

As is well-known the cross product can be written as a matrix-vector product

The matrix N has as general element

where εαβγ is the antisymmetric Levi-Civita tensor.
For further use we compute normalization constants of x and y,


and divide the two

When we rotate r over an angle φ around n, the component of r along n is unchanged, while the component x of r perpendicular to n becomes x′

Hence the rotated vector r′ is

We may introduce the dyadic product of the vector n with itself, which has the form of a 3 × 3 symmetric matrix, and write

Now,
![{\displaystyle \mathbf {r} '=\left[\cos \phi \;\mathbf {E} +{\frac {(1-\cos \phi )}{n^{2}}}\;{\big (}\mathbf {n} \otimes \mathbf {n} {\big )}+{\frac {1}{n}}\sin \phi \;\mathbf {N} \right]\mathbf {r} ,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f4327fbae77df335801e79619cf265e58a59c35a)
where E is the identity matrix. The quantity between square brackets is the matrix R that rotates r around n over an angle φ. This equation is very well-known and was first derived by Leonhard Euler [check].
A general element of R is

where the unit vector is

Consider now two non-parallel unit vectors f and t making an angle φ

We want to find the matrix that rotates f (the "from" vector) to t (the "to" vector). An obvious way is using the cross product f×t as a rotation axis and rotating f over φ. We can use the equation just derived by substituting

Hence

so that an element of the rotation matrix takes the form

which is Eq. (1) of Möller and Hughes (1999).