c---------------------------------------------------------------
      program hydrogen
c---------------------------------------------------------------
c
c	Solve the radial Sh. eq for hydogen-like atoms
c	using Numerov method on a radial logarithmic grid
c       solution matching at the classical turning point
c
      implicit none
      integer mshx
      parameter (mshx = 2001)

      integer mesh
      integer n, l, i
      real*8  zmesh, rmax, xmin, zeta, dx, e
      real*8 r(0:mshx), sqr(0:mshx), r2(0:mshx)
      real*8 y(0:mshx), vpot(0:mshx), work(0:mshx)
c
c initialize atomic charge
c
      zeta = 1.d0
      write (*,'(a,f10.6)') ' ATOMIC CHARGE = ', zeta
c
c initialize logarithmic mesh
c
      zmesh=  zeta 
      rmax = 100.d0
      xmin = -8.d0
      dx   =  0.01d0
      call do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
c
c initialize the potential
c
      call init_pot(zeta,r,mesh,vpot)
c
c open wfc file
c
      open (7,file='wfc.out',status='unknown',form='formatted')
c
c read number of nodes (stop if nodes .lt. 0)
c
    1 continue
      write (*,'(a,$)') ' n,l >> '
      read (*,*) n,l
      if (n.lt.0) stop

      if (n.lt.l+1) then
         write(*,*) 'error in main: n.lt.l -> wrong number of nodes '
         stop
      end if
c
c solve the schroedinger equation in radial coordinates by Numerov method
c
      call solve_sheq(n,l,e,mesh,dx,r,sqr,r2,vpot,zeta,y,work)
c
c write out the eigenvalue energy to be compared with the external potential
c
      write (*,*) 'eigenvalue =',e,',  eig*(n/zeta)^2 =',e*(n/zeta)**2
      do i=0,mesh
         write (7,*) r(i),y(i)/sqr(i), y(i)*sqr(i), e
      end do
      write (7,*) 0.d0,0.d0, 0.d0, e

      go to 1

      end 
c---------------------------------------------------------------------
      subroutine solve_sheq(n,l,e,mesh,dx,r,sqr,r2,vpot,zeta,y,f)
c---------------------------------------------------------------------
c
c solve the schroedinger equation in radial coordinates by Numerov method
c
      implicit none
      integer mesh, n,l,j
      integer i, icl, nodes, ncross, kkk
      real*8 dx, ddx12, sqlhf, x2l2, zeta
      real*8 r(0:mesh), sqr(0:mesh), r2(0:mesh), 
     +       vpot(0:mesh), y(0:mesh), f(0:mesh)
      real*8 ycusp, dfcusp, fac, norm, eup, elw, e, de

      ddx12=dx*dx/12.d0
      sqlhf = (l+0.5d0)**2
      x2l2 = 2*l+2
c
c set initial lower and upper bounds to the eigenvalue
c
      eup=vpot(mesh)
      elw=eup
      do  i=0,mesh
          elw=min(elw,sqlhf/r2(i)+vpot(i))
      end do
      if (eup-elw.lt.1.d-10) then
         write (*,*) elw, eup
         write (*,*) 'solve_sheq: lower and upper bounds are equal'
         stop
      end if
      e = 0.5d0 * (elw + eup)

      kkk = 0
  1   continue ! this is the entry point for the solution at fixed energy
      kkk = kkk + 1
c
c set up the f-function and determine the position of its last
c change of sign
c f < 0 (approximatively) means classically allowed   region
c f > 0         "           "        "      forbidden   "
c
      icl = 2
      f(0) = ddx12 *( sqlhf + r2(0) * (vpot(0)-e) )
      do i=1,mesh
         f(i) = ddx12 * ( sqlhf + r2(i) *(vpot(i)-e) )
         if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i
      end do

      if(icl.ge.mesh-2) then
         write (*,*) icl, mesh
         write (*,*) 'error in solve_sheq: last change of sign too far'
         stop
      end if

      do i=0,mesh
         f(i)=1.0d0-f(i)
         y(i)=0.d0
      end do 
c
c determination of the wave-function in the first two points 
c
      nodes = n - l - 1
      y(0) = r(0)**(l+1) *(1.d0 - 2.d0*zeta*r(0)/x2l2) / sqr(0)
      y(1) = r(1)**(l+1) *(1.d0 - 2.d0*zeta*r(1)/x2l2) / sqr(1)
c
c outward integration and count number of crossings
c
      ncross=0
      do i =1, icl-1
         y(i+1)=((12.d0-10.d0*f(i))*y(i)-f(i-1)*y(i-1))/f(i+1)
         if ( y(i) .ne. sign(y(i),y(i+1)) ) ncross=ncross+1
      end do
      fac = y(icl) 
c
c check number of crossings
c
      if (ncross.ne.nodes) then
         if (kkk.gt.100) then
            write (*,*) kkk, e, elw, eup, ncross, nodes, icl
            write (*,*) ' error in solve_sheq: too many iterations'
            stop
         end if
         if (ncross.gt.nodes) then
            eup = e
         else 
            elw = e
         end if
         e = 0.5d0 * (eup+elw)
         go to 1
      end if
c
c determination of the wave-function in the last two points 
c assuming y(mesh+1) = 0 and y(mesh) = dx
c
      y(mesh) = dx
      y(mesh-1) = (12.d0-10.d0*f(mesh))*y(mesh)/f(mesh-1) 
c
c inward integration 
c
      do i = mesh-1,icl+1,-1
         y(i-1)=((12.d0-10.d0*f(i))*y(i)-f(i+1)*y(i+1))/f(i-1)
         if (y(i-1).gt.1.d10) then
            do j=mesh,i-1,-1
               y(j) = y(j)/y(i-1)
            end do
         end if
      end do
c
c rescale function to match at the classical turning point (icl)
c
      fac = fac/y(icl)
      do i=icl,mesh
         y(i) = y(i)*fac
      end do
c
c normalize on the segment
c
      norm = 0.d0
      do i=1,mesh
         norm = norm + y(i)*y(i) * r2(i) * dx
      end do
      norm = sqrt(norm)
      do i=0,mesh
         y(i) = y(i) / norm
      end do
c
c find the value of the cusp at the matching point (icl)
c
      i = icl
      ycusp = (y(i-1)*f(i-1)+f(i+1)*y(i+1)+10.d0*f(i)*y(i)) / 12.d0
      dfcusp = f(i)*(y(i)/ycusp - 1.d0)
c
c eigenvalue update using perturbation theory
c
      de = dfcusp/ddx12 * ycusp*ycusp * dx 
      if (de.gt.0.d0) elw = e
      if (de.lt.0.d0) eup = e
      e = max( min (e+de,eup),elw)

      if (kkk.gt.100) then
         write (*,*) kkk, e, de
         write (*,*) ' error in solve_sheq: too many iterations'
         stop
      end if

      if (abs(de).gt.1.d-10) go to 1

c ---- convergenge has been acheived -----
      write (*,*) ' convergence acheived at iter #',kkk,
     ,            ' de = ', de

      return
      end
c--------------------------------------------------------------------
      subroutine do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
c--------------------------------------------------------------------
c
c initialize mesh
c
      implicit none
      integer mshx, mesh
      real*8 zmesh,xmin,dx,rmax,r(0:mshx), sqr(0:mshx), r2(0:mshx)

      integer i
      real*8 x

      mesh = (log(zmesh*rmax) - xmin)/dx
      if (mesh.gt.mshx) then
         write (*,*) mesh, mshx
         write (*,*) 'error in routine do_mesh: mesh.gt.mshx '
         stop
      end if

      do i=0,mesh
         x = xmin + dx * i
         r(i)  = exp(x)/zmesh
         sqr(i)= sqrt(r(i))
         r2(i) = r(i) * r(i)
      end do
      write(*,*) ' radial mesh information:'
      write(*,'(a,f10.6,a,f10.6,a,f10.6)') 
     +          '  dx =',dx,', xmin =', xmin,', zmesh =',zmesh
      write(*,'(a,i6,a,f10.6,a,f10.6)')
     +          '  mesh =',mesh,', r(0) =',r(0),', r(mesh) =', r(mesh)
      write(*,*) 
      
      return
      end
c--------------------------------------------------------------------
      subroutine init_pot(zeta,r,mesh,vpot)
c--------------------------------------------------------------------
c
c initialize potential
c
      implicit none

      integer mesh
      real*8 zeta, r(0:mesh), vpot(0:mesh)
      integer i

      open (7,file='pot.out',status='unknown',form='formatted')
      do i =0,mesh
         vpot (i) = - 2.d0*zeta/r(i)
         write (7,*) r(i),vpot(i)
      end do
      close(7)

      return
      end
