function R = ConstructRotation(axis, angle)
  % -----------------------------------------------------------------------
  % ConstructRotation - Construct the rotation matrix that rotates about
  %                     the given axis by the given angle
  %
  % Args: axis - 'x', 'y' or 'z': the axis of rotation
  %       angle - the angle of rotation
  %
  % Rtns: R - the 4 x 4 matrix of the rotation
  % -----------------------------------------------------------------------

  s = sin(angle); c = cos(angle);
  switch axis
  case 'x' % rotate about the x axis
    R = [
      1  0  0  0
      0  c -s  0
      0  s  c  0
      0  0  0  1
    ];
  case 'y' % rotate about the y axis
    R = [
       c  0  s  0
       0  1  0  0
      -s  0  c  0
       0  0  0  1
    ];
  case 'z' % rotate about the z axis
    R = [
       c -s  0  0
       s  c  0  0
       0  0  1  0
       0  0  0  1
    ];
  otherwise
    error('"%s" is not a valid axis', axis);
  end
end
