Validation of the Local QM Product Calculation Prediction Using Pauli Matrices
and Quaternions with 3D vectors, Based on Joy Christian 3-Sphere Model.
Created by Fred Diether Mar. 2022

<<"clifford.m"
Qcoordinates={1,i,j,k};
m=20000;
s1=ConstantArray[0,m];
s2=ConstantArray[0,m];
σs1=ConstantArray[0,m];
σs2=ConstantArray[0,m];
a1=ConstantArray[0,m];
b1=ConstantArray[0,m];
ra1=ConstantArray[0,m];
rb1=ConstantArray[0,m];
qA=ConstantArray[0,m];
qB=ConstantArray[0,m];
A=ConstantArray[0,m];
B=ConstantArray[0,m];
pc=ConstantArray[0,m];
plotAB=Table[{0,0},m];
I3=Pseudoscalar[3];

Do[s=RandomPoint[Sphere[]];  (*Uniform Unit 3D Vectors; Hidden Variable*)
s1[[h]]=s;
s2[[h]]=-s;
σs1[[h]]=PauliMatrix[1]*s[[1]]+PauliMatrix[2]*s[[2]]+PauliMatrix[3]*s[[3]];
σs2[[h]]=-(PauliMatrix[1]*s[[1]]+PauliMatrix[2]*s[[2]]+PauliMatrix[3]*s[[3]]),{h,m}]

Do[a=RandomPoint[Sphere[]];  (*Uniform Unit 3D Vectors*)
a1[[h]]=a;
σa=PauliMatrix[1]*a[[1]]+PauliMatrix[2]*a[[2]]+PauliMatrix[3]*a[[3]];
cosas1=Re[FullSimplify[Extract[Flatten[(1/2)((({{1, 0}})).σa.σs1[[h]].(({{1}, {0}}))+(({{0, 1}})).σa.σs1[[h]].(({{0}, {1}})))],1]+0.00001i]];
(*Particle - Detector interaction*)
ra=Cross[a,s1[[h]]];
ra1[[h]]=ra;
qA[[h]]={cosas1,ra[[1]],ra[[2]],ra[[3]]}.Qcoordinates;
A[[h]]=Sign[cosas1],{h,m}]

Do[b=RandomPoint[Sphere[]];  (*Uniform Unit 3D Vectors*)
b1[[h]]=b;
σb=PauliMatrix[1]*b[[1]]+PauliMatrix[2]*b[[2]]+PauliMatrix[3]*b[[3]];
cosbs2=Re[FullSimplify[Extract[Flatten[(1/2)((({{1, 0}})).σs2[[h]].σb.(({{1}, {0}}))+(({{0, 1}})).σs2[[h]].σb.(({{0}, {1}})))],1]+0.00001i]];
(*Particle - Detector interaction*)
rb=Cross[s2[[h]],b];
rb1[[h]]=rb;
qB[[h]]={cosbs2,rb[[1]],rb[[2]],rb[[3]]}.Qcoordinates;
B[[h]]=Sign[cosbs2],{h,m}]

Do[r0=Expand[({e[1], e[2],e[3]}).(Re[qA[[h]]]*Limit[Cross[s4,b1[[h]]]+0.000001,s4->Sign[Re[qB[[h]]]]b1[[h]]]+Re[qB[[h]]]*Limit[Cross[a1[[h]],s3]+0.000001,s3->Sign[Re[qA[[h]]]]a1[[h]]]-
Cross[Limit[Cross[a1[[h]],s3]+0.000001,s3->Sign[Re[qA[[h]]]]a1[[h]]],
Limit[Cross[s4,b1[[h]]]+0.000001,s4->Sign[Re[qB[[h]]]]b1[[h]]]])/
(Sin[ArcCos[a1[[h]].b1[[h]]]])];
Lr0=InnerProduct[I3,r0];
qpc=(Re[qA[[h]]]*Re[qB[[h]]]-ra1[[h]].rb1[[h]])+Lr0;  (*Product Calculation*)
pc[[h]]=qpc;
ϕa=ArcTan[a1[[h]][[1]],a1[[h]][[2]]];
ϕb=ArcTan[b1[[h]][[2]],b1[[h]][[1]]];
If[ϕa*ϕb>0,angle=ArcCos[a1[[h]].b1[[h]]]/Degree, 
angle=(2π-ArcCos[a1[[h]].b1[[h]]])/Degree];
plotAB[[h]]={angle,Re[qpc]},{h,m}]

simulation=ListPlot[plotAB,PlotMarkers->{Automatic,Small},AspectRatio->9/16,Ticks->{{{90,90°},{180,180°},{0,0°},{270,270°},{360,360°}},Automatic},
GridLines->Automatic,AxesOrigin->{0,-1.0}];
negcos=Plot[-Cos[x Degree],{x,0,360},PlotStyle->{Magenta}];
p1=Plot[-1+2x1 Degree/π,{x1,0,180},PlotStyle->{Gray,Dashed}];
p2=Plot[3-2x2 Degree/π,{x2,180,360},PlotStyle->{Gray,Dashed}];
Show[ simulation,p1,p2, negcos]

Blue is the  correlation data, magenta is the negative cosine curve for an exact match.

AveA=N[Total[A]/m];
AveB=N[Total[B]/m];
Print[" <A> = ",AveA,"  <B> = ",AveB];
meanpc=Expand[Mean[pc]];
Print["Cross products vanish, meanpc = ",meanpc];