function w = rotation_axis_vector ( axis, angle, v )

%*****************************************************************************80
%
%% rotation_axis_vector() rotates a vector around an axis vector in 3D.
%
%  Discussion:
%
%    Thanks to Cody Farnell for correcting some mistakes in an earlier
%    version of this routine.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 July 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real AXIS(3), the axis vector for the rotation.
%
%    real ANGLE, the angle, in radians, of the rotation.
%
%    real V(3), the vector to be rotated.
%
%  Output:
%
%    real W(3), the rotated vector.
%
  axis = axis(:);
  v = v(:);
%
%  Compute the length of the rotation axis.
%
  u = axis;
  axis_norm = norm ( u );

  if ( axis_norm == 0.0 )
    w = zeros ( 3, 1 );
    return
  end

  u = u / axis_norm;
%
%  Compute the dot product of the vector and the rotation axis.
%
  udotv = u' * v;
%
%  Compute the parallel component of the vector.
%
  parallel = udotv * u;
%
%  Compute the normal component of the vector.
%
  normal = v - parallel;

  normal_norm = norm ( normal );

  if ( normal_norm == 0.0 )
    w = parallel;
    return
  end

  normal = normal / normal_norm;
%
%  Compute a second vector, lying in the plane, perpendicular
%  to V, and forming a right-handed system, as the cross product
%  of the first two vectors.
%
  normal2 = zeros ( 3, 1 );

  normal2(1) = u(2) * normal(3) - u(3) * normal(2);
  normal2(2) = u(3) * normal(1) - u(1) * normal(3);
  normal2(3) = u(1) * normal(2) - u(2) * normal(1);

  normal2_norm = norm ( normal2 );

  normal2 = normal2 / normal2_norm;
%
%  Rotate the normal component by the angle.
%
  rot = normal_norm * ( ...
      cos ( angle ) * normal ...
    + sin ( angle ) * normal2 );
%
%  The rotated vector is the parallel component plus the rotated component.
%
  w = parallel + rot;

  return
end
