%Tests the Deflation Identity (quant-ph/0509111).

%_p = positive subscript +
%_n = negative subscript -
global_defs;

tL = 2*pi*randn();
%tL=pi/2
tR =  2*pi*randn();
bet1 =  2*pi*randn();
bet0 =  2*pi*randn();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%constructing U_in

cnot_tL = [id2,ze2;ze2, cos(tL)*id2 + i*sin(tL)*sigz];
cnot_tR = [id2,ze2;ze2, cos(tR)*id2 + i*sin(tR)*sigz];
exp_bet0 = cos(bet0)*id2 + i*sin(bet0)*sigy;
exp_bet1 = cos(bet1)*id2 + i*sin(bet1)*sigy;
U_in = cnot_tL*kron(exp_bet1, exp_bet0)*cnot_tR;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%constructing U_out

cos_xi_p = cos((tL+tR)/2)*cos(-bet1 + bet0);
cos_xi_n = cos((tL+tR)/2)*cos(-bet1 - bet0);
sin_xi_p = cos((tL-tR)/2)*sin(-bet1 + bet0);
sin_xi_n = cos((tL-tR)/2)*sin(-bet1 - bet0);

[cos_xi_p, sin_xi_p, p_p, xi_p] = normalize_cs_pair(cos_xi_p, sin_xi_p);
[cos_xi_n, sin_xi_n, p_n, xi_n] = normalize_cs_pair(cos_xi_n, sin_xi_n);

cos_eta_p =  sin( (tR-tL)/2)*sin(-bet1 + bet0);
cos_eta_n = -sin( (tR-tL)/2)*sin(-bet1 - bet0);
sin_eta_p =  sin(-(tL+tR)/2)*cos(-bet1 + bet0);
sin_eta_n = -sin(-(tL+tR)/2)*cos(-bet1 - bet0);

[cos_eta_p, sin_eta_p, q_p, eta_p] = normalize_cs_pair(cos_eta_p, sin_eta_p);
[cos_eta_n, sin_eta_n, q_n, eta_n] = normalize_cs_pair(cos_eta_n, sin_eta_n);

mu_p = atan2(q_p, p_p);
mu_n = atan2(q_n, p_n);

gam0L = (      eta_p - eta_n + xi_p - xi_n)/4;
gam1L = ( pi - eta_p - eta_n - xi_p - xi_n)/4;
gam0R = (    - eta_p + eta_n + xi_p - xi_n)/4;
gam1R = (-pi + eta_p + eta_n - xi_p - xi_n)/4;

exp_tL = cos(tL/2)*id2 + i*sin(tL/2)*sigz;
exp_tR = cos(tR/2)*id2 + i*sin(tR/2)*sigz;
exp0L = cos(gam0L)*id2 +i*sin(gam0L)*sigy;
exp1L = cos(gam1L)*id2 +i*sin(gam1L)*sigy;
exp0R = cos(gam0R)*id2 +i*sin(gam0R)*sigy;
exp1R = cos(gam1R)*id2 +i*sin(gam1R)*sigy;


mid_x = cos((mu_p+mu_n)/2)*eye(4) + i*sin((mu_p+mu_n)/2)*sigxx;
mid_z = cos((mu_p-mu_n)/2)*eye(4) + i*sin((mu_p-mu_n)/2)*sigzz;

U_out= kron(exp1L, exp_tL*exp0L)*mid_x*mid_z*kron(exp1R, exp0R*exp_tR);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
error = norm(U_out-U_in)