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

% c0   b0   a0         a0f   [U0]
% |    |    |    -->   |           
% c1   b1   a1         a1f   [U1]

%Let {sx,sy,sz} be the standard basis. Then 
%a0=a1=cos(lam)sx -sin(lam)sy
%b0=b1=sx
%c0=cos(phi)sz +sin(phi)sx
%c1=cos(phi)sz -sin(phi)sx

global_declarations;

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

c_lam = a0(1);
s_lam = -a0(2);
c_phi = c0(3);
s_phi = c0(1);

a0f = c_lam*c0 + s_lam*[0;1;0];
a1f = c_lam*c1 + s_lam*[0;1;0];

vecs_out = [a0f;a1f];

%%%%%%%%%%%%%%%%%% now find U1 and U0
phi = atan2(s_phi, c_phi);
lam = atan2(s_lam, c_lam);


alp = (pi/2 - phi)/2;
c_alp = cos(alp);
s_alp = sin(alp);

bet = (pi - 2*alp)/2;
c_bet = cos(bet);
s_bet = sin(bet);

U0 = (c_alp*sigx + s_alp*sigz)*(c_alp*sig(a0) + s_alp*sigz);
U1 = (c_bet*sigx + s_bet*sigz)*(c_bet*sig(a0) + s_bet*sigz);
 

%[U1,U0] = factor_SU2pow2_matrix(
%dr11(a1f,a0f)*dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0));