function M = rotmatrix(a,phi)
% MATLAB function computing the matrix representation of a rotation about the
% axis in direction of the vector \Blue{$\Va\in\bbR^3$} with angle \Blue{$\varphi$}.
% Note: \Red{Right hand rule} should be applied to fix the sense of the rotation, but,
% for the sake of simplicity, this is \textbf{not} done in this code. 
% 
% We know that rotations have a canonical matrix representation with respect to
% an orthonormal basis of \Blue{$\bbR^3$} with first vector in the direction of
% \Blue{$\Va$}. Thus the matrix representation of the rotation with respect to
% the standard Cartesian basis of unit vectors can be computed by using the
% formula (4.4.B) for the transformation of the matrix representation under a
% change of basis. Here the ``old basis'' is the special orthonormal basis, the
% ``new basis'' is the Cartesian basis. Hence the matrix \Blue{$\VS$} of
% formula (4.4.B) contains the coefficients of the Cartesian basis vectors with
% respect to the special orthonormal basis vectors. Since both bases are
% orthonormal, the transformation matrix will be \emph{orthogonal}. 

% First compute a special orthonormal basis with the normalized vector
% \Blue{$\Va$} as first element. This can be done by means of \Red{partial QR-decomposition}
% of the \Blue{$3\times 1$}-``matrix'' \Blue{$\Va$}. 
[Q,R] = qr(a); 
% Note that the \Blue{$\VQ$} obtained thus is the inverse of the basis
% transformation matrix \Blue{$\VS$} introduced above. However, this is not a
% big deal, because \Blue{$\VS = \VQ^{-1} = \VQ^T$} for the orthogonal matrix
% \Blue{$\VQ$}!

% This is the \textbf{matrix representation of the rotation} with respect to the
% specical orhtoonormal basis
c = cos(phi); s = sin(phi);
D = [1 , 0  , 0  ;...
     0 , c , -s  ;...
     0 , s , c  ];
% Finally compute the transformation of the matrix representation according to
% Formula (4.4.B). Note that a prime in MATLAB denotes transposition.
M = Q*D*Q';
