****************************** fig6 ******************************* * * computer code from R. M. Feenstra, Path distributions for describing * eigenstates of orbital angular momentum, arXiv:2308.02884 [quant-ph] * * FORTRAN computations to match Mathematica code for fig6 * * variable names here match those in the Mathematica code, with the * exception that names which start with an upper-case letter (such * as Theta0) are preceded here by 'cap' (capTheta0) since FORTRAN is * not case sensitive real*8 function theta(gammat,thetaf,phif,Phi0) implicit real*8(a-h,o-z) stf=dsin(thetaf) ctf=dcos(thetaf) theta=dacos(-stf*dcos(Phi0)*dsin(gammat)+ctf*dcos(gammat)) return end real*8 function phi(gammat,thetaf,phif,Phi0) implicit real*8(a-h,o-z) stf=dsin(thetaf) ctf=dcos(thetaf) spf=dsin(phif) cpf=dcos(phif) phi=datan2(cpf*dsin(Phi0)*dsin(gammat)+ & spf*(ctf*dcos(Phi0)*dsin(gammat)+stf*dcos(gammat)), & -spf*dsin(Phi0)*dsin(gammat)+ & cpf*(ctf*dcos(Phi0)*dsin(gammat)+stf*dcos(gammat))) return end complex*16 function prefact(el,T) implicit real*8(a-h,o-z) integer el pi=4.d0*datan(1.d0) prefact=(1/(dsqrt(2.d0)*pi*dcmplx(0,1)))* & (1/(dsqrt(2.d0*pi*T)*dcmplx(1/dsqrt(2.d0),1/dsqrt(2.d0))))* & cdexp(dcmplx(0,(el+0.5)**2*T/2)) return end complex*16 function Y(el,m,theta,phi) implicit real*8(a-h,o-z) integer el pi=4.d0*datan(1.d0) go to (100,200,300) el+1 go to 900 100 Y=1/dsqrt(4*pi) return 200 go to (210,220,230) m+el+1 go to 900 210 Y=dsqrt(3/(8*pi))*dsin(theta)*cdexp(dcmplx(0,-phi)) return 220 Y=dsqrt(3/(4*pi))*dcos(theta) return 230 Y=-dsqrt(3/(8*pi))*dsin(theta)*cdexp(dcmplx(0,phi)) return 300 go to (310,320,330,340,350) m+el+1 go to 900 310 Y=dsqrt(15/(32*pi))*dsin(theta)**2*cdexp(dcmplx(0,-2*phi)) return 320 Y=dsqrt(15/(8*pi))*dsin(theta)*dcos(theta)*cdexp(dcmplx(0,-phi)) return 330 Y=dsqrt(5/(16*pi))*(3*dcos(theta)**2-1) return 340 Y=-dsqrt(15/(8*pi))*dsin(theta)*dcos(theta)*cdexp(dcmplx(0,phi)) return 350 Y=dsqrt(15/(32*pi))*dsin(theta)**2*cdexp(dcmplx(0,2*phi)) return 900 Y=0 return end implicit real*8(a-h,o-z) real*8 L,Lmin,Lmax,Lp,Lpmin,Lpmax complex*16 Y,sum0,sum1,sum2,sum3,sum4,sum5,term,prefact integer el pi=4.d0*datan(1.d0) * parameters el=1 m=0 T=32*pi * nPhi0=64 nPhi0=8 * nalpha=64 nalpha=8 * nthetaf=32 nthetaf=8 nphif=1 * delL=0.001 delL=0.025 delPhi0=pi/nPhi0 delthetaf=pi/nthetaf delphif=2*pi/nphif Lpmin=-1/dsqrt(T) Lpmax=1/dsqrt(T) Lmin=-4.5 Lmax=4.5 nL=nint((Lmax-Lmin)/delL) open(unit=21,file='tab1.dat') open(unit=22,file='tab2.dat') open(unit=23,file='fig6.dat') * do the computations do 100 iL=1,nL L=Lmin+(iL-0.5)*delL if (iL.eq.1) then write(21,51) L 51 format(' {',F20.15,' ,') else if (iL.ne.nL) then write(21,52) L 52 format(' ',F20.15,' ,') else write(21,53) L 53 format(' ',F20.15,' }') end if 100 continue sum0=0.d0 do 900 iL=1,nL write(6,*) iL,nL L=Lmin+(iL-0.5)*delL nLp=max(ifloor(5*dsqrt(T)), & ifloor(3*dsqrt(T)*dabs(L))) delLp=(Lpmax-Lpmin)/nLp sum1=0.d0 do 800 iLp=1,nLp Lp=L+Lpmin+(iLp-0.5)*delLp gammat=-Lp*T nsav=n(gammat) gammasav=gamm(gammat) gamman=gammasav+2*pi*nsav sum2=0.d0 do 700 ithetaf=1,nthetaf thetaf=(ithetaf-0.5)*delthetaf sum3=0.d0 do 600 iphif=1,nphif phif=(iphif-0.5)*delphif sum4=0.d0 do 500 iPhi0=1,nPhi0 Phi0=(iPhi0-0.5)*delPhi0 sum5=0.d0 theta0=theta(gammat,thetaf,phif,Phi0) phi00=phi(gammat,thetaf,phif,Phi0) alphamin=gammasav alphamax=pi delalpha=(alphamax-alphamin)/nalpha do 400 ialpha=1,nalpha alpha=alphamin+(ialpha-0.5)*delalpha alphan=alpha+2*pi*nsav if (ialpha.ne.1) then term=delalpha*alphan* & cdexp(dcmplx(0.d0,alphan**2/(2*T)))/ & dsqrt(dcos(gamm(gammat))-dcos(alpha)) else alpha1n=gammasav+delalpha/4+2*pi*nsav term=alpha1n*2*dsqrt(delalpha)* & cdexp(dcmplx(0.,alpha1n**2/(2*T)))/ & dsqrt(dsin(gammasav+delalpha/8)) end if sum5=sum5+term 400 continue sum5=sum5*dconjg(Y(el,m,thetaf,phif))* & Y(el,m,theta0,phi00)*(-1)**nsav sum4=sum4+sum5 500 continue sum4=sum4*delPhi0 sum3=sum3+sum4 600 continue sum3=sum3*dsin(thetaf)*delphif sum2=sum2+sum3 700 continue sum2=sum2*dabs(dsin(gammat))*delthetaf sum1=sum1+sum2 800 continue sum1=sum1*delLp*dsqrt(T)*prefact(el,T)/2 if (iL.eq.1) then if (dimag(sum1).ge.0) then write(22,851) dreal(sum1),dabs(dimag(sum1)) 851 format('{',f20.15,' +',f20.15,'*I',' ,') else write(22,852) dreal(sum1),dabs(dimag(sum1)) 852 format('{',f20.15,' -',f20.15,'*I',' ,') end if else if (iL.ne.nL) then if (dimag(sum1).ge.0) then write(22,853) dreal(sum1),dabs(dimag(sum1)) 853 format(' ',f20.15,' +',f20.15,'*I',' ,') else write(22,854) dreal(sum1),dabs(dimag(sum1)) 854 format(' ',f20.15,' -',f20.15,'*I',' ,') end if else if (dimag(sum1).ge.0) then write(22,855) dreal(sum1),dabs(dimag(sum1)) 855 format(' ',f20.15,' +',f20.15,'*I',' }') else write(22,856) dreal(sum1),dabs(dimag(sum1)) 856 format(' ',f20.15,' -',f20.15,'*I',' }') end if end if write(23,*) Lmin+(iL-0.5)*delL,dreal(sum1),dimag(sum1) sum0=sum0+sum1 900 continue sum0=sum0*delL write(6,910) dreal(sum0),dimag(sum0) 910 format('integrated total = ',f9.5,2x,f9.5) write(6,*) 'press ENTER to exit' * read(5,*) stop end