(* Michel Fodje's Minkwe simulation
translated from Python to Mathematica by John Reed
13 Nov 2013 *)

(* Set run time parameters, initialize arrays *)

spin=1/2;
phase=2 π spin;
spin2=2 spin;
trials=800000;

aliceDeg=ConstantArray[0,trials];
bobDeg=ConstantArray[0,trials];
aliceDet=ConstantArray[0,trials];
bobDet=ConstantArray[0,trials];

nPP=ConstantArray[0,361];
nNN=ConstantArray[0,361];
nPN=ConstantArray[0,361];
nNP=ConstantArray[0,361];
nA=ConstantArray[0,361];
nB=ConstantArray[0,361];

(* Detector test function *)

test[angle_,e_,λ_]:=Module[{c,out},
c=(-1)^spin2 Cos[spin2(angle-e)];
If[λ≥Abs[c],out=0,out=Sign[c]];
out]

(* Generate particle data  *)

Do[
eVector=RandomReal[{0,2 π}];
λ=1/2 Sin[RandomReal[{0,π/2}]]^2;
eLeft=RandomReal[{0,2 π}];
eRight=eLeft+2 π spin;
aliceAngle=RandomReal[{0,2 π}];
aliceDeg[[i]]=aliceAngle/Degree;
bobAngle=RandomReal[{0,2 π}];
bobDeg[[i]]=bobAngle/Degree;
aliceDet[[i]]=test[aliceAngle,eLeft,λ];
bobDet[[i]]=test[bobAngle,eRight,λ],
{i,trials}]

(* statistical analysis of particle data  *)

Do[
θ=Ceiling[Abs[(aliceDeg[[i]]-bobDeg[[i]])]];
aliceD=aliceDet[[i]];bobD=bobDet[[i]]; 
If[aliceD==1,nA[[θ]]++];
If[bobD==1,nB[[θ]]++];
If[aliceD==1&&bobD==1,nPP[[θ]]++];
If[aliceD==1&&bobD==-1,nPN[[θ]]++];
If[aliceD==-1&&bobD==1,nNP[[θ]]++];
If[aliceD==-1&&bobD==-1,nNN[[θ]]++],
{i,trials}]

(* Calculate mean values and plot *)

pPP=0; pPN=0; pNP=0;pNN=0;

mean=ConstantArray[0,361];

Do[
sum=nPP[[i]]+nPN[[i]]+nNP[[i]]+nNN[[i]];
If[sum==0,Goto[jump],
{pPP=nPP[[i]]/sum;
       pNP=nNP[[i]]/sum;
       pPN=nPN[[i]]/sum;
       pNN=nNN[[i]]/sum;
mean[[i]]=pPP+pNN-pPN-pNP}];
Label[jump],
{i,361}]

simulation=ListPlot[mean]

(mean[[24]]+mean[[23]])/2//N

(mean[[68]]+mean[[69]])/2//N

cos=Plot[-Cos[x Degree],{x,0,360},PlotStyle->{Red,Thick}];

(* Compare mean values with Cosine 
the Cosine curve is off by one degree compared to simulation *)

Show[simulation,cos]