function [L, D, R] = odo(U)
%odo=orthogonal diagonal orthogonal
%Given a unitary matrix U,
%it finds special (det=1) orthogonal matrices L and R
%and a diagonal unitary matrix D
%such that L*D*R'=U
dim_U= rows(U);
if(norm(U*U'-eye(dim_U))>1e-8)
error("input matrix for odo decomposition is not unitary");
end
X = (U + conj(U))/2;%IMP: this is NOT equal to (U+U')/2
Y = (U - conj(U))/(2*i);
% U = X + i*Y; X, Y both real
[L, DX, DY, R] = simul_real_svd(X,Y);
D = DX + i*DY;