Ultimate local-realistic simulation of two level entangled state event by event.
 By Fred Diether with 2D vectors, also with Joy’s hidden variable and modified
 epr-simple HV by Michel Fodje.  Some parts by John Reed.  6-10-2021

<<Quaternions`;
β0=Quaternion[1,0,0,0];
β1=Quaternion[0,1,0,0];
β2=Quaternion[0,0,1,0];
β3=Quaternion[0,0,0,1];
Qcoordinates={β1,β2,β3};
m=5000000;
lista=Table[{0,0,0,0,0},m];
list1=Table[{0,0,0,0,0},m];
list2=Table[{0,0,0,0,0},m];
list3=Table[{0,0,0,0,0},m];
list4=Table[{0,0,0,0,0},m];
list5=Table[{0,0,0,0,0},m];
list6=Table[{0,0,0,0,0},m];
Do[vectorA=Flatten[{RandomPoint[Circle[]],0}];
vectorB=Flatten[{RandomPoint[Circle[]],0}];
vectorS=Flatten[{RandomPoint[Circle[]],0}]; (*Singlet spin vector*)
s1=ToSphericalCoordinates[vectorS][[3]];
λ2=(1/4)Cos[((s1*90)/2)]^2;  (*hidden variable*)
λ=RandomChoice[{-1,1}];  (*hidden variable*)
Da=vectorA.Qcoordinates;  (*Convert to quaternion coordinates*)
Db=vectorB.Qcoordinates;
Ls1=λ*vectorS.Qcoordinates; (*singlet spin quaternion*)
Ls2=-Ls1;
If[Abs[Re[Da**Ls1]]>λ2,A=Sign[Re[Da**Ls1]],A=Null];
If[(π/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤(π/2), Aa=-1,Aa=Null];
If[((7π)/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤((2π)/1), AA=1,AA=Null];
If[((3π)/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤((3π)/2), AA1=1,AA1=Null];
If[((0π)/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤(π/4), AA2=1,AA2=Null];
If[((0π)/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤(π/4), AA3=-1,AA3=Null];
If[((0π)/4)≤ArcTan[vectorA[[1]],vectorA[[2]]]+π≤(π/4), AA4=1,AA4=Null];
If[Abs[Re[Ls2**Db]]>λ2,B=Sign[Re[Db**Ls2]],B=Null ];
If[((0π)/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤(π/4), Bb=1,Bb=Null];
If[((0π)/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤(π/4), BB=-1,BB=Null];
If[((0π)/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤(π/4), BB1=1,BB1=Null];
If[((3π)/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤((3π)/2), BB2=1,BB2=Null];
If[(π/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤(π/2), BB3=1,BB3=Null];
If[((7π)/4)≤ArcTan[vectorB[[1]],vectorB[[2]]]+π≤((2π)/1), BB4=-1,BB4=Null];
AT=ArcTan[vectorA[[1]],vectorA[[2]]];
BT=ArcTan[vectorB[[1]],vectorB[[2]]];
lista[[i]]={AT,BT,A,B,i};
list1[[i]]={AT,BT,AA1,BB1,3m+i};
list2[[i]]={AT,BT,AA2,BB2,4m+i};
list3[[i]]={AT,BT,AA3,BB3,5m+i};
list4[[i]]={AT,BT,AA4,BB4,6m+i};
list6[[i]]={AT,BT,AA,BB,m+i};
list5[[i]]={AT,BT,Aa,Bb,2m+i},{i,m}]

listb=DeleteCases[lista,{_,_,Null, Null,_}];
listc=DeleteCases[listb,{_,_,_,Null,_}];
listd=DeleteCases[listc,{_,_,Null,_,_}];
neg=DeleteCases[list6,{_,_,Null,Null,_}];
neg1=DeleteCases[neg,{_,_,_,Null,_}];
neg2=DeleteCases[neg1,{_,_,Null,_,_}];
neg3=DeleteCases[list3,{_,_,Null,Null,_}];
neg4=DeleteCases[neg3,{_,_,_,Null,_}];
neg5=DeleteCases[neg4,{_,_,Null,_,_}];
neg6=DeleteCases[list4,{_,_,Null,Null,_}];
neg7=DeleteCases[neg6,{_,_,_,Null,_}];
neg8=DeleteCases[neg7,{_,_,Null,_,_}];
negneg=DeleteCases[list5,{_,_,Null, Null,_}];
negneg1=DeleteCases[negneg,{_,_,_,Null,_}];
negneg2=DeleteCases[negneg1,{_,_,Null,_,_}];
pos1=DeleteCases[list1,{_,_,Null, Null,_}];
pos2=DeleteCases[pos1,{_,_,_,Null,_}];
pos3=DeleteCases[pos2,{_,_,Null, _,_}];
pos4=DeleteCases[list2,{_,_,Null, Null,_}];
pos5=DeleteCases[pos4,{_,_,_, Null,_}];
pos6=DeleteCases[pos5,{_,_,Null, _,_}];

negneg3=Take[Catenate[{listd, neg2,negneg2,pos3,pos6,neg5,neg8}],m];
trials=Length[negneg3];
trialDeg=720;
a1=ConstantArray[0,trials];
b1=ConstantArray[0,trials];
A1=ConstantArray[0,trials];
B1=ConstantArray[0,trials];
nPP=ConstantArray[0,trialDeg];
nNN=ConstantArray[0,trialDeg];
nPN=ConstantArray[0,trialDeg];
nNP=ConstantArray[0,trialDeg];
nAP=ConstantArray[0,trialDeg];
nBP=ConstantArray[0,trialDeg];
nAN=ConstantArray[0,trialDeg];
nBN=ConstantArray[0,trialDeg];
a1=negneg3[[All,1]];b1=negneg3[[All,2]];A1=negneg3[[All,3]];B1=negneg3[[All,4]];

Do[a2=a1[[j]];
b2=b1[[j]];
t1=Round[(a2-b2)/Degree+360];
aliceD=A1[[j]];bobD=B1[[j]];
If[aliceD==1,nAP[[t1]]++];
If[bobD==1,nBP[[t1]]++];
If[aliceD==-1,nAN[[t1]]++];
If[bobD==-1,nBN[[t1]]++];
If[aliceD==1&&bobD==1,nPP[[t1]]++];
If[aliceD==1&&bobD==-1,nPN[[t1]]++];
If[aliceD==-1&&bobD==1,nNP[[t1]]++];
If[aliceD==-1&&bobD==-1,nNN[[t1]]++],{j,trials}]
pPP=0; pPN=0; pNP=0;pNN=0; 
mean=ConstantArray[0,trialDeg];
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,trialDeg}]
sim=ListPlot[(mean),PlotMarkers->{Automatic,Tiny}];
negcos=Plot[-Cos[x Degree],{x,0,720},PlotStyle->{Magenta},Ticks->{{{0,-360°},{90,-270°},{180,-180°},{270,-90°},{360,0°},{450,90°},{540,180°},{630,270°},{720,360°}},Automatic},GridLines->Automatic];
Show[negcos, sim]

AveA=N[Sum[A1[[i]], {i,trials}]/trials];
AveB=N[Sum[B1[[i]], {i, trials}]/trials];
Print["AveA = ",AveA]
Print["AveB = ",AveB]
totAB=Sum[nPP[[i]]+nNN[[i]]+nPN[[i]]+nNP[[i]], {i, trialDeg}]
CHSH=Abs[N[mean[[45]]]+N[mean[[45]]]-N[mean[[135]]]+N[mean[[45]]]];
Print["CHSH = ", CHSH]