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