!--------------------------------------------------------------------------------!
! - Solves the Schrodinger equation in one dimension, for a potential given      !
!   by the function 'potential(x)', for x in the range (-xmax,xmax).             !
! - The  integration of the wave functipn Psi(x) starts from boundary conditions !
!   applied to Psi(-xmax) and Psi(-xmax+h), where h is the integration step      !
!   (this boundary condition is set in the subroutine 'setinitcond').            !
! - The boundary condition at x=xmax is taken to be Psi(xmax)=0.                 !
! - Bisection is used to find the energy (preceded by a course search            !
!   for two energies bracketing a solution).                                     !
!--------------------------------------------------------------------------------!
 program schrodinger1d
!---------------------!
 implicit none

 real(8), parameter :: eps=1.d-6

 integer :: i,ni,nx
 real(8) :: xmax,de,b0,b1,b2,e0,e1,e2

 real(8), allocatable :: psi(:)

 write(*,*)'integrating between -xmax and xmax; give xmax:'
 read(*,*)xmax
 write(*,*)'Number of points on either side of 0'
 read(*,*)nx
 write(*,*)'searching for E>E0; give E0'
 read(*,*)e1
 write(*,*)'Delta-E in initial course search'
 read(*,*)de
 
 allocate(psi(-nx:nx))

 ni=0
 write(*,*)'Starting course search'
 call integrate(xmax,nx,e1,psi)
 b1=psi(nx)
 call writemessage(ni,e1,b1) 
 do 
   ni=ni+1
   e2=e1+de
   call integrate(xmax,nx,e2,psi)
   b2=psi(nx)
   call writemessage(ni,e2,b2) 
   if (b1*b2 < 0.d0) exit
   e1=e2
   b1=b2
 enddo

 ni=0
 write(*,*)'Starting bisection'
 do
   ni=ni+1
   e0=(e1+e2)/2.0d0
   call integrate(xmax,nx,e0,psi)
   b0=psi(nx)
   call writemessage(ni,e0,b0) 
   if (abs(b0) <= eps) exit
   if (b0*b1 <= 0.0d0) then
     e2=e0      
     b2=b0 
   else
     e1=e0
     b1=b0
   endif
 enddo

 open(1,file='psi.dat',status='replace')
 do i=-nx,nx
   write(1,*)dble(i)*xmax/dble(nx),psi(i)
 enddo
 close(1)

 end program schrodinger1d
!-------------------------!

!------------------------------------!
 subroutine integrate(xmax,nx,e1,psi)
!------------------------------------!
 implicit none

 real(8) :: q0,q1,p1,f1,ee,h2,h12
 common/block1/q0,q1,p1,f1,ee,h2,h12

 integer :: i,n,nx
 real(8) :: psi(-nx:nx),xmax,x,h,de,e1

 ee=e1
 h=xmax/dble(nx)      
 h2=h**2
 h12=h2/12.d0
 call setinitcond(xmax,h,psi(-nx),psi(-nx+1))
 do i=-nx+2,nx
   x=dble(i)*h
   call numerovstep(x)
   psi(i)=p1
 end do 
 call normalize(psi,nx,h)

 end subroutine integrate 
!------------------------!

!----------------------------------------------------------------!
!Integrates the wave function one step, using Numerov's method.  !
!On entry, the wave-function at the previous step is p1, on exit !
!it is the wave-function after the integration step. The modified! 
!wave function, Psi(x)*[1-(h**2)*f1], at the two previous steps  !
!are q0,q1 on entry, and f1=2.d0*(potential(x1)-energy).         !
!----------------------------------------------------------------!
 subroutine numerovstep(x)
!---------------------------!
 implicit none

 real(8) :: q0,q1,p1,f1,ee,h2,h12
 common/block1/q0,q1,p1,f1,ee,h2,h12
       
 real(8) :: x,q2,potential

 q2=h2*f1*p1+2.d0*q1-q0
 q0=q1; q1=q2
 f1=2.d0*(potential(x)-ee)        
 p1=q1/(1.d0-h12*f1)
         
 end subroutine numerovstep
!--------------------------!

!---------------------------------------------------!
!Potential energy as a function of the coordinate x !
!Three examples are given (two commented out).      !
!---------------------------------------------------!
 real(8) function potential(x)
!-----------------------------!
 implicit none

 real(8) :: x

!----- square well with infinite walls       
! potential=0.d0

!----- square well with soft walls V(|x|>1)=10
!if (abs(x) < 1.d0) then
!  potential=0.d0
!else
!  potential=10.d0
!endif 

!----- square well with hard walls internal and box between -0.5 and 0.5
 if (x > -0.5d0 .and. x < 0.5d0) then
   potential=5.d0
 else
   potential=0.d0
 endif 

 end function potential
!----------------------!

!----------------------------------------------------------------!
!Sets the boundary conditions Psi(-xmax)=0, Psi(-xmax+h)=0.0001  !
!- note that the boundary condition is strictly not correct for  !
!  a "sof" potential, but OK if xmax is large enough             !
!----------------------------------------------------------------!
 subroutine setinitcond(xmax,h,psi0,psi1)
!----------------------------------------!
 implicit none

 real(8) :: q0,q1,p1,f1,ee,h2,h12
 common/block1/q0,q1,p1,f1,ee,h2,h12

 real(8) :: h,vv,xmax,psi0,psi1,potential

 psi0=0.d0
 psi1=0.0001d0

 p1=psi0
 f1=2.d0*(potential(-xmax)-ee)        
 q0=psi0*(1.d0-h12*f1)
 f1=2.d0*(potential(-xmax+h)-ee)        
 q1=psi1*(1.d0-h12*f1)
	
 end subroutine setinitcond
!--------------------------!

!------------------------------!
 subroutine normalize(psi,nx,h)
!------------------------------!
 implicit none

 integer :: i,nx
 real(8) :: psi(-nx:nx),norm,h

 norm=psi(-nx)**2+psi(nx)**2
 do i=-nx+1,nx-3,2
   norm=norm+4.d0*psi(i)**2+2.d0*psi(i+1)**2
 end do
 norm=norm+4.d0*psi(nx-1)**2
 norm=1.d0/sqrt(norm*h/3.d0)
 do i=-nx,nx
   psi(i)=psi(i)*norm
 end do

 end subroutine normalize
!------------------------!

!------------------------------!
 subroutine writemessage(n,e,b) 
!------------------------------!
 implicit none

 integer :: n
 real(8) :: e,b

 write(*,1)n,e,b
 1 format(i4,'    E = ',f10.7,'    Boundary deviation ',f10.7)

 end subroutine writemessage
!---------------------------!