computer code from R. M. Feenstra, “Path distributions for describing eigenstates of orbital angular momentum”, arXiv:2308.02884 [quant-ph].

figs 7, 11, 12 - compute P(2) angular distribution

same comments here as for preceding fig 6 code, except replace nΦ0 in the comments there with nθg as applicable here

Clear["Global`*"]; 
el=1; m=0; T=N[32 π]; nθg=8; nα=8; nθf=8; nϕf=1; delL=0.025;
delθg=N[(π/2)/nθg];   delϕf=N[2π/nϕf];  Lpmin=-N[1/Sqrt[T]]; Lpmax=N[1/Sqrt[T]];
θ[γt_,θf_,ϕf_,Φ0_]:=N[(stf=Sin[θf]; ctf=Cos[θf]; ArcCos[-stf Cos[Φ0]Sin[γt]+ctf Cos[γt]])]; 
ϕ[γt_,θf_,ϕf_,Φ0_]:=N[(stf=Sin[θf]; ctf=Cos[θf]; spf=Sin[ϕf]; cpf=Cos[ϕf]; ArcTan[-spf Sin[Φ0] Sin[γt]+cpf(ctf Cos[Φ0] Sin[γt]+stf Cos[γt]),cpf  Sin[Φ0] Sin[γt]+spf(ctf  Cos[Φ0] Sin[γt]+stf Cos[γt])])];
n[γt_]:=N[Floor[γt/(2 π)]-(1+2Floor[γt/(2 π)])Mod[Floor[γt/π],2]];
γ[γt_]:=N[ Abs[Mod[γt-π,2π]-π]];prefact[T_]:=N[(1/(√2π I ))(((1/(2π I T)))^(1/2))Exp[I (((el+1/2))^2) T/2]];
tab[Lmin_,Lmax_,nL_]:=Table[θg=(iθg-0.5)delθg; θfmin=θg;  θfmax=N[π]-θg; delθf=N[(θfmax-θfmin)/nθf]; (Sqrt[T] prefact[T]  delL/2)  Total[Table[L=Lmin+(iL-0.5)delL;  nLp=Max[Floor[5 Sqrt[T]],Floor[3 Sqrt[T] Abs[L]]];  delLp=N[(Lpmax-Lpmin)/nLp];
delLp Total[Table[Lp=L+Lpmin+(iLp-0.5)delLp;  γt=-Lp T; Abs[Sin[γt]]
Total[Table[ θf=θfmin+(iθf-0.5) delθf; 
delϕf Total[Table[ ϕf=(iϕf-0.5) delϕf;
 Φ0=ArcTan[((((((Tan[θf]))^2)/((((Cos[θg]/Cos[θf]))^2)-1) -1))^(1/2))];  Φ0p=π-Φ0;
 If[iθf!=1&&iθf!=nθf,delθf Sin[θf](Sin[θg]/((((Cos[θf]))^2) ((((((Tan[θf]))^2)-((((Cos[θg]/Cos[θf]))^2)-1) ))^(1/2))((((((Cos[θg]/Cos[θf]))^2)-1))^(1/2)))),(((2 delθf Abs[Tan[θf]]))^(1/2))]
Conjugate[SphericalHarmonicY[el,m,θf,ϕf]](SphericalHarmonicY[el,m,θ[γt,θf,ϕf,Φ0],ϕ[γt,θf,ϕf,Φ0]]+SphericalHarmonicY[el,m,θ[γt,θf,ϕf,Φ0p],ϕ[γt,θf,ϕf,Φ0p]])(((-1))^(n[γt])) (αmin=γ[γt]; αmax=π; delα=(αmax-αmin)/nα; 
Total[Table[α=αmin+(iα-0.5)delα; αn=(α+2 N[π] n[γt]); 
If[iα!=1,delα αn Exp[I ((αn)^2)/(2 T)]/(((Cos[γ[γt]]-Cos[α]))^(1/2)),α1n=γ[γt]+delα/4+2 N[π]n[γt];α1n Exp[I ((α1n)^2)/(2 T)]2 ((√delα)/(√Sin[γ[γt]+delα/8]))],{iα,nα}]]),{iϕf,nϕf}]],{iθf,nθf}]],{iLp,nLp}]],{iL,nL}]],{iθg,nθg}];
tab1=Table[θg=(iθg-0.5)delθg; θg,{iθg,nθg}];
Lmin=0; Lmax=3(el+1/2);   nL=Round[(Lmax-Lmin)/delL];
tab2=tab[Lmin,Lmax,nL];
tabx1=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab2[[ iθg]];rad Cos[ang],{iθg,nθg}];
taby1=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab2[[ iθg]];rad Sin[ang],{iθg,nθg}];
tabx2=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab2[[ iθg]];-rad Cos[ang],{iθg,nθg,1,-1}];
taby2=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab2[[ iθg]];rad Sin[ang],{iθg,nθg,1,-1}];
Lmin=-3(el+1/2); Lmax=0;   nL=Round[(Lmax-Lmin)/delL];
tab3=tab[Lmin,Lmax,nL];
tabx3=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab3[[ iθg]];-rad Cos[ang],{iθg,nθg}];
taby3=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab3[[ iθg]];-rad Sin[ang],{iθg,nθg}];
tabx4=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab3[[ iθg]];rad Cos[ang],{iθg,nθg,1,-1}];
taby4=Table[θg=(iθg-0.5)delθg;ang=θg;rad=tab3[[ iθg]];-rad Sin[ang],{iθg,nθg,1,-1}];
tabx=Join[tabx1,tabx2,tabx3,tabx4,{tabx1[[1]]}];taby=Join[taby1,taby2,taby3,taby4,{taby1[[1]]}];
fname="tab1.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab1];  Close[stmp];
fname="tab2.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab2];  Close[stmp];
fname="tab3.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab3];  Close[stmp];
ListLinePlot[Transpose[{tabx,taby}],PlotRange->{{-1,1},{-1,1}},AspectRatio->1]