!     This version solves the Schroedinger equation for the 
!     one-dimensional 
!     harmonic oscillator using matrix diagonalization

!     declaration of constants used in the calculations
!     define variables, where we have chosen to use atomic units for 
!     m=c=hbar=k=1.  

MODULE parameters     
  USE constants
  !  the step
  REAL(DP), PUBLIC :: step
  !  r min
  REAL(DP), PUBLIC, PARAMETER :: rmin=0.0_dp
  !  r max
  REAL(DP), PUBLIC, PARAMETER :: rmax=20.0_dp
  !  number of steps from rmin to rmin to rmax
  INTEGER, PUBLIC, PARAMETER ::  max_step=1000
  !  size of matrix to be diagonalized, do not need points at rmin and rmax
  INTEGER, PUBLIC, PARAMETER ::  n=998
END MODULE parameters

!    This module contains all constants and declarations 
!    of variables read in by the function read_data. These
!    variables are used by many functions.

MODULE constants
  INTEGER,  PARAMETER :: dp = KIND(1.0D0)
  INTEGER, PARAMETER :: dpc = KIND((1.0D0,1.0D0))
END MODULE constants


!    F90 program library, adapted from Numerical Recipes
!    All functions have been translated to F90 from F77

MODULE F90library
  USE constants
CONTAINS


  !     determine eigenvalues and eigenvectors of a real symmetric
  !     tri-diagonal matrix, or a real, symmetric matrix previously
  !     reduced by function tred2 to tri-diagonal form. On input,
  !     d[] contains the diagonal element and e[] the sub-diagonal
  !     of the tri-diagonal matrix. On output d[] contains the
  !     eigenvalues and  e[] is destroyed. If eigenvectors are
  !     desired z[][] on input contains the identity matrix. If
  !     eigenvectors of a matrix reduced by tred2() are required,
  !     then z[][] on input is the matrix output from tred2().
  !     On output, the k'th column returns the normalized eigenvector
  !     corresponding to d[k]. 
  !     The function is modified from the version in Numerical recipe.

  SUBROUTINE tqli(d,e,n,z)
    IMPLICIT NONE
    INTEGER :: n 
    REAL(DP)  :: d(n),e(n),z(n,n)
    INTEGER :: i,iter,k,l,m
    REAL(DP)  :: b,c,dd,f,g,p,r,s,one

    DO i=2,n
       e(i-1)=e(i)
    ENDDO
    one=1.
    DO l=1,n
       iter=0
       ITERATE : DO
          DO m=l,n-1
             dd=ABS(d(m))+ABS(d(m+1))
             IF (ABS(e(m))+dd == dd) EXIT
          ENDDO
          IF(m == l) EXIT ITERATE
          IF(iter == 30) STOP 'too many iterations in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.*e(l))
          r=pythag(g,one)
          g=d(m)-d(l)+e(l)/(g+sign(r,g))
          s=1.
          c=1.
          p=0.
          DO i=m-1,l,-1
             f=s*e(i)
             b=c*e(i)
             r=pythag(f,g)
             e(i+1)=r
             IF(r == 0.) THEN
                d(i+1)=d(i+1)-p
                e(m)=0.
                CYCLE ITERATE
             ENDIF
             s=f/r
             c=g/r
             g=d(i+1)-p
             r=(d(i)-g)*s+2.*c*b
             p=s*r
             d(i+1)=g+p
             g=c*r-b
             !     Omit lines from here ...
             DO k=1,n
                f=z(k,i+1)
                z(k,i+1)=s*z(k,i)+c*f
                z(k,i)=c*z(k,i)-s*f
             ENDDO
             !     ... to here when finding only eigenvalues.
          ENDDO
          d(l)=d(l)-p
          e(l)=g
          e(m)=0.
       ENDDO ITERATE
    ENDDO

  END SUBROUTINE tqli


  REAL(DP) FUNCTION pythag(a,b)
    REAL(DP)  :: a,b
    REAL(DP)  :: absa,absb
    absa=ABS(a)
    absb=ABS(b)
    IF(absa > absb) THEN
       pythag=absa*sqrt(1.+(absb/absa)**2)
    ELSE
       IF(absb == 0.) THEN
          pythag=0.
       ELSE
          pythag=absb*sqrt(1.+(absa/absb)**2)
       ENDIF
    ENDIF

  END FUNCTION pythag

  !    perform a Housholder reduction of a real symmetric matrix
  !    a[][]. On output a[][] is replaced by the orthogonal matrix 
  !    effecting the transformation. d[] returns the diagonal elements
  !    of the tri-diagonal matrix, and e[] the off-diagonal elements, 
  !    with e[0] = 0.
  !    The function is modified from the version in Numerical recipes.

  SUBROUTINE tred2(a,n,d,e)
    IMPLICIT NONE
    INTEGER :: n
    REAL(DP) :: a(n,n),d(n),e(n)
    INTEGER :: i,j,k,l
    REAL(DP) ::  f,g,h,hh,scale

    DO i=n,2,-1
       l=i-1
       h=0.
       scale=0.
       IF (l > 1) THEN
          scale=SUM(abs(a(i,1:l)))
          IF (scale == 0.) THEN
             e(i)=a(i,l)
          ELSE
             a(i,1:l)=a(i,1:l)/scale
             h=sum(a(i,1:l)**2)
             f=a(i,l)
             g=-sign(sqrt(h),f)
             e(i)=scale*g
             h=h-f*g
             a(i,l)=f-g
             f=0.
             DO j=1,l
                !     Omit following line if finding only eigenvalues
                a(j,i)=a(i,j)/h
                g=0.
                DO k=1,j
                   g=g+a(j,k)*a(i,k)
                ENDDO
                DO  k=j+1,l
                   g=g+a(k,j)*a(i,k)
                ENDDO
                e(j)=g/h
                f=f+e(j)*a(i,j)
             ENDDO
             hh=f/(h+h)
             DO j=1,l
                f=a(i,j)
                g=e(j)-hh*f
                e(j)=g
                DO k=1,j
                   a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
                ENDDO
             ENDDO
          ENDIF
       ELSE
          e(i)=a(i,l)
       ENDIF
       d(i)=h
    ENDDO
    !     Omit following line if finding only eigenvalues.
    d(1)=0.
    e(1)=0.
    DO i=1,n
       !     Delete lines from here ...
       l=i-1
       IF (d(i) /= 0.) THEN
          DO j=1,l
             g=0.
             DO k=1,l
                g=g+a(i,k)*a(k,j)
             ENDDO
             DO k=1,l
                a(k,j)=a(k,j)-g*a(k,i)
             ENDDO
          ENDDO
       endif
       !     ... to here when finding only eigenvalues.
       d(i)=a(i,i)
       !     Also delete lines from here ...
       a(i,i)=1.
       DO j=1,l
          a(i,j)=0.
          a(j,i)=0.
       ENDDO
       !     ... to here when finding only eigenvalues.
    ENDDO

  END SUBROUTINE tred2


END MODULE F90library


!  main program starts here

PROGRAM schroedinger_equation
  USE constants
  USE parameters
  USE f90library
  IMPLICIT NONE
  REAL(DP) :: potential, const1, const2
  REAL(DP), DIMENSION(:), ALLOCATABLE :: w, r
  REAL(DP), ALLOCATABLE, DIMENSION(:) :: e, d
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: z
  INTEGER :: i

  !  reserve place in memory for the arrays w, r, e, d and z

  ALLOCATE ( e(n) , d(n))
  ALLOCATE ( w(0:max_step), r(0:max_step))

  !  define the step size

  step=(rmax-rmin)/FLOAT(max_step)

  !  define constants for the matrix to be diagonalized

  const1=1.0/(step*step)
  const2=-0.5/(step*step)

  !     set up r and the function w for energy =0
  !     w corresponds then to the potential
  !  values at 

  DO i=0, max_step
     r(i) = rmin+i*step
     w(i) = potential(r(i))
  ENDDO

  !     setup the diagonal d and the non-diagonal part e  of
  !     the  tri-diagonal matrix to be diagonalized

  d(1:n)= const1+w(1:n)  ;  e = const2

  !  allocate space for eigenvector info

  ALLOCATE ( z(n,n) )

  !  obtain the eigenvalues

  CALL tqli(d,e,n,z)

  !  sort eigenvalues as an ascending series 

  CALL eigenvalue_sort(d,n)

  !  write out 

  WRITE(6,*) ' Rmin and Rmax ', rmin, rmax
  WRITE(6,*) ' Number of steps -1 ', n
  WRITE(6,*) ' The 5 lowest eigenvalues '
  WRITE(6,'(5F12.6)') d(1:5)



  DEALLOCATE (z)
  DEALLOCATE ( w, r, d, e)

END PROGRAM schroedinger_equation

!  Harmonic oscillator potential

REAL(DP) FUNCTION potential(x)
  USE constants
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: x

  potential=0.5*x*x

END FUNCTION  potential


!
!           Sort the eigenvalues as
!           ascending series 
!

SUBROUTINE eigenvalue_sort(e,n)
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n
  INTEGER :: i, j
  REAL(DP), DIMENSION(n), INTENT(INOUT) :: e
  REAL(DP) :: temp

  DO i = 1, n
     DO j =i, n
        IF( ABS(e(i)) >   ABS(e(j)) ) THEN
           ! change the energy  
           temp=e(i)
           e(i)=e(j)
           e(j)=temp
        ENDIF
     ENDDO
  ENDDO

END SUBROUTINE eigenvalue_sort



