****************************** 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
