!     This program computes the wave function numerically for a 
!     one particle problem using Numerov's algorithm to solve
!     the homogeneous  Schroedinger equation for Hydrogen 
!     with atomic units. The wave function is compared with the 
!     theoretical radial part. 
!     For 1s this is R(r) =  (2/a_0^3/2)*exp(-Z*r/a_0)
!     For 2s it is   R(r) =  2/(2a_0)^3/2)*(1-r/2)*exp(-Z*r/2)
!     For 2p it is   R(r) =  1/(sqrt(3)*(2a_0)^3/2)*r/a_0*exp(-Z*r/2)
!     
!
!     The homogenous radial equation you solve  is rewritten  for u(r) = r*R(r)
!     The program solves 
!     [-d^2/dr^2 + l(l+1)/r*r  - 2Z/r -  energy] u(r)  = 0  using Numerov's method
!
!     You need to specify  the number of integration points, the orbital momentum -l-, 
!     the principal quantum number -n- and the maximum radius, in atomic units, for the integration.
!     Typically, the 1s wave vanishes at 5a_0 (a_0= Bohr radius = 0.05 nm), while the 2s vanishes at 15a_0.   
!
!     declaration of constants used in the calculations


MODULE parameters     
  USE constants
  !  r min
  REAL(DP), PUBLIC, PARAMETER :: rmin =  1.0E-10
  !  r max and the step length
  REAL(DP), PUBLIC :: rmax, step
  !  number of steps from rmin to rmax, charge, l and n
  INTEGER, PUBLIC ::  max_step, n, l, charge

END MODULE parameters

!  main program starts here

PROGRAM schroedinger_equation
  USE constants
  USE parameters
  USE f90library
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), ALLOCATABLE :: r
  REAL(DP) :: energy
  INTEGER :: i

  !  Read from screen  rmax, max_step, n, l and charge
!  READ(5,*) ; READ(5,*) rmax
!  READ(5,*); READ(5,*) max_step
!  READ(5,*); READ(5,*) n
!  READ(5,*); READ(5,*) l
!  READ(5,*); READ(5,*) charge

!  quick and dirty setup of parameters
  rmax = 30.  ;  max_step = 1000;  n = 1; l = 0; charge = 1

  IF ( n < 1) WRITE(6,*) 'n has a non-allowed value, choose n > 0'
  IF ( l < 0) WRITE(6,*) 'l has a non-allowed value, choose l >= 0'
  !  reserve place in memory for the array r
  ALLOCATE (r(0:max_step))
  !  define the step size
  step=(rmax-rmin)/max_step
  !     set up r 
  DO i=0, max_step
     r(i) = rmin+i*step
  ENDDO
  energy = -1.0_dp/(2.0_dp*n*n) 
  CALL wavefunction(energy, r )
  DEALLOCATE (r)

END PROGRAM schroedinger_equation

!  This function computes  the wave function for a specific energy in atomic units
!  The homogeneous ODE [-d^2/dr^2 + l(l+1)/r*r  - 2Z/r -  energy] u(r)  = 0  
!  with boundary conditions u(0) = u(infty) = 0  is rewritten
!  [-d^2/dr^2 + k(r)] u(r)  = 0, with k(r) =  l(l+1)/r*r  - 2Z/r -  energy
!  It employ's Numerov's method to integrate the differential equation
!  You will need to replace the simple matching approach with the calculation
!  of the Green's function.

SUBROUTINE wavefunction(energy, r )
  USE constants
  USE parameters
  IMPLICIT NONE
  INTEGER :: i, match
  REAL(DP), INTENT(IN) :: energy
  REAL(DP), DIMENSION(0:max_step), INTENT(IN) :: r
  REAL(DP), DIMENSION(:), ALLOCATABLE :: k, wf 
  REAL(DP) :: wwf, fac, norm

  ALLOCATE ( k(0:max_step), wf(0:max_step))
  wf = 0.0_dp
  !  add the chosen energy to the potential and set up the homogeneous part
  k = -(l*(l+1)/(r*r)-2*energy-2*charge/r)*step*step
  !  Boundary conditions
  !  value of the wave function at max step
  wf(max_step)=0.0_dp
  wf(max_step-1)  = 1.0E-10
  !  start from right and move towards rmin
  OUTER_TO_INNER : DO
     DO i=max_step-2, 1, -1
        wf(i)=      (2*(1-5*k(i+1)/12)* wf(i+1)- (1+k(i+2)/12)* wf(i+2))/(1+k(i)/12)
        !  find point were the outer part starts turning in order to determine matching point
        IF ( ABS(wf(i)) <=  ABS(wf(i+1)) ) EXIT OUTER_TO_INNER  
     ENDDO
  ENDDO OUTER_TO_INNER
  match = i+1
  wwf=wf(match)
  !  start from rmin  and move towards the matching point
  wf(0) = 0.0_dp
  wf(1) = r(1)**(l+1)
  DO i=2, match
     wf(i)=(2*(1-5*k(i-1)/12)* wf(i-1)- (1+k(i-2)/12)* wf(i-2))/(1+k(i)/12)
  ENDDO

  !  scale the wave function below the matching point
  fac=wwf/wf(match)
  write(6,*) 'matching point', match, wf(match+1), fac
  wf(0:match)=wf(0:match)*fac
  !  normalize the wave function
  norm = step*SUM(wf*wf)
  IF ( norm > 0.0_dp ) norm = 1.0_dp/SQRT(norm)
  wf = wf*norm
  !  write out of numerically calculated and analytical expression for the wave functions, here 1s
  WRITE(6,'(E12.6,2X,E12.6,2X,E12.6)')( r(i), wf(i), 2*r(i)*EXP(-charge*r(i)),i=0, max_step)
  DEALLOCATE ( wf, k) 

END SUBROUTINE wavefunction

