      PROGRAM tdse
      IMPLICIT NONE
      REAL Im(0:1100), Imold(0:1100), Re(0:1100), P(0:1100), 
     +       V0, a, dt, dx, dx2, k0, t, width, x0, xmax, xmin
      INTEGER IC1_, IC2_, IC3_, IC4_, IR_, i, n
      DATA Im, Imold, Re, P/4404*0.0/
      LOGICAL TBkeyinput
      DATA IC1_/11/
      DATA IC2_/12/
      DATA IC3_/13/
      DATA IC4_/14/
      CALL GWopen(IR_, 0)
      CALL GWmode(IR_, 3, 0)
*     motion of a wavepacket incident on a potential
      CALL parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt)
      CALL initial_packet(Re,Im,x0,k0,width,xmin,n,dx,dt)
      CALL set_up_windows(xmin,xmax,V0,IC1_,IC2_,IC3_,IC4_)
      CALL draw_potential(V0,a,xmin,xmax,IC4_)
      CALL GWanchor(IR_, 1)
      DO WHILE(.NOT.TBkeyinput())
         DO i = 1, 100
            CALL evolve(Re,Im,Imold,t,V0,a,dx,dx2,dt,xmin,n)
         ENDDO
         CALL plots(Re,Im,Imold,P,xmin,n,dx,IC1_,IC2_,IC3_,IC4_)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt)
      IMPLICIT NONE
      REAL x0, k0, width, V0, a, xmin, xmax, dx, dx2, dt
      INTEGER n
      x0 = -15                              ! initial position of packet
      width = 1                             ! width of wave packet in x space
      k0 = 2                                ! group velocity of packet
      xmax = 20                             ! maximum value of x
      xmin = -xmax                          ! minimum value of x
      V0 = 2
      a = 1                                 ! half-width of potential well
      dx = 0.1                              ! grid size
      dx2 = dx*dx
      n = (xmax - xmin)/dx
      dt = 0.005                            ! time interval
      END
C
      SUBROUTINE initial_packet(Re,Im,x0,k0,width,xmin,n,dx,dt)
      IMPLICIT NONE
      REAL A, arg, arg2, b, delta2, e, e2, pi, x, Re(0:1100),
     +       Im(0:1100), x0, k0, width, xmin, dx, dt
      PARAMETER(pi = 3.141592)
      INTEGER i, n
*     initial Gaussian wavepacket
      delta2 = width*width
      A = (2*pi*delta2)**(-0.25)
      b = 0.5*k0*dt
      DO i = 1, n
         x = xmin + (i-1)*dx
         arg = 0.25*(x - x0)**2/delta2
         arg2 = 0.25*(x - x0 - 0.5*k0*dt)**2/delta2
         e = exp(-arg)
         e2 = exp(-arg2)
         Re(i) = A*cos(k0*(x - x0))*e       ! value at t = 0
         Im(i) = A*sin(k0*(x - x0 - 0.5*b))*e2 ! value at t = dt/2
      END DO
      END
C
      SUBROUTINE set_up_windows(xmin,xmax,V0,IC1_,IC2_,IC3_,IC4_)
      IMPLICIT NONE
      REAL GWaspect, xmin, xmax, V0
      INTEGER IR_, IC1_, IC2_, IC3_, IC4_
      CALL GWvport(IR_, 0.0, 0.75, GWaspect(1), 1.0)    ! probability
      CALL GWindow(IR_, xmin, -0.1, xmax, 0.5)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWline(IR_, xmin, 0.0, xmax, 0.0)
      CALL GWsetpen(IR_, 16, -1, -1, -1);
      CALL GWsavevp(IR_, IC1_)
      CALL GWvport(IR_, 0.0, 0.5, GWaspect(1), 0.75)    ! real part
      CALL GWindow(IR_, xmin, -1.0, xmax, 1.0)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWline(IR_, xmin, 0.0, xmax, 0.0)
      CALL GWsetpen(IR_, 16, -1, -1, -1);
      CALL GWsavevp(IR_, IC2_)
      CALL GWvport(IR_, 0.0, 0.25, GWaspect(1), 0.5)    ! imaginary part
      CALL GWindow(IR_, xmin, -1.0, xmax, 1.0)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWline(IR_, xmin, 0.0, xmax, 0.0)
      CALL GWsetpen(IR_, 16, -1, -1, -1);
      CALL GWsavevp(IR_, IC3_)
      CALL GWvport(IR_, 0.0, 0.0, GWaspect(1), 0.22)    ! potential
      CALL GWindow(IR_, xmin, -V0, xmax, V0)
      CALL GWsavevp(IR_, IC4_)
      END
C
      SUBROUTINE draw_potential(V0,a,xmin,xmax,IC4_)
      IMPLICIT NONE
      REAL V0, a, xmin, xmax
      INTEGER IR_, IC4_
      CALL GWselvp(IR_, IC4_)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWmove2(IR_, xmin, 0.0)
      CALL GWline2(IR_, a, 0.0)
      CALL GWline2(IR_, a, V0)
      CALL GWline2(IR_, xmax, V0)
      CALL GWindow(IR_, 0.0, 20.0, 80.0, 0.0)
      CALL GWsavevp(IR_, IC4_)
      END
C
      SUBROUTINE evolve(Re,Im,Imold,t,V0,a,dx,dx2,dt,xmin,n)
      IMPLICIT NONE
      REAL HIm, HRe, x, Re(0:1100), Im(0:1100), Imold(0:1100), t, V0,
     +       a, dx, dx2, dt, xmin, V
      INTEGER i, n
      DO i = 1, n
         x = xmin + (i-1)*dx
         HIm = V(x,V0,a)*Im(i) -0.5*(Im(i+1) -2*Im(i) +Im(i-1))/dx2
*        real part defined at multiples of dt
         Re(i) = Re(i) + HIm*dt
      END DO
      DO i = 1, n
         x = xmin + (i-1)*dx
*        dt/2 earlier than real part
         Imold(i) = Im(i)
         HRe = V(x,V0,a)*Re(i) -0.5*(Re(i+1) -2*Re(i) +Re(i-1))/dx2
*        dt/2 later than real part
         Im(i) = Im(i) - HRe*dt
      END DO
      t = t + dt                            ! time of real part
      END
C
      SUBROUTINE plots(Re,Im,Imold,P,xmin,n,dx,IC1_,IC2_,IC3_,IC4_)
      IMPLICIT NONE
      CHARACTER Stmp1_*80
      REAL Psum, x, Re(0:1100), Im(0:1100), Imold(0:1100), P(0:1100),
     +       xmin, dx
      INTEGER IR_, i, n, IC1_, IC2_, IC3_, IC4_
      DATA Psum/0.0/
      CALL GWselvp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, -IC2_)
      CALL GWplot1(IR_, -1, n, xmin, xmin+n*dx, 0.0,0.0,1.0,0.0, Re(1))
      CALL GWflush(IR_, -IC2_)
      CALL GWselvp(IR_, IC3_)
      CALL GWsetogn(IR_, 0, -IC3_)
      CALL GWplot1(IR_, -1, n, xmin, xmin+n*dx, dx/2, 0.0,1.0,0.0,Im(1))
      CALL GWflush(IR_, -IC3_)
      Psum = 0.0
      DO i = 1, n
         x = xmin + (i - 1)*dx
         P(i) = Re(i)*Re(i) + Im(i)*Imold(i)
         Psum = Psum + P(i)*dx
      END DO
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, -IC1_)
      CALL GWplot1(IR_, -1, n, xmin, xmin+n*dx, dx/2, 0.0,1.0,0.0, P(1))
      CALL GWflush(IR_, -IC1_)
      CALL GWselvp(IR_, IC4_)
      CALL GWsetogn(IR_, 0, -IC4_)
      WRITE(Stmp1_,*) 'total probability =', Psum
      CALL GWputtxt(IR_, 1.0, 1.0, Stmp1_)  ! total probability should be unity
      CALL GWflush(IR_, -IC4_)
      CALL GWsetogn(IR_, 0, 1)
      END
C
      FUNCTION V(x,V0,a)
      IMPLICIT NONE
      REAL V, x, V0, a
*     step potential
      IF(x .GT. a) THEN
         V = V0
      ELSE
         V = 0
      END IF
      END
