function 	vecs_out = diag_ckt_invar2(ckt_invar)

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

lam2r= real(trace(ckt_invar)/4);
Del = ckt_invar - lam2r*eye(4);
L2r = real(Gamma_rep((Del + Del')/2));
L2i = real(Gamma_rep((Del - Del')/(2i)));

c1_c0 = lam2r;

fy1 = get_unit_vec(L2r(1:3,1));
fy0 = get_unit_vec(L2r(1,1:3)');
s1_s0 = - fy1'*L2r*fy0;

hy1 = fy1;
hz1 = get_normal_unit_vec(hy1,[0;0;0]);
hx1 = cross(hy1,hz1);

hy0 = fy0;
hz0 = get_normal_unit_vec(hy0,[0;0;0]);
hx0 = cross(hy0,hz0);

M = [hz1'*L2i*hx0,hz1'*L2i*hz0; hx1'*L2i*hx0,hx1'*L2i*hz0];

[s1_c0, c1_s0, fx0, fz0, fx1, fz1]= diag_ckt_invar2_aux(M, hx0, hz0, hx1, hz1);


c_sum = c1_c0 - s1_s0;
c_dif = c1_c0 + s1_s0;
s_sum = s1_c0 + c1_s0;
s_dif = s1_c0 - c1_s0;
alp_sum = atan2(s_sum,c_sum);
alp_dif = atan2(s_dif,c_dif);
alp1 = (alp_sum + alp_dif)/2;
alp0 = (alp_sum - alp_dif)/2;

b0f= fx0;
a0f = fx0*cos(alp0)- fy0*sin(alp0);

b1f= fx1;
a1f = fx1*cos(alp1)- fy1*sin(alp1);

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