implicit real*8(a-h,o-z)
      parameter(Lbiggest=120000,MaxPts=5000,ActualPts=300)
      dimension qint(MaxPts),qintp(MaxPts),x(MaxPts),Delta(0:Lbiggest)
      dimension indexSurvey(10),rSurvey(10),StepSize(10)
      common/KmuL/wavenum,ReducedMass,L
      data Pi/3.14159265358979d0/
      external yfcn
            character*9 today
            character*8 now
            call date (today)
            call time (now)   

c     write(6,*) '%Screened Coulomb Potential only', today, ' ',now
      write(6,*) '%Azzi at large r, screened at small r', today, ' ',now
c     write(6,*) '%Azzi and Slaman potential', today, ' ',now
c     write(6,*) '%Aziz extrapolated to all r', today, ' ',now

c     write(40.*) '%Screened Coulomb Potential only', today, ' ',now
      write(40,*) '%Azzi -large r, screened -small r',today, ' ',now
c     write(40,*) '%Azzi and Slaman potential', today, ' ',now
c     write(40,*) '%Aziz extrapolated to all r', today, ' ',now

      write(6,*) 'eMinEV,eMaxEV,nE'
      read(5,*) eMinEV,eMaxEV,nE
      write(6,*)'%energyEV, Qtot, Qv, Delta(Lmax), Lmax'
      write(40,*)'%energyEV, Qtot, Qv, Delta(Lmax), Lmax'

      nq = ActualPts
      eStartAU=eMinEV*3.674872d-2
      eFinishAU=eMaxEV*3.674872d-2
      xStart=dsqrt(eStartAU)
      xFinish=dsqrt(eFinishAU)
      xStep=0.d0
c     IF(nE>1) THEN
      IF(nE .GT .1) THEN
       xStep=(xFinish-xStart)/(nE-1.d0)
      ENDIF
      DO iEnergy=1,nE
        xE=xStart+(iEnergy-1)*xStep
        energyAU=xE*xE
        energyEV=energyAU/3.674872d-2

c      Lmax=(10.d0**2.85d0)*energyEV**0.41d0
       Lmax=(10.d0**3)*energyEV**0.5d0
       Lmax=2*(Lmax/2)

      write(6,*) 'Lmax,energyEV=',Lmax,energyEV

      DO iL=0,Lmax,2
       L=iL

      ReducedMass= 32782.41335975446
      C6=64.3d0 
c *** 35.967546282d0 AMU for this Ar isotope!
      wavenum=dsqrt(2.d0*energyAU*ReducedMass)
     
      Rmin=1.d-06
      Qmin=dsqrt(Rmin)
      Rmax=((2.d0*ReducedMass*C6/wavenum)**0.2d0) *1.5d0
      IF(Rmax.lt.8.d0*ZeroRoot) Rmax=8.d0*ZeroRoot
      Qmax=dsqrt(Rmax)
      Rlast=0.d0
      Nsurvey=1501
      hStep=(Qmax-Qmin)/(Nsurvey-1.d0)
      XkinLast=-1.d0
      iCount = 0

      DO iS=1,Nsurvey
      q = Qmin+(iS-1)*hStep
      r=q*q
      XkinNew=yfcn(r)
      IF(XkinNew*XkinLast .LT. 0.d0) THEN
          iCount = iCount+1
          indexSurvey(iCount)=iS
          rSurvey(iCount)=(r+Rlast)/2.d0
          StepSize(iCount)=r-Rlast
c         write(6,*) 'iCount ',iCount,rSurvey(iCount),StepSize(iCount)
      ENDIF
      rLast=r
      XkinLast=XkinNew
      ENDDO

      guessRoot=(L+0.5d0)/wavenum
      ZeroRoot=guessRoot
c     write(6,*) 'ZeroRoot,L,E=',ZeroRoot,L,energyAU

      h1=-0.5d0
      iroute=1
      error=1.d-12
      h1=StepSize(1)/3.d0
      call root(rSurvey(1),h1,iroute,error,GoodRoot,ErrorKin,yfcn)
c     write(6,*) 'GoodRoot,ErrorKin=',GoodRoot,ErrorKin
     
      IF(ZeroRoot-GoodRoot.GT.0.d0) THEN
      q0=0.d0
      qf= dsqrt(ZeroRoot-GoodRoot)
      hq=(qf-q0)/(nq-1.d0)
     
      do iq=1,nq
       q=q0+(iq-1)*hq
       r=GoodRoot + q*q
       qint(iq) = dsqrt(dabs(yfcn(r)))  *q*2.d0
      enddo

      qint(1)=0.d0

      ELSE
      q0=0.d0
      qf= dsqrt(GoodRoot-ZeroRoot)
      hq=(qf-q0)/(nq-1.d0)
     
      do iq=1,nq
       q=q0+(iq-1)*hq
       r=ZeroRoot + q*q
      qint(iq)=-dsqrt(dabs(wavenum*wavenum-(L+0.5d0)**2/(r*r)))*q*2.d0
      enddo

      qint(1)=0.d0
      ENDIF


c***  rint(f,na,nb,nq,h)
      FirstIntegral=rint(qint,1,nq,10,hq)

      Rmax=((2.d0*ReducedMass*C6/wavenum)**0.2d0) *1.5d0
      IF(Rmax.lt.8.d0*ZeroRoot) Rmax=8.d0*ZeroRoot

      IF(ZeroRoot-GoodRoot.GT.0.d0) THEN
      q0=0.d0
      qf= dsqrt(Rmax-ZeroRoot)
      hq2=(qf-q0)/(nq-1.d0)
      do iq=1,nq
       q=q0+(iq-1)*hq2
       r=ZeroRoot + q*q
      qint(iq) = dsqrt(dabs(yfcn(r))) 
      qint(iq)=qint(iq)-dsqrt(dabs(wavenum*wavenum-(L+0.5d0)**2/(r*r)))
      qint(iq) = qint(iq)   *q*2.d0

      enddo

      ELSE
      q0=0.d0
      qf= dsqrt(Rmax-GoodRoot)
      hq2=(qf-q0)/(nq-1.d0)
      do iq=1,nq
       q=q0+(iq-1)*hq2
       r=GoodRoot + q*q
      qint(iq) = dsqrt(dabs(yfcn(r))) 
      qint(iq)=qint(iq)-dsqrt(dabs(wavenum*wavenum-(L+0.5d0)**2/(r*r)))
      qint(iq) = qint(iq)   *q*2.d0

      enddo
      ENDIF

c***  rint(f,na,nb,nq,h)
      SecondIntegral=rint(qint,1,nq,10,hq2)

      phaseshift=FirstIntegral+SecondIntegral
      Delta(L)=phaseshift
     
      write(40,*) L,energyEV,delta(L)/Pi

       ENDDO

      call SigmaCalc(energyEV,Lmax,Delta,SigmaTOT,SigmaV)
      write(6,156) energyEV,SigmaTOT,SigmaV,Delta(Lmax),Lmax
      write(140,156) energyEV,SigmaTOT,SigmaV,Delta(Lmax),Lmax
156      format(4(1x,1pd12.5),(1x,I7))

      ENDDO

      stop
      end

c***
      Double Precision Function yfcn(r)
      implicit real*8(a-h,o-z)
      common/KmuL/wavenum,ReducedMass,L
 
      yfcn = wavenum**2-(L+0.5d0)**2/(r*r) -2.d0*ReducedMass*Vpot(r)
c***  write(6,*) 'r,Vpot,yfcn=',r,yfcn,Vpot(r)
      return
      end

c***
      Double Precision Function Vpot(r)
      implicit real*8(a-h,o-z)
      double precision eps,D,c6,c8,c10,Astar,alpha,beta
      double precision rmin,Fscale,vscale,vlr
      double precision a1,a2,a3,a4,a5,alpha2,series1
      double precision za,zb,r1,r2,b1,b2,b3,series2,beta2
      beta = 1.72832d0

c convert from angstroms to a.u.
      rmin = 3.761d0/0.529177249d0

c convert from Kelvin to a.u.
      eps = 143.25d0*3.166829d-6

      astar = 1.13211845d5
      alpha = 9.00053441d0
      beta = -2.60270226d0
      c6 = 1.09971113d0
      c8 = 0.54511632d0
      c10 = 0.39278653d0
      D = 1.0400d0
      a1 = 9.6379d0
      a2 = 32.7048d0
      a3 = 827.11d0
      a4 = -3310.0d0
      a5 = 10417d0
      alpha2 = 14.880d0
      za = 18.0d0
      zb = 18.0d0
      r1 = 0.40d0
      r2 = 3.20d0
      b1 = -1.78799d0
      b2 = 1.03663d0
      b3 = -0.14975d0
      beta2 = 1.72832d0
      rp = 4.25d0
      ascr = 0.59516d0
      z1 = 18.d0
      z2 = 18.d0

      xscale = r/rmin

      vscale = Astar*dexp(-alpha*xscale + beta*xscale*xscale)

      vlr = c6/(xscale**6) + c8/(xscale**8) + c10/(xscale**10)

c     series1 corrected 5/6/99 for error pointed out by c. Greene.

      series1 = 1d0+a1*r+a2*r**2+a3*r**3+a4*r**4+a5*r**5
 
      series2 = 1d0+b1*r+b2*r**2+b3*r**3

      if (xscale .gt. D)then
       Fscale = 1.0d0
      else
       Fscale = dexp(-(D/xscale-1.0d0)*(D/xscale-1.0d0))
      endif 

c     the following is the potential at large radii from Aziz and 
c     Slaman (1990) and used by Greene in his initial calculations.
c     Note that a lot of the constants defined above are not used here.

c       Vpot = eps*(vscale - Fscale*vlr)

c     the following is the potential at all radii from Aziz and 
c     Slaman (1990) using an extension to small radii from 
c     Pathak and Thakkar (1987)

c      if (xscale .ge .(r2/rmin)) then
c        Vpot = eps*(vscale - Fscale*vlr)
c      else
c       if (xscale .le .(r1/rmin))then
c        Vpot = za*zb/xscale/rmin*series1*dexp(-alpha2*xscale*rmin)
c       else 
c        Vpot = za*zb/(xscale*rmin)*series2*dexp(-beta2*xscale*rmin)
c       endif
c      endif

c     the following is the potential at large radii from Aziz and 
c     Slaman (1990) with an extension to small radii using a 
c     screened Coulomb formula fitted at r = rp

       if (r .ge .rp) then
         Vpot = eps*(vscale - Fscale*vlr)
       else
         Vpot = z1*z2/r*dexp(-r/ascr)
       endif

c     The following is a pure screened Coulomb potential

c       Vpot = z1*z2/r*dexp(-r/ascr)



      return
      end


      subroutine SigmaCalc(EnergyEV,Lmax,Defect,SigmaTOT,SigmaV)
      implicit real*8(a-h,o-z)
      dimension Defect(0:Lmax)
      data pi/3.14159265358979d0/

c     Defect is the array containing the scattering phaseshifts for each L.
c
      energy=EnergyEV
      reducedMass = 32782.41336d0
      energyAU=energy/27.2113961d0
      wavenum=dsqrt(2.d0*reducedMass*energyAU)
      sum0=0.d0
      sum2=0.d0
      altsum=0.d0
      altprod=1.d0

c -- the factor 'conversion' takes us to m^2 units!

      conversion = (5.29178d-9*5.29178d-9)/100.d0**2

c -- Note: change to 4.d0 instead of 2.d0 for non-identical particles

      altfront=(4.d0*pi/wavenum**2)*conversion
      total=0.d0

c -- This version of the code is for scattering of identical bosons

      do ll=0,Lmax-2,2
     
      altprod=(ll+1.d0)*(ll+2.d0)/(ll+1.5d0)
      altprod=altprod*dsin(Defect(ll+2)-Defect(ll))**2
      altsum=altsum+altprod*altfront

      factor=altfront*2*dsin(Defect(ll))**2 *(2.d0*ll+1)
      total=total+factor

      enddo
      SigmaTOT=total
      SigmaV=altsum
      RETURN
      end