function [UL1, UL0, vecs_out, UR1, UR0] = breach(vecs_in)

%This function takes a 2-qubit circuit
%with 4 DC-NOTs and it returns
%an equivalent circuit which has a breach
%in the lower wire, between its two middle DC-NOTs.
%conventions:
%L, R refer to left and right hand side
%0=top and 1=bottom qubit
%f denotes final quantities.

% qL0   pL0   pR0   qR0       [UL0]   qL0f  pL0f  pR0f  qR0f  [UR0]
% |     |     |     |    -->           |     |     |     |  
% qL1   pL1   pR1   qR1       [UL1]   qL1f  pL1f  pR1f  qR1f  [UR1]

% where pL1f = pR1f

%IMP: qL1, qL0, pL1, pL0, pR1, pR0, qR1, qR0 should all be unit vectors
check_dcnots(vecs_in);

qL0=vecs_in(1:3,1);pL0=vecs_in(1:3,2);pR0=vecs_in(1:3,3);qR0=vecs_in(1:3,4);
qL1=vecs_in(4:6,1);pL1=vecs_in(4:6,2);pR1=vecs_in(4:6,3);qR1=vecs_in(4:6,4);

cross_qL1_pL1 = cross(qL1, pL1);
cross_qR1_pR1 = cross(qR1, pR1);

vecL= cross(cross_qL1_pL1, pL1);
vecR= cross(cross_qR1_pR1, pR1);
t1=get_normal_unit_vec(vecL,vecR);

hxL1= pL1; 
norm_cross_qL1_pL1 = norm(cross_qL1_pL1);
if(norm_cross_qL1_pL1>1e-10)
	hzL1 = cross_qL1_pL1/norm_cross_qL1_pL1;
else
	error("qL1==pL1 so breach already exixts")
end
hyL1= cross(hzL1, hxL1);
cos_the_L1= qL1'*hxL1; sin_the_L1 = -qL1'*hyL1;
cos_phi_L1= t1'*hzL1; sin_phi_L1 = t1'*hxL1;

hxR1= pR1;
norm_cross_qR1_pR1 = norm(cross_qR1_pR1);
if(norm_cross_qR1_pR1>1e-10)
	hzR1 = cross_qR1_pR1/norm_cross_qR1_pR1;
else
	error("qR1==pR1 so breach already exixts")
end
hyR1= cross(hzR1, hxR1);
cos_the_R1= qR1'*hxR1; sin_the_R1 = -qR1'*hyR1;
cos_phi_R1= t1'*hzR1; sin_phi_R1 = t1'*hxR1;


vecL= (cos_phi_L1*(qL0'*pL0))*cross(qL0,pL0) - (sin_the_L1*cos_the_L1*sin_phi_L1)*pL0;
vecR= (cos_phi_R1*(qR0'*pR0))*cross(qR0,pR0) - (sin_the_R1*cos_the_R1*sin_phi_R1)*pR0;
t0=get_normal_unit_vec(vecL,vecR);


vin = [t0,pR0,qR0;t1,pR1,qR1];
[vout, UR1, UR0] = dr_3to2(vin,1);
pR0f=vout(1:3,1);qR0f=vout(1:3,2);
pR1f=vout(4:6,1);qR1f=vout(4:6,2);

vin = [t0,pL0,qL0;t1,pL1,qL1];
[vout, UL1, UL0] = dr_3to2(vin,1);
pL0f=vout(1:3,1);qL0f=vout(1:3,2);
pL1f=vout(4:6,1);qL1f=vout(4:6,2);


UL0 = UL0';
UL1 = UL1';

vecs_out=[qL0f,pL0f,pR0f,qR0f;qL1f,pL1f,pR1f,qR1f];