Simulation Based on Michel Fodje's epr-simple simulation translated from 
Python to Mathematica by John Reed 13 Nov 2013 and Quaternions
Modified by Fred Diether for Completely Local-Realistic Sep 2021
Some parts by Bill Nelson.  Includes Joy’s (S^3) Quaternion Model.

Load Quaternion Package, Set Run Time Parameters, Initialize Arrays and Tables

<<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=10000;
trialDeg=721;
Ls1=ConstantArray[0,m];
Ls2=ConstantArray[0,m];
λ1=ConstantArray[0,m];
λ2=ConstantArray[0,m];
λ3=ConstantArray[0,m];
qa=ConstantArray[0,m];
aq1=ConstantArray[0,m];
qa1=ConstantArray[0,m];
qb=ConstantArray[0,m];
bq1=ConstantArray[0,m];
qb1=ConstantArray[0,m];
aa1=ConstantArray[0,m];
bb1=ConstantArray[0,m];
outA1=Table[{0,0,0},m];
outA2=Table[{0,0,0},m];
outB1=Table[{0,0,0},m];
outB2=Table[{0,0,0},m];
a2=ConstantArray[0,m];
b2=ConstantArray[0,m];
A=ConstantArray[0, m];
B=ConstantArray[0, m];
ϕ=3; β=0.32; ξ=0.939;  (*Adustable parameters for fine tuning*)

Generating Particle Data with Three Independent Do-Loops

Do[θ=RandomReal[{-179,180}];  (*Singlet vector angle*)   (*Hidden Variable*)
λ1[[i]]=β(Cos[(θ/ϕ) ]^2);
λ2[[i]]=ξ(Cos[(θ/ϕ) ]^2);
λ3[[i]]=Sign[θ];
θθ=N[Flatten[{FromPolarCoordinates[{1,θ *π/180}],0}]];Ls1[[i]]=λ3[[i]]*θθ.Qcoordinates;  (*Convert to quaternion coordinates*)
Ls2[[i]]=-λ3[[i]]*θθ.Qcoordinates,{i,m}]

Do[a=RandomChoice[{0,90}];  (*Detector vector angle 1 degree increments*)
aa1[[i]]=a;
aa=N[Flatten[{FromPolarCoordinates[{1,a *π/180}],0}]];
Da=aa.Qcoordinates;  (*Convert to quaternion coordinates*)
qa=Da**Ls1[[i]];
qa1[[i]]=qa;
aq=-(-Da**Ls1[[i]]);
aq1[[i]]=aq;
If[Abs[Re[qa]]>λ1[[i]],Aa1=Sign[Re[qa]],Aa1=Re[aq]];
outA1[[i]]={a,Aa1,qa},{i,m}]
Do[If[Abs[Re[qa1[[i]]]]<λ2[[i]],Aa2=-Sign[qa1[[i]][[4]]],Aa2=Re[aq1[[i]]]];outA2[[i]]={aa1[[i]],Aa2,qa1[[i]]},{i,m}]
outA=Catenate[{outA1,outA2}];
a2=outA[[All,1]];
A=outA[[All,2]];

Do[b=RandomChoice[{45,135}];  (*Detector vector angle 1 degree increments*)
bb1[[i]]=b;
bb=N[Flatten[{FromPolarCoordinates[{1,b *π/180}],0}]];
Db=bb.Qcoordinates;  (*Convert to quaternion coordinates*)
qb=Ls2[[i]]**Db;
qb1[[i]]=qb;
bq=-(-Ls2[[i]]**Db);
bq1[[i]]=bq;
If[Abs[Re[qb]]>λ1[[i]],Bb1=Sign[Re[qb]],Bb1=Re[bq]];
outB1[[i]]={b,Bb1,qb},{i,m}]
Do[If[Abs[Re[qb1[[i]]]]<λ2[[i]],Bb2=Sign[qb1[[i]][[4]]],Bb2=Re[bq1[[i]]]];outB2[[i]]={bb1[[i]],Bb2,qb1[[i]]},{i,m}]
outB=Catenate[{outB1,outB2}];
b2=outB[[All,1]];
B=outB[[All,2]];

CHSH Analysis of the Particle Data Received from Alice and Bob

nP1=0; nN1=0; nP2=0; nN2=0; nP3=0; nN3=0;nP4=0;  nN4=0;
Do[a1=a2[[i]];b1=b2[[i]];
aliceD=A[[i]];bobD=B[[i]];
If[(b1==45)&&(a1-b1==-45)&&aliceD*bobD==1,nP1++];
If[(b1==45)&&(a1-b1==-45)&&aliceD*bobD==-1,nN1++];If[(a1-b1)==-135&&aliceD*bobD==1,nP2++];
If[(a1-b1)==-135&&aliceD*bobD==-1,nN2++];If[(a1-b1)==45&&aliceD*bobD==1,nP3++];
If[(a1-b1)==45&&aliceD*bobD==-1,nN3++];
If[a1==90&&(a1-b1)==-45&&aliceD*bobD==1,nP4++];
If[a1==90&&(a1-b1)==-45&&aliceD*bobD==-1,nN4++],{i,2m}]
E1=N[(nP1-nN1)/(nP1+nN1)];
E2=N[(nP2-nN2)/(nP2+nN2)];
E3=N[(nP3-nN3)/(nP3+nN3)];
E4=N[(nP4-nN4)/(nP4+nN4)];
CHSH=N[Abs[E1-E2+E3+E4]];
Print["CHSH = ",CHSH]