****************************** fig7 ******************************* * * 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 fig7, 11, 12 * * 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,400,500,600) 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 400 go to (410,420,430,440,450,460,470) m+el+1 go to 900 410 Y=dsqrt(35/(64*pi))*dsin(theta)**3 & *cdexp(dcmplx(0,-3*phi)) return 420 Y=dsqrt(105/(32*pi))*dsin(theta)**2*dcos(theta) & *cdexp(dcmplx(0,-2*phi)) return 430 Y=dsqrt(21/(64*pi))*dsin(theta)*(5*dcos(theta)**2-1) & *cdexp(dcmplx(0,-phi)) return 440 Y=dsqrt(7/(16*pi))*(5*dcos(theta)**3-3*dcos(theta)) return 450 Y=-dsqrt(21/(64*pi))*dsin(theta)*(5*dcos(theta)**2-1) & *cdexp(dcmplx(0,phi)) return 460 Y=dsqrt(105/(32*pi))*dsin(theta)**2*dcos(theta) & *cdexp(dcmplx(0,2*phi)) return 470 Y=-dsqrt(35/(64*pi))*dsin(theta)**3 & *cdexp(dcmplx(0,3*phi)) return 500 go to (510,520,530,540,550,560,570,580,590) m+el+1 go to 900 510 Y=dsqrt(315/(512*pi))*dsin(theta)**4 & *cdexp(dcmplx(0,-4*phi)) return 520 Y=dsqrt(315/(64*pi))*dsin(theta)**3*dcos(theta) & *cdexp(dcmplx(0,-3*phi)) return 530 Y=dsqrt(45/(128*pi))*dsin(theta)**2*(7*dcos(theta)**2-1) & *cdexp(dcmplx(0,-2*phi)) return 540 Y=dsqrt(45/(64*pi))*dsin(theta)*(7*dcos(theta)**3-3*dcos(theta)) & *cdexp(dcmplx(0,-phi)) return 550 Y=dsqrt(9/(256*pi))*(35*dcos(theta)**4-30*dcos(theta)**2+3) return 560 Y=-dsqrt(45/(64*pi))*dsin(theta)*(7*dcos(theta)**3-3*dcos(theta)) & *cdexp(dcmplx(0,phi)) return 570 Y=dsqrt(45/(128*pi))*dsin(theta)**2*(7*dcos(theta)**2-1) & *cdexp(dcmplx(0,2*phi)) return 580 Y=-dsqrt(315/(64*pi))*dsin(theta)**3*dcos(theta) & *cdexp(dcmplx(0,3*phi)) return 590 Y=dsqrt(315/(512*pi))*dsin(theta)**4 & *cdexp(dcmplx(0,4*phi)) return 600 if (m.ne.0) go to 900 Y=dsqrt(11/(256*pi))*(63*dcos(theta)**5-70*dcos(theta)**3 & +15*dcos(theta)) return 900 Y=0 return end implicit real*8(a-h,o-z) real*8 L,Lmin,Lmax,Lp,Lpmin,Lpmax,kappa complex*16 Y,sum0,sum1,sum2,sum3,sum4,sum5,term,prefact, &tab2(64),tab3(64),tabx1(64),taby1(64),tabx2(64),taby2(64), &tabx3(64),taby3(64),tabx4(64),taby4(64),tabx(256),taby(256),rad integer el pi=4.d0*datan(1.d0) * parameters el=1 m=0 T=32*pi * nthetaG=64 nthetaG=8 * nalpha=64 nalpha=8 * nthetaf=32 nthetaf=8 nphif=1 * delL=0.001 delL=0.025 kappa=3 delthetaG=(pi/2)/nthetaG delthetaf=pi/nthetaf delphif=2*pi/nphif Lpmin=-1/dsqrt(T) Lpmax=1/dsqrt(T) open(unit=21,file='tab1.dat') open(unit=22,file='tab2.dat') open(unit=23,file='tab3.dat') open(unit=24,file='fig7a.dat') open(unit=25,file='fig7b.dat') open(unit=26,file='tab4a.dat') open(unit=27,file='tab4b.dat') open(unit=28,file='fig7c.dat') open(unit=29,file='fig7d.dat') * do the computations do 100 ithetaG=1,nthetaG thetaG=(ithetaG-0.5)*delthetaG if (ithetaG.eq.1) then write(21,51) thetaG 51 format(' {',F20.15,' ,') else if (ithetaG.ne.nthetaG) then write(21,52) thetaG 52 format(' ',F20.15,' ,') else write(21,53) thetaG 53 format(' ',F20.15,' }') end if 100 continue sum0=0.d0 do 905 itab=0,1 if (itab.eq.0) then Lmin=0.d0 Lmax=kappa*(el+0.5d0) else Lmin=-kappa*(el+0.5d0) Lmax=0.d0 end if nL=nint((Lmax-Lmin)/delL) do 900 ithetaG=1,nthetaG write(6,*) ithetaG,nthetaG thetaG=(ithetaG-0.5)*delthetaG thetafmin=thetaG thetafmax=pi-thetaG delthetaf=(thetafmax-thetafmin)/nthetaf sum1=0.d0 do 800 iL=1,nL L=Lmin+(iL-0.5)*delL nLp=max(ifloor(5*dsqrt(T)), & ifloor(3*dsqrt(T)*dabs(L))) delLp=(Lpmax-Lpmin)/nLp sum2=0.d0 do 700 iLp=1,nLp Lp=L+Lpmin+(iLp-0.5)*delLp gammat=-Lp*T nsav=n(gammat) gammasav=gamm(gammat) gamman=gammasav+2*pi*nsav sum3=0.d0 do 600 ithetaf=1,nthetaf thetaf=thetafmin+(ithetaf-0.5)*delthetaf sum4=0.d0 do 500 iphif=1,nphif phif=(iphif-0.5)*delphif sum5=0.d0 Phi0=datan(dsqrt(tan(thetaf)**2/ & ((dcos(thetaG)/dcos(thetaf))**2-1.d0)-1.d0)) Phi0p=pi-Phi0 theta0=theta(gammat,thetaf,phif,Phi0) theta0p=theta(gammat,thetaf,phif,Phi0p) phi00=phi(gammat,thetaf,phif,Phi0) phi00p=phi(gammat,thetaf,phif,Phi0p) alphamin=gamm(gammat) 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 if (ithetaf.ne.1.and.ithetaf.ne.nthetaf) then sum5=sum5*(dsin(thetaG)/(dcos(thetaf)**2* & dsqrt(dtan(thetaf)**2- & ((dcos(thetaG)/dcos(thetaf))**2-1.d0))* & dsqrt((dcos(thetaG)/dcos(thetaf))**2-1.d0)))* & delthetaf*dsin(thetaf) else sum5=sum5*dsqrt(2*delthetaf*dabs(dtan(thetaf))) end if sum5=sum5*dconjg(Y(el,m,thetaf,phif))* & (Y(el,m,theta0,phi00)+Y(el,m,theta0p,phi00p)) & *(-1)**nsav sum4=sum4+sum5 500 continue sum4=sum4*delphif sum3=sum3+sum4 600 continue sum3=sum3*dabs(dsin(gammat)) sum2=sum2+sum3 700 continue sum2=sum2*delLp sum1=sum1+sum2 800 continue sum1=sum1*delL*dsqrt(T)*prefact(el,T)/2 if (itab.eq.0) then tab2(ithetaG)=sum1 else tab3(ithetaG)=sum1 end if if (ithetaG.eq.1) then if (dimag(sum1).ge.0) then write(22+itab,851) dreal(sum1),dabs(dimag(sum1)) 851 format('{',f20.15,' +',f20.15,'*I',' ,') else write(22+itab,852) dreal(sum1),dabs(dimag(sum1)) 852 format('{',f20.15,' -',f20.15,'*I',' ,') end if else if (ithetaG.ne.nthetaG) then if (dimag(sum1).ge.0) then write(22+itab,853) dreal(sum1),dabs(dimag(sum1)) 853 format(' ',f20.15,' +',f20.15,'*I',' ,') else write(22+itab,854) dreal(sum1),dabs(dimag(sum1)) 854 format(' ',f20.15,' -',f20.15,'*I',' ,') end if else if (dimag(sum1).ge.0) then write(22+itab,855) dreal(sum1),dabs(dimag(sum1)) 855 format(' ',f20.15,' +',f20.15,'*I',' }') else write(22+itab,856) dreal(sum1),dabs(dimag(sum1)) 856 format(' ',f20.15,' -',f20.15,'*I',' }') end if end if write(24+itab,*) thetaG,dreal(sum1),dimag(sum1) sum0=sum0+sum1*dcos(thetaG) 900 continue 905 continue sum0=sum0*delthetaG do 910 ithetaG=1,nthetaG thetaG=(ithetaG-0.5)*delthetaG ang=thetaG rad=tab2(ithetaG) tabx1(ithetaG)=rad*dcos(ang) taby1(ithetaG)=rad*dsin(ang) tabx2(nthetaG-ithetaG+1)=-rad*dcos(ang) taby2(nthetaG-ithetaG+1)=rad*dsin(ang) rad=tab3(ithetaG) tabx3(ithetaG)=-rad*dcos(ang) taby3(ithetaG)=-rad*dsin(ang) tabx4(nthetaG-ithetaG+1)=rad*dcos(ang) taby4(nthetaG-ithetaG+1)=-rad*dsin(ang) 910 continue do 920 ithetaG=1,nthetaG tabx(ithetaG)=tabx1(ithetaG) tabx(ithetaG+nthetaG)=tabx2(ithetaG) tabx(ithetaG+2*nthetaG)=tabx3(ithetaG) tabx(ithetaG+3*nthetaG)=tabx4(ithetaG) taby(ithetaG)=taby1(ithetaG) taby(ithetaG+nthetaG)=taby2(ithetaG) taby(ithetaG+2*nthetaG)=taby3(ithetaG) taby(ithetaG+3*nthetaG)=taby4(ithetaG) 920 continue do 960 i=1,4*nthetaG+1 if (i.le.4*nthetaG) then write(28,950) dreal(tabx(i)),dreal(taby(i)) 950 format(f20.15,' ',f20.15) else write(28,950) dreal(tabx(1)),dreal(taby(1)) end if if (i.eq.1) then if (dreal(taby(i)).ge.0) then write(26,951) dreal(tabx(i)),dabs(dreal(taby(i))) 951 format('{',f20.15,' +',f20.15,'*I',' ,') else write(26,952) dreal(tabx(i)),dabs(dreal(taby(i))) 952 format('{',f20.15,' -',f20.15,'*I',' ,') end if else if (i.lt.4*nthetaG) then if (dreal(taby(i)).ge.0) then write(26,951) dreal(tabx(i)),dabs(dreal(taby(i))) else write(26,952) dreal(tabx(i)),dabs(dreal(taby(i))) end if else if (i.eq.4*nthetaG) then if (dreal(taby(i)).ge.0) then write(26,951) dreal(tabx(i)),dabs(dreal(taby(i))) else write(26,952) dreal(tabx(i)),dabs(dreal(taby(i))) end if else if (dreal(taby(1)).ge.0) then write(26,951) dreal(tabx(1)),dabs(dreal(taby(1))) else write(26,952) dreal(tabx(1)),dabs(dreal(taby(1))) end if end if 960 continue do 970 i=1,4*nthetaG+1 if (i.le.4*nthetaG) then write(29,950) dimag(tabx(i)),dimag(taby(i)) else write(29,950) dimag(tabx(1)),dimag(taby(1)) end if if (i.le.4*nthetaG) then write(27,950) dimag(tabx(i)),dimag(taby(i)) else write(27,950) dimag(tabx(1)),dimag(taby(1)) end if if (i.eq.1) then if (dimag(taby(i)).ge.0) then write(27,951) dimag(tabx(i)),dabs(dimag(taby(i))) else write(27,952) dimag(tabx(i)),dabs(dimag(taby(i))) end if else if (i.lt.4*nthetaG) then if (dimag(taby(i)).ge.0) then write(27,951) dimag(tabx(i)),dabs(dimag(taby(i))) else write(27,952) dimag(tabx(i)),dabs(dimag(taby(i))) end if else if (i.eq.4*nthetaG) then if (dimag(taby(i)).ge.0) then write(27,951) dimag(tabx(i)),dabs(dimag(taby(i))) else write(27,952) dimag(tabx(i)),dabs(dimag(taby(i))) end if else if (dimag(taby(1)).ge.0) then write(27,951) dimag(tabx(1)),dabs(dimag(taby(1))) else write(27,952) dimag(tabx(1)),dabs(dimag(taby(1))) end if end if 970 continue write(6,980) dreal(sum0),dimag(sum0) 980 format('integrated total = ',f9.5,2x,f9.5) write(6,*) 'press ENTER to exit' * read(5,*) stop end