function  [vecs_out, U1, U0] = swap_t3(vecs_in)



% c0   b0   a0         a0f   [U0]

% |    |    |    -->   |           

% c1   b1   a1         a1f   [U1]

%where a0 is perp to span(b0, c0)

%and   a1 is perp to span(b1, c1)



check_dcnots(vecs_in);



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



a0f = a0;

a1f = a1;



vecs_out = [a0f;a1f];



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



cross_ab0 = cross(a0,b0);

eta0 = norm(cross_ab0);

alp0 = atan2(-c0'*cross_ab0/eta0^2, c0'*b0);

sin_halp0 = sin(alp0/2);

cos_halp0 = cos(alp0/2);



cross_ab1 = cross(a1,b1);

eta1 = norm(cross_ab1);

alp1 = atan2(-c1'*cross_ab1/eta1^2, c1'*b1);

sin_halp1 = sin(alp1/2);

cos_halp1 = cos(alp1/2);



U0 = (cos_halp0*eye(2) + i*sin_halp0*sig(a0))*(cos_halp1*eye(2) - i*sin_halp1*sig(b0));



U1 = (cos_halp1*eye(2) + i*sin_halp1*sig(a1))*(cos_halp0*eye(2) - i*sin_halp0*sig(b1));

 



%[U1,U0] = factor_SU2pow2_matrix(

%dr11(a1f,a0f)*dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0));