function [theta] = rotangle(H)
%
if det(H)<0
   error('The matrix is not a product of rotations')
end
n = size(H,1);
theta = [];
for i = 1 : n - 1
   for j = i + 1 : n
      theta_k = atan2(-H(j,i),H(i,i));
      theta = [theta_k, theta];
      c = cos(theta_k); s = sin(theta_k);
      R = eye(n);
      R(i,i) = c; R(j,j) = c;
      R(i,j) = -s; R(j,i) = s;
      H = R * H;
   end
end

