%test breach()

%conventions:
%L, R refer to left and right 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

global_defs;

example = 1
switch example
	case (1)
		qL0= get_unit_vec(randn(3,1));
		pL0= get_unit_vec(randn(3,1));
		pR0= get_unit_vec(randn(3,1));
		qR0= get_unit_vec(randn(3,1));

		qL1= get_unit_vec(randn(3,1));
		pL1= get_unit_vec(randn(3,1));
		pR1= get_unit_vec(randn(3,1));
		qR1= get_unit_vec(randn(3,1));
	otherwise
		error("example number is out of range");
end
vecs_in=[qL0,pL0,pR0,qR0;qL1,pL1,pR1,qR1];
[UL1, UL0, vecs_out, UR1, UR0] = breach(vecs_in);
qL0f=vecs_out(1:3,1);pL0f=vecs_out(1:3,2);pR0f=vecs_out(1:3,3);qR0f=vecs_out(1:3,4);
qL1f=vecs_out(4:6,1);pL1f=vecs_out(4:6,2);pR1f=vecs_out(4:6,3);qR1f=vecs_out(4:6,4);

error1 = norm( kron(UL1,UL0)*dr11(qL1f,qL0f)*dr11(pL1f,pL0f)*dr11(pR1f,pR0f)*dr11(qR1f,qR0f)*kron(UR1,UR0) -dr11(qL1,qL0)*dr11(pL1,pL0)*dr11(pR1,pR0)*dr11(qR1,qR0) )

error2 = norm(pL1f-pR1f)