      program hydrogen
*
*	Risolve l'equazione radiale per atomi idrogenoidi
*	utilizzando il metodo di Numerov su una griglia logaritmica.
*	Matching delle soluzioni a un punto di inversione.
*

      implicit none
*          Numero massimo di punti nella griglia
      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)
*	Il numero atomico e' fissato qui:
      parameter (zeta = 1.d0)
      character*80 fileout

      print '(/,a,/)','Hydrogen'
*
*	Inizializzazione griglia logaritmica
*
      zmesh=  zeta 
      rmax = 100.d0
      xmin = -8.d0
      dx   =  0.01d0
      call do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
*
*	Inizializzazione potenziale
*
      call init_pot(zeta,r,mesh,vpot)
*
*	Definizione file di output
*
      print '($,a)','Nome file output = '
      read '(a)',fileout
      open (7,file=fileout,status='unknown',form='formatted')
*
*	Legge n e l (da cui il numero di nodi (n-l-1)):
*
    1 continue
      print '($,a)', 'n,l (n<0 per uscire) ? '
      read (*,*) n,l
      if (n.lt.0) stop

      if (n-1.lt.l) then
         print*, 'l > n-1 e'' impossibile.'
         stop
      end if
*
*	Risoluzione dell'equazione radiale col metodo di Numerov
*
      call solve_sheq(n,l,e,mesh,dx,r,sqr,r2,vpot,zeta,y,work)
*
*	Autovalore, e verifica [dovrebbe essere e=-(zeta/n)^2 ]
*
      print '(a,f15.10,a,f15.10)',
     $      'e =',e,' ,  eig*(n/zeta)^2 =',e*(n/zeta)**2
*
*	Output su file: r, R(r), Chi(r)=rR(r), |Chi(r)|^2
*	V(r), V(r)+l(l+1)/r^2, E
*
      do i=0,mesh
         write (7,'(f9.4, e14.6, f11.7, f11.8, 2e12.4, f10.6)') 
     $      r(i), y(i)/sqr(i), y(i)*sqr(i), y(i)*y(i)*r(i),
     $      vpot(i), vpot(i) + l*(l+1)/r2(i), e
      end do

      go to 1

      end 

*---------------------------------------------------------------------
      subroutine solve_sheq(n,l,e,mesh,dx,r,sqr,r2,vpot,zeta,y,f)
*---------------------------------------------------------------------
*
*	Risoluzione dell'equazione di Schroedinger in coordinate
*	radiali col metodo di Numerov.
*
      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, yicl, norm, eup, elw, e, de

      ddx12=dx*dx/12.d0
      sqlhf = (l+0.5d0)**2
      x2l2 = dble(2*l+2)
*
*       Definizione limiti all'autovalore
*
      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
         print*, elw, eup
         print '(a)',
     $       'solve_sheq: i limiti inferiore e superiore sono uguali'
         stop
      end if
      e = 0.5d0 * (elw + eup)

      kkk = 0
*
*       Punto di entrata per calcolo soluzione a energia e
*
  1   continue
      kkk = kkk + 1
*
*       Definizione funzione g usata nell'algoritmo di Numerov
*
      icl = 2
      f(0) = ddx12 *( sqlhf + r2(0) * (vpot(0)-e) )
      do i=1,mesh
*               f <= 0  regione permessa classicamente
*               f > 0   regione proibita classicamente
         f(i) = ddx12 * ( sqlhf + r2(i) *(vpot(i)-e) )
*               Nessun f(i) puo' essere zero altrimenti
*               il cambio di segno non viene rilevato.
*               La linea seguente e' un trucco per evitare
*               di imbattersi in f(i) nulli:
         if( f(i) .eq. 0.d0 ) f(i)=1.d-20
*               Registra l'indice a cui e' avvenuto 
*               l'ultimo cambio di segno:
         if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i
      end do

      if(icl.ge.mesh-2) then
         print*, icl, mesh
         print '(a)', 
     $       'Errore in solve_sheq: ultimo cambio segno troppo lontano'
         stop
      end if

*
*       f richiesta dal metodo di Numerov:
*
      do i=0,mesh
         f(i) = 1.0d0 - f(i)
         y(i) = 0.d0
      end do 
*
*       Determinazione funzione d'onda nei primi due punti
*
      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)
*
*       Integrazione verso l'esterno fermandosi all'ultimo punto
*       di inversione, e conteggio dei cambi di segno
*
      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
      yicl = y(icl) 
*
*       Verifica il numero di cambi di segno
*
      if (ncross.ne.nodes) then
*           Numero errato di cambi di segno
         if (kkk.gt.100) then
            print*, kkk, e, elw, eup, ncross, nodes, icl
            print '(a)', 'Errore in solve_sheq: troppe iterazioni'
            stop
         end if
         if (ncross.gt.nodes) then
*              Avvenuti piu' cambi di segno del previsto.  
*              L'energia e' troppo alta e va abbassata 
*              (--> semiintervallo inferiore).
            eup = e
         else 
*              Avvenuti meno cambi di segno del previsto.
*              L'energia e' troppo bassa e va alzata 
*              (--> semiintervallo superiore).
            elw = e
         end if
*            Nuovo valore da provare:
         e = 0.5d0 * (eup+elw)
*            Riiteriamo comunque (il numero di nodi era errato)
         go to 1
      end if
*
*       Il numero di nodi e' giusto.  Completiamo l'integrazione
*       e utilizziamo il matching della derivata al punto di
*       giunzione come criterio per la ricerca dell'autovalore.
*
*	Determinazione della funzione d'onda negli ultimi due punti
*	assumendo y(mesh+1) = 0 e y(mesh) = dx
*
      y(mesh) = dx
      y(mesh-1) = (12.d0-10.d0*f(mesh))*y(mesh)/f(mesh-1) 
*
*       Integrazione dall'esterno verso l'interno
*
      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
*
*       Riscala per ottenere matching a icl:
*
      yicl = yicl/y(icl)
      do i=icl,mesh
         y(i) = y(i)*yicl
      end do
*
*       Rinormalizzazione della funzione d'onda in modo che
*       int r^2 |chi|^2 dr = 1:  
*
      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
*
*       Calcolo discontinuita' della derivata prima al punto di
*       raccordo
*
      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)
*
*	Aggiornamento autovalore usando teoria delle perturbazioni
*
      de = dfcusp/ddx12 * ycusp*ycusp * dx 
*         Il nuovo semiintervallo e' quello superiore?
      if (de.gt.0.d0) elw = e
*         Il nuovo semiintervallo e' quello inferiore?
      if (de.lt.0.d0) eup = e
*         La linea seguente definisce e=e+de, con "clipping" a
*         elw e eup nel caso che e+de esca da [elw,eup]
      e = max(min(e+de,eup),elw)

      if (kkk.gt.100) then
         print*, kkk, e, de
         print '(a)', 'Errore in solve_sheq: troppe iterazioni'
         stop
      end if

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

*
*       La convergenza e' stata raggiunta
*
      print '(a,i4,a,e12.4)', 'Convergenza raggiunta a iterazione #',
     $     kkk, ' , de = ', de

      end


*--------------------------------------------------------------------
      subroutine do_mesh(zmesh,xmin,dx,rmax,mshx,mesh,r,sqr,r2)
*--------------------------------------------------------------------
*
*	Inizializzazione della griglia logaritmica
*
      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
         print*, mesh, mshx
         print*, 'Errore in routine do_mesh: mesh > 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
      print '(a)','Informazioni sulla griglia radiale:'
      print '(a,f7.3)','Z =',zmesh
      print '(a,f10.6,a,f10.6,a,f10.6)',
     $          '  dx =',dx,', xmin =', xmin,', zmesh =',zmesh
      print '(a,i6,a,f10.6,a,f10.6)',
     $          '  mesh =',mesh,', r(0) =',r(0),', r(mesh) =', r(mesh)
      print*
      
      end

*--------------------------------------------------------------------
      subroutine init_pot(zeta,r,mesh,vpot)
*--------------------------------------------------------------------
*
*	Inizializzazione del potenziale coulombiano
*
      implicit none

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

      do i =0,mesh
         vpot (i) = - 2.d0*zeta/r(i)
      end do
      end
