function 	vecs_out = diag_ckt_invar3(ckt_invar)

%Given a matrix ckt_invar which is known to be
%the LO-RHS invariant for a 2-qubit circuit with
%3 DCNOTs,
%find a (non-unique) set of defining vectors for the 
%DC-NOTs of a circuit with that invariant.


lam3 = trace(ckt_invar)/4;
lam3r = real(lam3);
lam3i = imag(lam3);
Del = ckt_invar - lam3*eye(4);
L3r = real(Gamma_rep((Del + Del')/2));
L3i = real(Gamma_rep((Del - Del')/(2i)));

xxo = -lam3r;
yyo = -lam3i;

[L, DX, DY, R] = simul_real_svd(L3r, L3i);

nu = [DX(1,1), DX(2,2), DX(3,3)];

mu = [DY(1,1), DY(2,2), DY(3,3)];

gx1 = L(1:3,1);
gy1 = L(1:3,2);
gz1 = L(1:3,3); 

gy0 = R(1:3,1);
gz0 = R(1:3,2);
gx0 = R(1:3,3);

%error_L3r = norm(L3r - nu(1)*gx1*gy0' - nu(2)*gy1*gz0' - nu(3)*gz1*gx0')
%error_L3i = norm(L3i - mu(1)*gx1*gy0' - mu(2)*gy1*gz0' - mu(3)*gz1*gx0')

if(mu(3)*nu(2)<0)
	xi=-1;
else
	xi=1;
end
xi2 = -xi;	%-xi*xi2=1

proot = sqrt(mu(3)^2+xxo^2);
sin_bet = xi*mu(3)/proot;
cos_bet = xi*xxo/proot;

%%%%%%%%%%%%%%%%%finding bet1, bet2

if(abs(cos_bet)>abs(sin_bet))
	s1_s2 = xi*xxo/cos_bet;
	c1_c2 = nu(3)/cos_bet;
	s1_c2 = -xi2*mu(1)/cos_bet;
	c1_s2 = -mu(2)/cos_bet;
else
	s1_s2 = xi*mu(3)/sin_bet;
	c1_c2 = yyo/sin_bet;
	s1_c2 = xi2*nu(2)/sin_bet;
	c1_s2 = nu(1)/sin_bet;
endif

c_sum = c1_c2 - s1_s2;
c_dif = c1_c2 + s1_s2;
s_sum = s1_c2 + c1_s2;
s_dif = s1_c2 - c1_s2;
bet_sum = atan2(s_sum,c_sum);
bet_dif = atan2(s_dif,c_dif);
bet1 = (bet_sum + bet_dif)/2;
bet2 = (bet_sum - bet_dif)/2;

%Must get sin(bet1) and sin(bet2) positive.
%At this point s1_s2 is positive, but
%it's possible that both
%sin(bet1) and sin(bet2) are negative

bad_signs1=0;
bad_signs2=0;
if(sin(bet1)<0)
	bet1=-bet1;
	bet2=-bet2;
	bad_signs1=1;
endif

%At this point
%it's possible that xi2*cos(bet2)<0
if(xi2*cos(bet2)<0)
	bet1=pi-bet1;
	bet2=pi-bet2;
	bad_signs2=1;
endif

%If both bad_signs1 and bad_signs2 are
%true but not both, do the following changes.
%The reason we don't do anything
%when they are both true is that
%we should do the changes twice,
%but that is equivalent to no changes
if(bad_signs1+bad_signs2==1)
	gx1=-gx1;
	gy1=-gy1;
	nu(1)=-nu(1);
	nu(2)=-nu(2);
	mu(1)=-mu(1);
	mu(2)=-mu(2);
endif


%%%%%%%%%%%%%%%%%%%%%%%%%tests

%sin1 = sin(bet1); 
%cos1 = cos(bet1);
%sin2 = sin(bet2);
%cos2 = cos(bet2);
%c2_xi2 = cos2*xi2;

%alt_mu = [-cos_bet*sin1*cos2*xi2,-cos_bet*cos1*sin2,xi*sin_bet*sin1*sin2];
%err_mu = norm(mu-alt_mu)
%alt_nu = [sin_bet*cos1*sin2,sin_bet*sin1*cos2*xi2,cos_bet*cos1*cos2];
%err_nu = norm(nu-alt_nu)
%alt_xxo = cos_bet*sin1*sin2*xi;
%err_xxo = abs(xxo-alt_xxo)
%alt_yyo = sin_bet*cos1*cos2;
%err_yyo = abs(yyo-alt_yyo)
%%%%%%%%%%%%%%%%%%%%%%%%%


	c0f = gx0;
	b0f = cos(bet2)*gx0 + sin(bet2)*gz0;
	a0f = cos(bet2 - xi2*bet1)*gx0 + sin(bet2 - xi2*bet1)*gz0;

	c1f = gx1;
	b1f = gz1;
	a1f = cos_bet*gx1 + sin_bet*gy1;

vecs_out=[c0f,b0f,a0f;c1f,b1f,a1f];