function [vecs_out, U1, U0] = dr_3to2(vecs_in, persistence)

%dr= dressed cnot
%%%%%%%%%%%%%%%%%%
% c0   b0   a0         b0f   a0f   [U0]
% |    |    |    -->   |     |           *exp(i*pi/4)
% c1   b1   a1         b1f   a1f   [U1]

% where U0, U1 are SU(2) matrices.
% Must satisfy [a1,b1,b1)'*c1=0 (right angle at b1).
% If persistence==true, b1f=c1 (i.e.,``c1 is persistent"),
% When persistence==true, checks persit_cond 
% (a sufficient condition for persistence) before proceding.
%%%%%%%%%%%%%%%%%%

c0=vecs_in(1:3,1);b0=vecs_in(1:3,2);a0=vecs_in(1:3,3);
c1=vecs_in(4:6,1);b1=vecs_in(4:6,2);a1=vecs_in(4:6,3);

kx1= b1;
kz1= get_unit_vec(cross(a1, b1));
ky1= cross(kz1, kx1);
cos_lam1= a1'*kx1; sin_lam1 = -a1'*ky1;
cos_phi1= c1'*kz1; sin_phi1 = c1'*kx1;

%%%%%%%%%%%%%%%%%%shorthands

cross_ab1= cross(a1,b1);
cross_bc1 = cross(b1,c1);
cross_abc1= cross(cross_ab1,c1);
cross_abcc1 = cross(cross_abc1,c1);
cross_ab1 = cross(cross_ab1,b1);
cross_abbc1 = cross(cross_ab1,c1);
cross_abbcc1 = cross(cross_abbc1,c1);
cross_bcc1 = cross(cross_bc1,c1);

cross_ab0= cross(a0,b0);
cross_bc0 = cross(b0,c0);
cross_abc0= cross(cross_ab0,c0);
cross_abcc0 = cross(cross_abc0,c0);
cross_ab0f = cross(cross_ab0,b0);
cross_abbc0 = cross(cross_ab0f,c0);
cross_abbcc0 = cross(cross_abbc0,c0);
cross_bcc0 = cross(cross_bc0,c0);

vol_0 = cross_ab0'*c0;
vol_1 = cross_ab1'*c1;

%%%%%%%%%%%%%%%%%%special _conditions
rangle_cond = cross_ab1'*c1;
persist_cond = (cos_phi1*(a0'*b0)*cross_ab0 - sin_lam1*cos_lam1*sin_phi1*b0)'*c0;

if(abs(rangle_cond)>1e-10)
	error("right angle condition not satisfied");
end
if(persistence==1 & (abs(persist_cond)>1e-10))
	error("persistence condition not satisfied");
end
%%%%%%%%%%%%%%%%%%

c1_c0 = (a0'*b0)*(b0'*c0)*sin_lam1*cos_phi1 + vol_0*cos_lam1*sin_phi1;

if(sin_lam1 > 1e-10)
	fy1 = cross_abbc1/sin_lam1;
else
	error("can't determine fy1");
end

h = (sin_lam1*sin_phi1*(a0'*b0))*cross_bcc0  -(cos_lam1*cos_phi1)*cross_abcc0  +(sin_lam1*(cross_ab0f'*c0))*c0;
norm_h = norm(h);
if(norm_h>1e-10)
	fy0 = h/norm_h;
else
	error("can't determine fy0");
end

s1_s0 =  norm_h;

v1 = -(cos_lam1*(a0'*b0))*c0  +(sin_lam1*cos_phi1*(a0'*b0))*cross_bc0   +(cos_lam1*sin_phi1)*cross_abc0;
norm_v1 = norm(v1);

v2=(sin_lam1*sin_phi1*(a0'*b0)*(b0'*c0) -cos_lam1*cos_phi1*vol_0)*c0 +sin_lam1*cross_abbcc0;
norm_v2 = norm(v2);

if(persistence==1)
	fx1 = c1; 
	fz1 = ky1; 
	
	fz0 = v1/norm_v1; 
	fx0 = v2/norm_v2; 
	
	c1_s0 = norm_v1;
	s1_c0 = norm_v2;
	%(fx1,fy1,fz1) always right handed 3d basis
	%but (fx0,fy0,fz0) may not be right handed at this point
	hand0 = cross(fz0,fx0)'*fy0;
	if(hand0<0)
		fy0 = -fy0;
		s1_s0 = -s1_s0;
	end
else %(persistence==0)
	L2i = c1*v1' + ky1*v2';
	hx1=c1;
	hy1=fy1;
	hz1=ky1;
	
	hy0 = fy0;
	hz0 = get_normal_unit_vec(fy0,[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);

endif

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];  

%%%%%%%%%%%%%%%%%% now find U1 and U0

[U1,U0]= factor_SU2pow2_matrix(
exp(-i*pi/4)*dr11(a1f,a0f)*dr11(b1f,b0f)*dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0));