	program harmonic1
*
*	Solve the quantum harmonic oscillator problem.
*	Forward and Backward integration by the Numerov method.
*	Eigenvalue search by the 'shooting method' aand bisection
*	Solution matching at the calssical turning point.
*
	implicit none

*	   Maximum dimension of the grid 
	integer mshx
	parameter (mshx = 2000)

	real*8 pi
	parameter (pi=3.14159265358979323846d0)

	integer mesh, i, icl
	integer nodes, hnodes, ncross, kkk, iterate
	real*8 xmax, dx, ddx12, xmcl, norm, yicl, djump
	real*8 eup, elw, e
	real*8 x(0:mshx), y(0:mshx), p(0:mshx), vpot(0:mshx), f(0:mshx)
	character*80 fileout

*
*	Input and grid initialization
*
	print '(a,$)','maximum x (typical value: 10) > '
	read*, xmax
	print '(a,i5,a,$)',
     $	    'Number of grid points (max:',mshx,') > '
	read*,mesh
	if (mesh.gt.mshx) stop
	dx =  xmax/mesh 
	ddx12 = dx*dx/12.d0
*
*	Definition of the potential (symmetric around x=0)
*
	do i = 0, mesh
	   x(i) = dfloat(i) * dx
	   vpot(i) = 0.5d0 * x(i)*x(i)
	enddo

	print '(a,$)','Output file name = '
	read '(a)',fileout
	open(1,file=fileout,form="formatted",status="new")

*
*	Entry point of the eigenvalue search loop
*
 999  continue
*
*	Input number of nodes (to terminate input <0)
*
	print '(a,$)', 'Number of nodes (-1=terminate) > '
	read*,nodes
	if (nodes.lt.0) then
	   close(1)
	   stop 
	endif
*
*	Define upper and lower bounds for the eigenvalue
*
	eup = vpot(mesh)
	elw = eup
	do  i=0,mesh
	   elw = min(elw,vpot(i))
	   eup = max(eup,vpot(i))
	enddo
*
*	Define trial energy 
*
	print '(a,$)','Trial energy (0=bisction search ) > '
	read*,e
	if ( e .eq. 0.d0 ) then
*	   Automatic search mode
	   e = 0.5d0 * (elw + eup)
	   iterate = 1
	else
*	   Test of a trial energy mode
	   iterate = 0
	endif

	kkk = 0
*
*	Entry point for the sh.eq. integration at fixed energy e
*
  1	continue
	kkk = kkk + 1
	djump = 0.d0
*
*	Define function g(x) used in the Numerov method
*
	f(0) = ddx12*(2.d0*(vpot(0)-e))
	do i=1,mesh
*		f <= 0  classically allowed region
*		f > 0   classically forbidden region
	   f(i)=ddx12*2.d0*(vpot(i)-e)
*		No f(i) can be exactly zero: sign change would not be detected.
*		Next line avoids the case that f(i) is exactly 0:
	   if( f(i) .eq. 0.d0 ) f(i)=1.d-20
*		Record the index of the last sign change
	   if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i
	enddo

	if(icl.ge.mesh-2) then
	   print*,'last sign change is too close to the edge.'
	   stop
	endif

*
*	f required by the Numerov method.
*
	do i=0,mesh
	   f(i) = 1.0d0 - f(i)
	enddo 

	do i=0,mesh
	   y(i) = 0.d0
	enddo
*
*	Determine the wavefunction in the first two points.
*
	hnodes = nodes/2
	if (2*hnodes.eq.nodes) then
	   y(0) = 1.0d0
	   y(1) = 0.5d0*(12.d0-10.d0*f(0))*y(0)/f(1)
	else
	   y(0) = 0.d0 
	   y(1) = dx
	endif
*
*	Outward integration up to the classical turning point,
*       counting the number of sign changes.
*
	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
	enddo
	yicl = y(icl)
*
*	Verify the number of sign changes
*
	if ( iterate .ne. 0 ) then
	   if (ncross.ne.hnodes) then
	      print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e
*		Wrong number of sign changes
	      if (ncross.gt.hnodes) then
*		   More sign changes than desired.
*		   Energy is too high: update upperbound
	         eup = e
	      else 
*                  Less sign changes that desired.
*		   Energy is too low: update lowerbound
	         elw = e
	      endif
*		   New trial energy
	      e = 0.5d0 * (eup+elw)
*	           Iterate again.
	      go to 1
	   endif
	endif
*
*	The number of nodes is correct. let's comlete the integration
*	and use left and right derivative mismatch to update the trial energy.
*
	y(mesh) = dx
	y(mesh-1) = (12.d0-10.d0*f(mesh))*y(mesh)/f(mesh-1)
*
*	Inward integration
*
	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)
	enddo
*
*	Riscale and ...
*
	yicl = yicl/y(icl)
	do i=icl,mesh
	   y(i) = y(i)*yicl
	enddo
*
*	... Rinormalize the wavefunction so that 
*	int |psi|^2 dx = 1:  
*
	norm = 0.d0
	do i=1,mesh
	   norm = norm + y(i)*y(i)
	enddo
	norm = dx * (2.d0 * norm + y(0)*y(0))
	norm = sqrt(norm)
	do i=0,mesh
	   y(i) = y(i) / norm
	enddo
*
*	Derivative mismatch at the classical turning point
*
	if ( iterate .ne. 0 ) then
	   i = icl
*	      Compute y'(i;RIGHT) - y'(i;LEFT)
	   djump = ( y(i+1) + y(i-1) - (14.d0-12.d0*f(i))*y(i) ) / dx
	   print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e
	   if (djump*y(i).gt.0.d0) then
*	         Too high energy : update upperbound
	      eup = e
	   else
*	         Too low energy : update lowerbound
	      elw = e
	   endif
	   e = 0.5d0 * (eup+elw)
	   if (eup-elw .gt. 1.d-10) go to 1
	endif
	print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e

*
*	Convergence has been acheived (or a single integrationn was required)
*

*
*	cassical density calculation (for comparison):
*
	xmcl = sqrt(2.d0*e)
	do i=icl,mesh
	   p(i) = 0.d0
	enddo
	do i=0,icl-1
	   p(i) = 1.d0/sqrt(xmcl**2 - x(i)**2)/pi
	enddo

*	fix the negative semiaxis:
	do i=mesh,1,-1
	   write (1,'(f7.3,3e16.8,f12.6)') 
     $        -x(i),(-1)**nodes*y(i),y(i)*y(i),p(i),vpot(i)
	enddo
*	positive semiaxis
	do i=0,mesh
	   write (1,'(f7.3,3e16.8,f12.6)') 
     %         x(i),y(i),y(i)*y(i),p(i),vpot(i)
	enddo

	go to 999

	end
