function	[vecs_out, U1, U0] = pass3(vecs_in, perp)

% e0  c0  b0  a0       c0f  e0f  b0f   a0f    [U0]
% |   |   |   |   -->  |    |    |     |
% |   c1  b1  a1       c1f  |    b1f   a1f    [U1]
% |                         |
% e2                        e2f               [1]

%Both LHS and RHS circuits contain
%4 DC-NOTs acting on 3 qubits.
%One of the bit0-bit1 DC-NOTs 
%"passes through" the only bit0-bit2 DC-NOT.
%We check that 
%(1)b1 perp span(a1,c1)
%(2){c0,b0,a0} is a dependent set
%(3)-(cross(a0,b0)'*cross(b0,c0))*(b0'*c0)>0
%We find a c1f that is 
%perpendicular to the input perp.


check_dcnots(vecs_in);

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

cross_ab1 = cross(a1,b1);
cross_bc1 = cross(b1,c1);
%cross_abc1= cross(cross_ab1,c1);
%cross_abcc1 = cross(cross_abc1,c1);
%cross_abb1 = cross(cross_ab1,b1);
%cross_abbc1 = cross(cross_abb1,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_abb0 = cross(cross_ab0,b0);
cross_abbc0 = cross(cross_abb0,c0);
%cross_abbcc0 = cross(cross_abbc0,c0);
cross_bcc0 = cross(cross_bc0,c0);

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

ab1 = a1'*b1;
bc1 = b1'*c1;
ac1 = a1'*c1;
ab0 = a0'*b0;
bc0 = b0'*c0;

eps = 1e-10;
if(norm(ab1)>eps | norm(bc1)>eps | vol_0>eps)
	error("first 3 dcnots are not in diagonalized form");
end
if(-cross_ab0'*cross_bc0*bc0<0)
	error("wrong orientation");
end

sin_bet2 = norm(cross_bc0);
eta = norm(cross_abbc0);
xxo = ac1*(cross_abb0'*c0);
yyo = ab0*bc0*vol_1;

mu = [-ac1*eta, -ab0*ac1*sin_bet2, (cross_abb0'*c0)*vol_1];

nu = [ab0*vol_1*sin_bet2, vol_1*eta, ab0*bc0*ac1];

gx0= c0;
gy0= cross_bc0/sin_bet2;
gz0= -cross_bcc0/sin_bet2;

gx1= c1;
gy1= cross_bc1;
gz1= b1;

Mmu = mu(1)*gx1*gy0' +mu(2)*gy1*gz0' +mu(3)*gz1*gx0';
Mnu = nu(1)*gx1*gy0' +nu(2)*gy1*gz0' +nu(3)*gz1*gx0';

d0 = e0;
Mmu_d0 = Mmu*d0;
Mnu_d0 = Mnu*d0;
d1=get_normal_unit_vec(Mmu_d0,perp);
Mmu_t_d1 = Mmu'*d1;
Mnu_t_d1 = Mnu'*d1;
x0 = cross(Mmu_t_d1,d0);
x1 = cross(Mmu_d0,d1);
%y0 = cross(Mnu_t_d1,d0);
%y1 = cross(Mnu_d0,d1);
norm_x0 = norm(x0);
norm_x1 = norm(x1);	

c1_c0 = d1'*Mnu_d0;

fy0 = d0+x0/xxo;
mag0 = norm(fy0);
fy0 = fy0/mag0;

fy1 = d1+x1/xxo;
mag1= norm(fy1);
fy1 = fy1/mag1;

s1_s0 = xxo*mag1*mag0;

hz1 = d1-x1*xxo/norm_x1^2;
mag = norm(hz1);
hz1 = hz1/mag;

hz0 = d0-x0*xxo/norm_x0^2;
mag = norm(hz0);
hz0 = hz0/mag;

if(xxo>0)
	hx0 = cross(x0,d0)/norm_x0;
	hx1 = cross(x1,d1)/norm_x1;
else
	hx0 = -cross(x0,d0)/norm_x0;
	hx1 = -cross(x1,d1)/norm_x1;
end

rho0= sqrt(1 + (xxo/norm_x0)^2);
rho1= sqrt(1 + (xxo/norm_x1)^2);
M = [Mnu_t_d1'*cross(Mmu_t_d1,d0)*rho1/norm_x0,-yyo*rho1*rho0
;-(Mmu_d0)'*Mmu*Mmu_t_d1/(norm_x0*norm_x1),(Mnu_d0)'*cross(Mmu_d0,d1)*rho0/norm_x1];
if(xxo<0)
	M(1,1)=-M(1,1);
	M(2,2)=-M(2,2);
end

[s1_c0, c1_s0, fx0, fz0, fx1, fz1]= diag_ckt_invar2_aux(M, hx0, hz0, hx1, hz1);

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

c0f=e0;
c1f=d1;
e0f=e0;
e2f=e2;
vecs_out=[c0f,e0f,b0f,a0f;c1f,e2f,b1f,a1f];

%%%%%%%%%%%%%%%%%% now find U1 and U0
[U1,U0] = factor_SU2pow2_matrix(
dr11(a1f,a0f)*dr11(b1f,b0f)*dr11(d1,d0)*dr11(c1,c0)*dr11(b1,b0)*dr11(a1,a0));