MODULE radialsch
  USE GlobalMath
  USE atomdata

  IMPLICIT NONE

CONTAINS

  SUBROUTINE boundsch(Grid,Pot,Orbit,l,start,nroot,emin,ierr)
    !  pgm to solve radial schroedinger equation for nroot bound state
    !    energies and wavefunctions for angular momentum l
    !    with potential rv/r, given in uniform mesh of n points
    !   r=i*h, i=1,...n-1 ;assuming p(r)=C*r**(l+1)*polynomial(r) for r==0;
    !                               p((n+1)*h)=0
    !  nz=nuclear charge
    !  emin=is estimate of lowest eigenvalue; used if nz=0
    !     otherwise, set to the value of -(nz/(l+1))**2
    !
    !  It is assumed that the wavefunction has np-l-1 nodes, where
    !    np is the principle quantum number-- np=1,2,..nroot
    !
    !  uses Noumerov algorithm
    !
    !  For l=0,1 corrections are needed to approximate wfn(r=0)
    !     These depend upon:
    !         e0 (current guess of energy eigenvalue)
    !         l,nz
    !         v(0) == v0 electronic potential at r=0
    !         v'(0) == v0p derivative of electronic potential at r=0
    !
    !  Corrections are also needed for r>n*h, depending on:
    !         e0 (current guess of energy eigenvalue
    !         the extrapolated value of rv == r * v
    !
    ! ierr=an nroot digit number indicating status of each root
    !   a digit of 1 indicates success in converging root
    !              2 indicates near success in converging root
    !              9 indicates that root not found
    !
    ! first check how many roots expected =  ntroot (returned as argument)
    !
    TYPE(GridInfo), INTENT(IN) :: Grid
    TYPE(PotentialInfo), INTENT(IN) :: Pot
    TYPE(OrbitInfo), INTENT(INOUT) :: Orbit
    INTEGER, INTENT(IN) :: l,start,nroot
    INTEGER, INTENT(INOUT) :: ierr
    REAL(8), INTENT(INOUT) :: emin

    REAL(8), PARAMETER :: convre=1.d-10,vlrg=1.d30
    INTEGER, PARAMETER :: niter=1000

    REAL(8), POINTER :: rv(:),eig(:),wfn(:,:)
    REAL(8), ALLOCATABLE :: p1(:),p2(:),a(:),b(:)
    INTEGER :: nz,n
    REAL(8) :: h,v0,v0p
    REAL(8) :: err,convrez,angm,energy,b0p0,c1,c2,c3,fact,eh2
    REAL(8) :: ssd, ssl,scale,emax,best,rout,ppp
    REAL(8) :: arg,r,r2,veff,pppp1,rin,dele,x,rvp1,pnp1,bnp1
    INTEGER :: iter,i,j,k,node,match,mxroot,ntroot,ir,iroot
    INTEGER :: least,many,ifac
    LOGICAL :: ok

    n=Grid%n
    h=Grid%h
    wfn=>Orbit%wfn
    eig=>Orbit%eig
    ALLOCATE(p1(n),p2(n),a(n),b(n),stat=i)
    IF (i/=0) THEN
       WRITE(6,*) ' Error in boundsch allocation ',i,n
       STOP
    ENDIF

    nz=Pot%nz
    v0=Pot%v0
    v0p=Pot%v0p
    rv=>Pot%rv
    err=n*nz*(h**4)
    convrez=convre
    IF (nz.GT.0) convrez=convre*nz
    !     write(6,*) 'expected error = ',err
    ierr=0

    angm=l*(l+1)
    DO  i=1,n
       a(i)=angm/(i*i)+h*rv(i)/i
    ENDDO
    DO  i=1,n
       b(i)=-1.2d0+0.1d0*a(i)
    ENDDO
    ssl=0.1d0
    ssd=1.d0
    DO  i=1,n
       a(i)=2.4d0+a(i)
    ENDDO
    WRITE(6,*) 'z , l = ',nz,l
    ! check how many roots expected by integration outward at
    !   energy = 0
    energy = 0
    !
    !  start outward integration
    !    correct behavior near r=0
    b0p0=0
    IF (l<=1) THEN
       c1=-nz/(l+1.d0)
       c2=((v0-energy)-2*nz*c1)/(4*l+6.d0)
       c3=(v0p+(v0-energy)*c1-2*nz*c2)/(6*l+12.d0)
       fact=1.d0+h*(c1+h*(c2+h*c3))
       IF (l.EQ.0) THEN
          b0p0=-0.2d0*h*nz/fact
       ENDIF
       IF (l.EQ.1) THEN
          b0p0=0.2d0/fact
       ENDIF
    ENDIF
    eh2=0
    node=0
    p1(1)=h**(l+1)
    p1(2)=(a(1)+b0p0-eh2*ssd)*p1(1)/(eh2*ssl-b(2))
    DO j=3,n
       p1(j)=((a(j-1)-eh2*ssd)*p1(j-1)+                                 &
            (b(j-2)-eh2*ssl)*p1(j-2))/(eh2*ssl-b(j))
       IF ((j.LT.n-10).AND.p1(j)*p1(j-1).LE.0.d0) node=node+1
       !  if wavefunction too large, renormalize
       IF (ABS(p1(j)).GT.vlrg) THEN
          scale=1.d0/vlrg
          p1(1:j)=p1(1:j)*scale 
       ENDIF
    ENDDO
    WRITE(6,*) ' nodes at e=0  ', node
    mxroot=node+1
    ntroot=node
    IF (mxroot.LT.nroot) THEN
       WRITE(6,*)'error in boundsch - for l = ',l
       WRITE(6,*) nroot,' states requested but only',mxroot,' possible'
       DO ir=mxroot+1,nroot
          ierr=ierr+9*(10**(ir-1))
       ENDDO
    ENDIF
    mxroot=min0(mxroot,nroot)
    !
    IF (nz.EQ.0) energy=-ABS(emin)
    IF (nz.NE.0) energy=-(nz/(l+1.d0))**2
    emin=energy-err
    emax=0.d0

    DO iroot=1,mxroot
       best=1.d10; dele=1.d10
       energy=emin+err
       IF (energy.LT.emin) energy=emin
       IF (energy.GT.emax) energy=emax
       ok=.FALSE.
            !write(6,*) 'iter,iroot,energy',iter,iroot,energy
            !write(6,*) 'emin,max',emin,emax
       BigIter: DO iter=1,niter
          !write(6,*) 'In iter with energy', iter,energy,niter,l,iroot
          eh2=energy*h*h
          !  start inward integration
          p2(1:n)=0.d0 
          !  start integration at n
          least=min0(10,n/4)
          r=n*h
          r2=r*r
          veff=rv(n)/n+rv(n-1)/(n-1)
          veff=0.5d0*veff/h+angm/r2
          ppp=SQRT(ABS(veff-energy))
          j=n
          p2(j)=1.d0/SQRT(ppp)
          rvp1=2*rv(j)-rv(j-1)
          bnp1=-1.2d0+0.1d0*(h*rvp1/(j+1)+angm/(j+1)**2)
          r=r+h
          r2=r*r
          pppp1=SQRT(ABS(rvp1/r+angm/r2-energy))
          arg=0.5d0*(ppp+pppp1)*h
          pnp1=EXP(-arg)/SQRT(pppp1)
          p2(j-1)=((a(j)-eh2*ssd)*p2(j)+                                    &
               (bnp1-eh2*ssl)*pnp1)/(eh2*ssl-b(j-1))
          j=j-1

          DO
             j=j-1
             p2(j)=((a(j+1)-eh2*ssd)*p2(j+1)+                                  &
                  (b(j+2)-eh2*ssl)*p2(j+2))/(eh2*ssl-b(j))
             !  if wavefunction too large, renormalize
             IF (ABS(p2(j)).GT.vlrg) THEN
                many=n-j+1
                scale=1.d0/vlrg
                p2(j:j+many-1)=p2(j:j+many-1)*scale
             ENDIF
             !  step inward until p2 stops increasing or j=least
             IF (.NOT.(p2(j).GT.p2(j+1).AND.j.GT.least)) EXIT
          ENDDO

          match=j+1
          rin=(p2(match+1)-p2(match-1))/(2*h*p2(match))
          !     write(6,*) ' match point = ',match,rin,p2(match)
          !  start outward integration
          !    correct behavior near r=0
          b0p0=0
          IF (l<=1) THEN
             c1=-nz/(l+1.d0)
             c2=((v0-energy)-2*nz*c1)/(4*l+6.d0)
             c3=(v0p+(v0-energy)*c1-2*nz*c2)/(6*l+12.d0)
             fact=1.d0+h*(c1+h*(c2+h*c3))
             IF (l.EQ.0) THEN
                b0p0=-0.2d0*h*nz/fact
             ENDIF
             IF (l.EQ.1) THEN
                b0p0=0.2d0/fact
             ENDIF
          ENDIF
          node=0
          p1(1)=h**(l+1)
          p1(2)=(a(1)+b0p0-eh2*ssd)*p1(1)/(eh2*ssl-b(2))
          DO j=3,match+1
             p1(j)=((a(j-1)-eh2*ssd)*p1(j-1)+                                 &
                  (b(j-2)-eh2*ssl)*p1(j-2))/(eh2*ssl-b(j))
             IF (p1(j)*p1(j-1).LE.0.d0) node=node+1
          ENDDO
          rout=(p1(match+1)-p1(match-1))/(2*h*p1(match))
          !   write(6,*) 'node,match,rin,rout',node,(iroot-1),match,rin,rout
          ! check whether node = (iroot-1)
          !   not enough nodes -- raise energy
          IF (node.LT.iroot-1) THEN
             emin=MAX(emin,energy)-err
             energy=emax-(emax-energy)*ranx()
             ifac=9
             !   too many nodes -- lower energy
          ELSEIF (node.GT.iroot-1) THEN
             IF (energy.LE.emin) THEN
                ierr=ierr+9*(10**(iroot-1))
                WRITE(6,*) 'boundsch error -- emin too high',l,nz,emin,energy
                STOP
             ENDIF
             emax=MIN(emax,energy+err)
             energy=emin+(energy-emin)*ranx()
             !   correct number of nodes -- estimate correction
          ELSEIF (node.EQ.iroot-1) THEN
             DO j=1,match
                p1(j)=p1(j)/p1(match)
                !         write(6,*) 'j,p1',j,p1(j)
             ENDDO
             DO j=match,n
                p1(j)=p2(j)/p2(match)
                !         write(6,*) 'j,p2',j,p1(j)
             ENDDO
             scale=1.d0/overlap(n,h,p1,p1)
             dele=(rout-rin)*scale
             !     write(6,*) 'energy,dele,scale',energy,dele,scale
             x=ABS(dele)
             IF (x.LT.best) THEN
                scale=SQRT(scale)
                p1(1:n)=p1(1:n)*scale   !call dscal(n,scale,p1,1)
                k=start+iroot-1
                wfn(1:n,k)=p1(1:n) !call dcopy(n,p1,1,wfn(1,iroot),1)
                eig(k)=energy
                !write(6,*) 'root',l,iroot,eig(k),emin,emax
                best=x
             ENDIF
             IF (ABS(dele).LE.convrez) THEN
                ok=.TRUE.
                !  eigenvalue found
                ierr=ierr+10**(iroot-1)
                IF (iroot+1.LE.mxroot) THEN
                   emin=energy+err
                   emax=0
                   energy=(emin+emax)/2
                   IF (energy.LT.emin) energy=emin
                   IF (energy.GT.emax) energy=emax
                   best=1.d10
                ENDIF
                EXIT BigIter
             EndIf
             If (ABS(dele).GT.convrez) THEN
                !write(6,*) 'iter with dele' , iter,dele
                energy=energy+dele
                ! if energy is out of range, pick random energy in correct range
                IF (emin-energy.GT.convrez.OR.energy-emax.GT.convrez)         &
                     energy=emin+(emax-emin)*ranx()
                ifac=2
                !write(6,*) 'continuing with iter dele', iter,dele
             ENDIF
          ENDIF
       ENDDO BigIter !iter
       IF (.NOT.ok) THEN
          ierr=ierr+ifac*(10**(iroot-1))
          WRITE(6,*) 'no convergence in boundsch',iroot,l,dele,energy
          WRITE(6,*) ' best guess of eig, dele = ',eig(iroot),best
          IF (iroot.LT.mxroot) THEN
             DO ir=iroot+1,mxroot
                ierr=ierr+9*(10**(ir-1))
             ENDDO
          ENDIF
       ENDIF
    ENDDO !iroot

   deallocate(p1,p2,a,b)
    !         write(6,*) 'returning from boundsch -- ierr=',ierr
  END SUBROUTINE Boundsch

END MODULE radialsch