function xyz_rotated = rotation_xyz ( xyz, axis, theta )

%*****************************************************************************80
%
%% rotation_xyz() rotates (x,y,z) through theta radians around an axis.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    12 December 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John D Cook,
%    Advantages of Redundant Coordinates,
%    11 December 2019,
%    https://www.johndcook.com/blog/
%
%  Input:
%
%    real xyz(3,1): the point to be rotated.
%
%    real axis(3,1): the axis vector.
%
%    real theta: the rotation angle in radians.
%
%  Output:
%
%    real xyz_rotated(3,1): the rotated point.
%
  p = [ 0; xyz ];
  th = theta / 2.0;
  q = cos ( th ) * [ 1.0; axis ];
  q_inverse = q8_inverse ( q );
  qp = q8_multiply ( q, p );
  qpq = q8_multiply ( qp, q_inverse );
  xyz_rotated = qpq(2:4);

  return
end

