c program for calculating a 1-dimensional electron wave
c atomic units used throughout
c everything on an equally spaced grid x_j=x_0+dx*j
c
      program onedschreq
      implicit real*8 (a-h,o-z)
      parameter(nptfin=10000)
c vpot will contain the potential + diagonal term from kinetic energy
      dimension vpot(nptfin)
c
c psi will contain the wave function
c phi will contain the inhomogeneous term
c diag will contain the diagonal term of the matrix
c offd will contain the offdiagonal term of the matrix, since
c      this does not depend on the index for this problem it
c      has been changed to a number
c the matrix is assumed to be symmetric
c
      complex*16 psi(0:nptfin),phi(0:nptfin),diag(0:nptfin),
     1  offd

      open(unit=25,file='onedtimsch1.dat',status='old')
c
c read in the data for a specific run
c
      read(25,*)ak,x0,xf,xm,wi,npts,dt,numtime,nmod
c
c ak = initial momentum ; x0 = lowest value of x
c xf = largest value of x ; xm = middle of packet
c wi = width of initial wave function
c npts = number of spatial points ; dt = time step size
c numtime = number of time steps to do propagation
c nmod = how many steps to skip before printing out wave function
c
      if(npts .ge. nptfin) stop 'increase the size of nptfin'
c
c initialize the wave function
c
      dx=(xf-x0)/dfloat(npts+1)
      tfin=dt*dfloat(numtime)
      iout=numtime/nmod
      write(6,*)'dx=',dx,' tfin= ',tfin,
     1  ' number of output wave functions = ',iout+1
       dx2m1=1.d0/(dx*dx)
      sum=0.d0
      psi(0)=(0.,0.)
      psi(npts+1)=(0.,0.)
       do 1 j = 1,npts
       x=x0+dx*dfloat(j)
       s=(x-xm)/wi
       psi(j)=exp(-s*s+(0.,1.)*ak*x)
c
c if the potential is time dependent need to put this in the 4 loop
c
        potential = -x/1.d2
        vpot(j)=potential+dx2m1
c
       sum = sum + abs(psi(j))**2
 1     continue
c
c normalize the wave function
c
      xnorm=1.d0/dsqrt(sum)
       do 2 j = 0,npts+1
       psi(j)=psi(j)*xnorm
 2     continue

c
c output the initial wave function
c
       do 202 j = 0,npts+1,10
       x=x0+dx*dfloat(j)
       write(100,900)x,abs(psi(j))**2
 202   continue
       close(100)

c
c the 3 loop does the time propagation
c
       pref=dt*dx2m1/4.d0
c
c set the off diagonal term
c
       offd=(0.,-1.)*pref
c
       do 3 itime = 1,numtime
       t=dt*dfloat(itime)
c
c calculate phi and initialize the diagonal and off diagonal elements
c
        do 4 j = 1,npts
        x=x0+dx*dfloat(j)
        phi(j)=(1.-(0.,.5)*dt*vpot(j))*psi(j)+(0.,1.)*pref*
     1   (psi(j+1)+psi(j-1))
        diag(j)=(1.+(0.,.5)*dt*vpot(j))
 4      continue
c
c sweep down the index to do the LU decomposition
c
        do 5 j = 2,npts
        phi(j)=phi(j)-phi(j-1)*offd/diag(j-1)
        diag(j)=diag(j)-offd*offd/diag(j-1)
 5      continue
c
c sweep up the index to do the back substitution
c
        do 6 j = npts,1,-1
        psi(j)=(phi(j)-offd*psi(j+1))/diag(j)
 6      continue
c
c at the appropriate times output various parts of the wave function
c
       ites=mod(itime,nmod)
        if(ites .eq. 0) then
        igo=itime/nmod
         do 7 j = 0,npts+1,10
         x=x0+dfloat(j)*dx
         tes=abs(psi(j))**2
         write(100+igo,900)x,tes
 7       continue
        close(100+igo)
        end if
 3     continue
 900   format(2f20.10)
       stop
       end