program sqwelltime

      implicit none

      INTERFACE

         FUNCTION f(x)
            REAL*8 f
            REAL*8, INTENT(in) :: x
         END FUNCTION f

         FUNCTION trapezoid(B,n,dx)
            complex*16 , dimension(:), intent(in) :: B
            integer, intent(in) :: n
            real*8, intent(in) :: dx
            real*8 trapezoid
         END FUNCTION trapezoid

      END INTERFACE

      integer, parameter :: n=400
      integer, parameter :: DP = SELECTED_REAL_KIND(14)
      real (kind=DP) :: pi
      integer :: i,j,INFO
      integer , dimension(n) :: IPIV
      complex*16 , dimension(n,n) :: T, Tstar
      complex*16 , dimension(n) :: B
      character*15 :: filename
      real (kind=DP) :: x, h, h2, p, alpha, bsum, dt

      pi = DATAN(1._DP)*4._DP
      h = 1._DP/dble(n+1)
      h2 = h*h
      dt = .45_DP*h2
      alpha = 2._DP*h*h/dt

      T = (0._DP, 0._DP)
      x = h
      T(1,1) = DCMPLX(2._DP + f(x)*h2 , -alpha)
      T(1,2) = (-1._DP, 0._DP)
      T(n,n-1) = (-1._DP, 0._DP)
      x = 1._DP-h
      T(n,n) = DCMPLX(2._DP + f(x)*h2, -alpha)

      do i=2,n-1
         x = dble(i)*h
         T(i,i-1) = (-1._DP, 0._DP)
         T(i,i) = DCMPLX(2._DP + f(x)*h2, -alpha)
         T(i,i+1) = (-1._DP, 0._DP)
      enddo

      Tstar = CONJG(T)

      bsum=0._DP
      do i=1,n
        x=dble(i)*h
        B(i) = dexp(-(x-0.25_DP)**2/0.002) 
        x= 39.d0*pi*(x-0.5d0)
        B(i) = B(i)*DCMPLX( dcos(x), dsin(x) )
        bsum=bsum+ zabs(B(i))**2
      enddo
      bsum = 1._DP/sqrt(h*bsum)
      do i=1,n
        B(i) = B(i) * bsum
      enddo

!      write(filename,'("frame",i6.6,".dat")')0
      open(20,file='initial_frame.dat',status='unknown')
      write(20,*)"# time step ",0
      write(20,*)"# position   real  imag   prob"
      write(20,*)0,    0,       0,      0
      do i=1,n
        x=real(i)*h
        write(20,*)x, REALPART(B(i)), IMAGPART(B(i)), ZABS( B(i)*CONJG(B(i)) )
      enddo
      write(20,*)1,    0,       0,      0
      close(20)

!      open(20,file='movie.dat',status='unknown')
!      write(20,*)"# time step ",0
!      write(20,*)"# position   real  imag   prob"
!      write(20,*)0,    0,       0,      0
!      do i=1,n
!        x=real(i)*h
!        write(20,*)x, REALPART(B(i)), IMAGPART(B(i)), ZABS( B(i)*CONJG(B(i)) )
!      enddo
!      write(20,*)1,    0,       0,      0
!      write(20,*)


! Now propagate foward in time!
      open(unit=1,file="XvsTime.dat",status='unknown')

      do j=1,2000
        B= MATMUL(Tstar,-B)
        call ZGESV(n,1,T,n,IPIV,B,n,INFO)

!       B is now one time step ahead
!       need to define function TRAPEZOID to find average x
        write(1,*) dble(j)*dt, trapezoid(B,n,h)

        T = CONJG(Tstar)

!        if(mod(j,20).eq.0)then
!           write(20,*)
!           write(20,*)"# time step ",j
!           write(20,*)"# position   real  imag   prob"
!           write(20,*)0,    0,       0,      0
!           do i=1,n
!              x=real(i)*h
!              write(20,*)x, REALPART(B(i)), IMAGPART(B(i)), ZABS( B(i)*CONJG(B(i)) )
!           enddo
!           write(20,*)1,    0,       0,      0
!           write(20,*)
!        endif
      enddo


      close(20)
      end program sqwelltime

      FUNCTION f(x)
      implicit none
      real*8 , intent(in) :: x
      real*8 :: f

      f = 0.

      END FUNCTION f


      FUNCTION trapezoid(B,n,dx)
      implicit none
      complex*16 , dimension(:), intent(in) :: B
      integer, intent(in) :: n
      real*8, intent(in) :: dx
      integer :: i
      real*8 :: trapezoid, x

      trapezoid = 0.d0
      do i = 1, n
        x = dble(i)*dx
	trapezoid = 0.d0
      enddo
      trapezoid=trapezoid*dx

      END FUNCTION trapezoid