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



% b0   a0         b0f   a0f   [U0]

% |    |    -->   |     |           

% b1   a1         b1f   a1f   [U1]

%where angle(b0,a0)->angle(b1f,a1f)

%and   angle(b1,a1)->angle(b0f,a0f)



check_dcnots(vecs_in);



b0=vecs_in(1:3,1);a0=vecs_in(1:3,2);

b1=vecs_in(4:6,1);a1=vecs_in(4:6,2);



cos_alp0 = a0'*b0;

sin_alp0 = sqrt(1-cos_alp0^2);

cos_alp1 = a1'*b1;

sin_alp1 = sqrt(1-cos_alp1^2);



fx0 = b0;

fy0 = (cos_alp0*fx0 - a0)/sin_alp0;

fz0 = cross(fx0, fy0);

fx1 = b1;

fy1 = (cos_alp1*fx1 - a1)/sin_alp1;

fz1 = cross(fx1, fy1);



b0f = fz0;

b1f = fz1;

a0f = cos_alp1*fz0 + sin_alp1*fy0;

a1f = cos_alp0*fz1 + sin_alp0*fy1;

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



sin_halp0 = sqrt((1-cos_alp0)/2);

cos_halp0 = sqrt((1+cos_alp0)/2);

sin_halp1 = sqrt((1-cos_alp1)/2);

cos_halp1 = sqrt((1+cos_alp1)/2);





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



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