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