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

fig 6, compute P(1) distribution as a function of characteristic angular momentum Lc

in variable names below, an appended ‘p’ means prime, so e.g.  Lpmin and Lpmax are min and max values of Lp = L’ (L-prime);  an appended ‘t’ means tilde, so e.g. γt = γ-tilde (γ with tilde mark over it)

coarse grid is shown below, with nΦ0=8; nα=8; nθf=8; nϕf=1; delL=0.025; requires about 15 minutes to run using Mathematica

fine grid is nΦ0=64; nα=64; nθf=32; nϕf=1; delL=0.001; requires about 1 day to run using FORTRAN code, with that code running about 100x faster than Mathematica code having the same parameters  (although part of that difference arises from the faster processor employed for the FORTRAN run)

Clear["Global`*"]; 
el=1; m=0; T=N[32 π]; nΦ0=8;nα=8; nθf=8; nϕf=1; delL=0.025; 
delΦ0=N[π/nΦ0];  delθf=N[π/nθf];  delϕf=N[2π/nϕf];  Lpmin=-N[1/Sqrt[T]]; Lpmax=N[1/Sqrt[T]];
Lmin=-4.5; Lmax=4.5;   nL=Round[(Lmax-Lmin)/delL];
θ[γ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_]:=(Sqrt[T]/2)N[(1/(√2π I ))(((1/(2π I T)))^(1/2))Exp[I (((el+1/2))^2) T/2]];
tab1=Table[L=Lmin+(iL-0.5)delL; L,{iL,nL}];
tab2=prefact[T]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]]delθf Total[Table[ θf=(iθf-0.5) delθf; Sin[θf] delϕf Total[Table[ ϕf=(iϕf-0.5) delϕf;
delΦ0 Total[Table[Φ0=(iΦ0-0.5)delΦ0;
Conjugate[SphericalHarmonicY[el,m,θf,ϕf]]SphericalHarmonicY[el,m,θ[γt,θf,ϕf,Φ0],ϕ[γt,θf,ϕf,Φ0]](((-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Φ0,nΦ0}]],{iϕf,nϕf}]],{iθf,nθf}]],{iLp,nLp}]],{iL,nL}];
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];
Show[{ListPlot[Transpose[{tab1,Re[tab2]}],PlotRange->{{-4.5,4.5},{-3,5}},PlotStyle->{ColorData[97,"ColorList"][[1]]}],ListPlot[Transpose[{tab1,Im[tab2]}],PlotStyle->{ColorData[97,"ColorList"][[2]]}]},ImageSize->Medium]