program     tdse_fd
  implicit none
  real(8), allocatable :: repsi(:,:), impsi(:,:)
  real(8), allocatable :: V(:), rho(:)
  integer, parameter   :: new = 2, old = 1
  integer              :: tskip
  real(8)              :: k0, sigma, sigma2
  real(8)              :: xmax, tmax, x, t, dt, dx, a, amax, x0, norm, pexp, delpsi
  integer              :: Nx, Nt
  integer              :: ix, it
  integer              :: chk
  real(8), parameter   :: tinypsi = 1.0D-20
  real(8), parameter   :: PI      = atan2(1.0D0,0.0D0) * 2.0D0
  integer, parameter   :: f_psi   = 11 , f_rho = 12
  !-------------------------------------------------------------------------------
  print *, '# Enter x0, sigma, k0, dx, a, xmax, tmax, tskip:'
  read  *,          x0, sigma, k0, dx, a, xmax, tmax, tskip
  !-------------------------------------------------------------------------------
  Nx = INT(2.0D0*xmax/dx) + 1
  dx =     2.0D0*xmax/(Nx-1)
  !-------------------------------------------------------------------------------
  allocate(repsi(2,Nx), impsi(2,Nx), stat=chk)
  if( chk /= 0 ) stop 'Allocation error (1)'
  allocate(V      (Nx), rho    (Nx), stat=chk)
  if( chk /= 0 ) stop 'Allocation error (2)'
  !-------------------------------------------------------------------------------
  !Potential:
  V        = 0.0D0
  do ix    = 1, Nx
     x     = -xmax + (ix-1) * dx
     V(ix) = 0.5D0 * x * x
  end do
  !-------------------------------------------------------------------------------
  amax     = 1.0D0/(2.0D0 + 0.5D0 * xmax * xmax * dx * dx)
  !-------------------------------------------------------------------------------
  dt       = 2.0D0 * a * dx * dx
  Nt       = INT(tmax/dt) + 1
  dt       =     tmax/(Nt-1)
  a        = 0.5D0 * dt / (dx*dx)
  sigma2   = sigma * sigma
  !-------------------------------------------------------------------------------
  print *, '#----------------------------------------------------------'
  print *, '# x0, sigma, k0'
  print *,    x0, sigma, k0
  print *, '# xmax, tmax, tskip:'
  print *,    xmax, tmax, tskip
  print *, '# Nx, Nt, dx, dt:'
  print *,    Nx, Nt, dx, dt
  print *, '# a, amax:'
  print *,    a, amax
  if(a>amax) print *, '# WARNING: a>amax'
  print *, '#----------------------------------------------------------'
  !-------------------------------------------------------------------------------
  open(unit=f_rho, file='rho')
  open(unit=f_psi, file='psi')
  !-------------------------------------------------------------------------------
  !Initial condition:
  t             = 0.0D0
  rho           = 0.0D0
  repsi(old,1)  = 0.0D0  ;  impsi(old,1 ) = 0.0D0
  repsi(old,Nx) = 0.0D0  ;  impsi(old,Nx) = 0.0D0
  norm          = (2.0D0 * PI * sigma2) ** (-0.25D0)
  do ix   = 2, Nx - 1
     x    = -xmax + (ix-1)*dx
     pexp = 0.25D0 * (x-x0) * (x-x0) / sigma2
     if( pexp < -log(tinypsi) )then
        repsi(old,ix) = norm * exp( -pexp) * cos( k0 * (x-x0) )
        impsi(old,ix) = norm * exp( -pexp) * sin( k0 * (x-x0) )
     else
        repsi(old,ix) = 0.0D0
        impsi(old,ix) = 0.0D0
     end if
     rho     (    ix) = repsi(old,ix) ** 2 + impsi(old,ix)**2
  end do
  norm                = SUM(rho) * dx  ! = int rho(x) dx
  !-------------------------------------------------------------------------------
  !Output initial state:
  it = 0 ; t = 0.0D0
  do ix = 1, Nx
     x  = -xmax + (ix-1) * dx
     write(f_psi, 101) it, t, x, rho(ix), repsi(old,ix), impsi(old,ix)
  end do
  write(f_psi,*) ''
  !-------------------------------------------------------------------------------
  !Half step for impsi:
  do ix = 2, Nx-1
     delpsi = -a*(repsi(old,ix+1)-2.0D0*repsi(old,ix)+repsi(old,ix-1)) -V(ix) * dt * repsi(old,ix)
     delpsi = 0.5D0 * delpsi
     if(abs(delpsi) > tinypsi) impsi(old,ix) = impsi(old,ix) + delpsi
  end do
  !-------------------------------------------------------------------------------
  !Time evolution
  do it = 1, Nt
     t  = it*dt

     !update repsi:
     do ix     = 2, Nx - 1
        delpsi = -a*(impsi(old,ix+1)-2.0D0*impsi(old,ix)+impsi(old,ix-1))+V(ix)*dt*impsi(old,ix)
        if(abs(delpsi) > tinypsi) repsi(new,ix) = repsi(old,ix) + delpsi
        rho(ix)=repsi(old,ix)*repsi(new,ix) + impsi(old,ix)*impsi(old,ix)
     end do

     !update impsi:
     do ix     = 2, Nx - 1
        delpsi =  a*(repsi(new,ix+1)-2.0D0*repsi(new,ix)+repsi(new,ix-1))-V(ix)*dt*repsi(new,ix)
        if(abs(delpsi) > tinypsi) impsi(new,ix) = impsi(old,ix) + delpsi
     end do

     !print int rho(x) dx:
     norm      = sqrt(SUM(rho) * dx)
     write(f_rho,*) it, t, norm

     !print psi:
     if(MOD(it,tskip) == 0)then
        do ix = 1, Nx
           x  = -xmax + (ix-1)*dx
           write(f_psi,101) it, t, x, rho(ix), repsi(new,ix), impsi(new,ix)
        end do
        write(f_psi,*)''
     end if

     !interchange new <-> old:
     repsi(old,:) = repsi(new,:)
     impsi(old,:) = impsi(new,:)
     
  end do    ! do it = 1, Nt

  close(f_rho)
  close(f_psi)

101 FORMAT(I20,1000E30.15)
end program tdse_fd
