% computes the angle and axis for a given rotation matrix

function [axis,varargout] = angleaxis(R)

alpha = acos((trace(R)-1)/2);
a = [R(3,2)-R(2,3); R(1,3)-R(3,1); R(2,1)-R(1,2)];
nn = norm(a);
if nn > 1e-10
 a=a./norm(a);
end

if nargout == 2 
 axis = a;
 varargout(1) = {alpha};
 varargout(2) = {axis};
else
 axis = a*alpha;
end; 
