!===========================================================
!
! File: sch.f90
!
! Integrate 1d Schrodinger equation from xmin to xmax. 
! Determine energy eigenvalue and eigenfunction by matching 
! evolving solutions from xmin and from xmax at a point xm. 
! Matching done by equating values of functions and their 
! derivatives at xm. The point xm chosen at the left most 
! turning point of the potential at any given value of the
! energy. The potential and boundary conditions chosen in 
! different file.
! ----------------------------------------------------------
! Input:  energy: Trial value of energy
!         de: energy step, if matching fails de -> e+de, if
!             logderivative changes sign     de -> -de/2
!         xmin, xmax, Nx
! ----------------------------------------------------------
! Output: Final value of energy, number of nodes of 
!    wavefunction in stdout
!    Final eigenfunction in file psi.dat
!    All trial functions and energies in file all.dat
! ----------------------------------------------------------
! Compile:
!    gfortran schInfSq.f90 sch_util.f90  sch.f90 -o s
!===========================================================
program schroedinger_equation_1D
 use                   :: schV
 use                   :: schutil
 implicit none 
 real(dp), allocatable :: psi(:),psip(:)
 integer  ::  Nx,NxL,NxR
 real(dp) ::  dx
 real(dp) ::  xmin,xmax,xm       !left/right/matching points
 real(dp) ::  psixmin,psipxmin,psixmax,psipxmax
 real(dp) ::  psileft ,psiright ,psistep,psinorm
 real(dp) ::  psipleft,psipright,psipstep
 real(dp) ::  de,epsilon
 real(dp) ::  matchlogd,matchold,psiold,norm,x
 integer  ::  iter,i,imatch,nodes
 integer  ::  uout               !file unit for output
!---------- Input:
 print *,'# Enter energy,de,xmin,xmax,Nx'
 read  *,         energy,de,xmin,xmax,Nx
!--- need even intervals for normalization integration
 if( xmin      >= xmax) stop 'Error: xmin >= xmax'
 if( mod(Nx,2) == 0   ) Nx=Nx+1 
 allocate(psi(Nx),psip(Nx))
 dx      = (xmax - xmin)/(Nx-1)
 epsilon = 1.0e-6_dp
 call boundary(xmin,xmax,psixmin,psipxmin,psixmax,psipxmax)
 print *,'# #######################################'
 print *,'# Estart= ',energy, ' de= ',de
 print *,'# Nx=  ',Nx ,' eps= ',epsilon
 print *,'# xmin= ',xmin,' xmax= ',xmax, ' dx= ',dx
 print *,'# psi(xmin)= ',psixmin,' psip(xmin)= ',psipxmin
 print *,'# psi(xmax)= ',psixmax,' psip(xmax)= ',psipxmax
 print *,'# #######################################'
!----- Calculate:
 open(newunit=uout,file='all.dat')
 matchold   =  0.0_dp
 do iter    =  1, 10000
!----- Determine matching point at turning point from the left:
  imatch    = -1
  do i      =  1, Nx
   x        = xmin + (i-1)*dx
   if( imatch < 0    .and. (energy-V(x)) >  0.0_dp) imatch = i 
  enddo
  if ( imatch <= 100 .or.   imatch       >= Nx-100) imatch = Nx/5
  xm        = xmin + (imatch-1)*dx
  NxL       = imatch
  NxR       = Nx-imatch+1
!----- Evolve wavefunction from the left:
  psi  (1)  = psixmin
  psip (1)  = psipxmin
  psistep   = psixmin
  psipstep  = psipxmin
  do i=2,NxL
   x        = xmin + (i-2)*dx !this is x before the step
   call RKSTEP(x,psistep,psipstep, dx)
   psi (i)  = psistep
   psip(i)  = psipstep
  enddo
! use this to normalize eigenfunction to match at xm
  psinorm   = psistep      
  psipleft  = psipstep
!----- Evolve wavefunction from the right:
  psi (Nx)  = psixmax
  psip(Nx)  = psipxmax
  psistep   = psixmax
  psipstep  = psipxmax
  do i=2,NxR
   x        = xmax - (i-2)*dx
   call RKSTEP(x,psistep,psipstep,-dx)
   psi (Nx-i+1) = psistep
   psip(Nx-i+1) = psipstep
  enddo
  psinorm       = psistep/psinorm
  psipright     = psipstep
!----- Renormalize psil so that psil(xm)=psir(xm)
  psi (1:NxL-1) = psinorm * psi (1:NxL-1)
  psip(1:NxL-1) = psinorm * psip(1:NxL-1)
  psipleft      = psinorm * psipleft
!----- print current solution:
  do i          = 1, Nx
   x            = xmin + (i-1)*dx
   write(uout,*)iter,energy,x,psi(i),psip(i)
  enddo
!----- matching using derivatives:
!Careful: this can fail if psi'(xm) = 0 !! (use also |de|<1e-6
!criterion)
  matchlogd     = &
       (psipright-psipleft)/(ABS(psipright)+ABS(psipleft))
  print *,'# iter,energy,de,xm,logd: ',&
       iter,energy,de,xm,matchlogd
!----- Exit condition:
  if(ABS(matchlogd) <= epsilon.or.ABS(de/energy) < 1.0e-12_dp)&
       exit
  if( matchlogd * matchold < 0.0_dp) de = -0.5_dp*de
  energy        = energy + de
  matchold      = matchlogd
 enddo ! do iter=1,10000
 close(uout)
!---------------------------------------------------
!----- Solution has been found and now it is stored:
 norm     = integrate(psi,dx,Nx)
 norm     = 1.0D0/sqrt(norm)
 psi(:)   = norm*psi(:)
!----- Cound number of zeroes, add one and get energy level:
 nodes    = 1
 psiold   = psi(1)
 do i=2,Nx-1
  if(  ABS(  psi(i))  > epsilon)then !should be 0 within epsilon
   if( psiold*psi(i)  < 0.0_dp )nodes = nodes+1
   psiold =   psi(i)
  endif
 enddo !i=2,Nx-1
!------- Print final solution:
 open(newunit=uout,file='psi.dat')
 print *,'Final result: E= ',energy,' n= ',nodes,&
      ' norm = '      ,norm
 if( ABS(matchlogd) .gt. epsilon) print *&
      ,'Final result: SOS: logd>epsilon. logd= ',matchlogd
 write (uout,*)'# E= '         ,energy,' n= ',nodes,&
      ' norm = '      ,norm
 do i     = 1, Nx
  x       = xmin + (i-1)*dx
  write(uout,*) x,psi(i)
 enddo
 close (uout)
end program schroedinger_equation_1D
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------
