function   [A1, A0, class_vec, B1, B0]=kak1(U)

%Calculates a special case (KAK1) 
%of Cartan's KAK Decomposition.
%For an input 4d unitary matrix U,
%it finds [A1, A0, class_vec, B1, B0]
%such that U = kron(A1,A0)*exp(i*M)*kron(B1,B0)
%where matrix M is a function of
%class_vec(1:4) 

%If 
%k=class_vec
%and 
%M=k(4)+sigxx*k(1)+sigyy*k(2)+sigzz*k(3)
%then
%kron(A1,A0)*exp(i*M)*kron(B1,B0) = U

global_declarations;

do_tests = false;

if(rows(U)!=4 | columns(U)!=4)
	error("input matrix for kak decomposition is not 4x4");
end
if(norm(U*U'-eye(4))>1e-8)
	error("input matrix for kak decomposition is not unitary");
end

%normalize U so that it has
%unit determinant. Store
%det(U)^(1/4) for later use
root4th = det(U)^(1/4);
if(abs(root4th-1)>1e-8)
	U = U/root4th;
end

U_bell = uni2ort'*U*uni2ort;

[L, D, R] = odo(U_bell);

if(do_tests)
	err1=norm(L*D*R'-U_bell)
	err2=norm(L*transpose(L)-eye(4))
	err3=norm(R*transpose(R)-eye(4))
	err4=norm(D*D'-eye(4))
end


class_vec = diag_to_class_vec(diag(D));
class_vec(4) = class_vec(4) + angle(root4th);

[A1,A0] = SO4_to_SU2xSU2(L);
[B1,B0] = SO4_to_SU2xSU2(R');

if(do_tests)
	det(A1)
	det(A0)
	det(A1)*det(A0)
	det(B1)
	det(B0)
	det(B1)*det(B0)
end