************************************************************************
       PROGRAM PATH 
************************************************************************
*      Calculates the time evolution of wavefunction in real time
*      using the Feyman Path Integral Method.
*      Integration is approximated by Riemann summation
*      Time evolution is converted into repeated matrix multiplications  
*      of the small time step propagator
************************************************************************
*    Variables:
*
*    K_dt(N,N) short time propagator matrix with N discretization points
*           
*
*
************************************************************************
       IMPLICIT NONE
       
       INTEGER I,J,JP,K,L,ND
       REAL*8 DX, DT, PI, TWOPI, ALPHA, A_HALF, WF_norm, X0, XD
       REAL*8 DX_3, X_start, X_av, PSI_norm, W0, W1, U0, U1
       REAL*8 K_dt, K_prop, K_w, PSI, PSI_1, C0, SUM0, SUM1
       DIMENSION K_dt(0:1,0:1000,0:1000), K_prop(0:1,0:1000,0:1000)
       DIMENSION K_w(0:1,0:1000,0:1000)
       DIMENSION PSI(0:1,0:1000), PSI_1(0:1,0:1000)
       DATA PI/3.141592653589793D0/
       
       ND       =  600
       X0       = -4.0D0
       XD       =  4.0D0
       ALPHA    =  2.0D0
       X_start  =  0.75D0
       A_HALF   =  ALPHA/2.0D0
       DX       =  (XD-X0)/DFLOAT(ND)
       WF_norm  =  (ALPHA/PI)**0.25D0
       DX_3     =  DX*DX*DX
      
       TWOPI =  2.0D0*PI 
       DT    =  TWOPI/128.0D0 
       
       C0    =  1.0D0/DSQRT( 4.0D0*PI*DSIN(DT) )
 
       DO K = 0,ND 
         PSI(0,K) = WF_norm*DEXP(-A_HALF*(X0+DX*DFLOAT(K)-X_start)**2)	 
	 PSI(1,K) = 0.0D0
       ENDDO	 

       DO K = 0,ND
       DO L = 0,ND  
         W0 = ( ( (X0+DFLOAT(K)*DX)*(X0+DFLOAT(K)*DX) + 
     >    (X0+DX*DFLOAT(L))*(X0+DX*DFLOAT(L)) ) -
     >    2.0D0*(X0+DX*DFLOAT(K))*(X0+DX*DFLOAT(L)) )/(2.0D0*DT) -
     >    (  (X0+DFLOAT(K)*DX)*(X0+DFLOAT(K)*DX) + 
     >    (X0+DX*DFLOAT(L))*(X0+DX*DFLOAT(L))  +
     >    2.0D0*(X0+DX*DFLOAT(K))*(X0+DX*DFLOAT(L)) )*DT*0.125D0	 
                 
         K_dt(0,K,L) = C0*(DSIN(W0) + DCOS(W0))
         K_dt(1,K,L) = C0*(DSIN(W0) - DCOS(W0))
       ENDDO
       ENDDO 
       
       X_av     = 0.0D0 
       PSI_norm = 0.0D0
            
       DO K=0,ND
          PSI_norm = PSI_norm + PSI(0,K)*PSI(0,K) + PSI(1,K)*PSI(1,K)
          X_av     = X_av + (X0+DFLOAT(K)*DX)*( PSI(0,K)*PSI(0,K) + 
     >                                          PSI(1,K)*PSI(1,K) )	
       ENDDO 
          PSI_norm = DX*PSI_norm
	  X_av     = DX*X_av
       
       print *, PSI_norm, X_av
       

       DO K  = 0,ND
        DO L  = 0,ND
	
	 K_w(0,K,L) = 0.0D0
	 K_w(1,K,L) = 0.0D0
	 
         DO I  = 0,ND 
            W0 =  K_dt(0,K,I)
	    W1 =  K_dt(1,K,I)
            U0 =  K_dt(0,I,L)
	    U1 =  K_dt(1,I,L)	 	         
            K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
	    K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
         ENDDO 	
	 
        ENDDO
       ENDDO 
              
       DO K=0,ND
       DO L=0,ND
         K_prop(0,K,L) = DX*K_w(0,K,L)
	 K_prop(1,K,L) = DX*K_w(1,K,L)
       ENDDO
       ENDDO           
       

       DO K  = 0,ND
        DO L  = 0,ND
	
	 K_w(0,K,L) = 0.0D0
	 K_w(1,K,L) = 0.0D0
	 
         DO I  = 0,ND 
            W0 =  K_dt(0,K,I)
	    W1 =  K_dt(1,K,I)
            U0 =  K_prop(0,I,L)
	    U1 =  K_prop(1,I,L)	 	         
            K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
	    K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
         ENDDO 	
	 
        ENDDO
       ENDDO 
              
       DO K=0,ND
       DO L=0,ND
         K_prop(0,K,L) = DX*K_w(0,K,L)
	 K_prop(1,K,L) = DX*K_w(1,K,L)
       ENDDO
       ENDDO


       DO K  = 0,ND
        DO L  = 0,ND
	
	 K_w(0,K,L) = 0.0D0
	 K_w(1,K,L) = 0.0D0
	 
         DO I  = 0,ND 
            W0 =  K_dt(0,K,I)
	    W1 =  K_dt(1,K,I)
            U0 =  K_prop(0,I,L)
	    U1 =  K_prop(1,I,L)	 	         
            K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
	    K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
         ENDDO 	
	 
        ENDDO
       ENDDO 
              
       DO K=0,ND
       DO L=0,ND
         K_prop(0,K,L) = DX*K_w(0,K,L)
	 K_prop(1,K,L) = DX*K_w(1,K,L)
       ENDDO
       ENDDO
       
       DO K=0,ND

       SUM0 = 0.0D0
       SUM1 = 0.0D0
       
       DO J=0,ND
       DO I=0,ND
       
        W0 =  K_prop(0,I,J)
	W1 = -K_prop(1,I,J)
        U0 =  K_prop(0,I,K)
	U1 =  K_prop(1,I,K)
	
	SUM0 = SUM0 + DX*DX*(W0*U0 - W1*U1)
	SUM1 = SUM1 + DX*DX*(W0*U1 + W1*U0)

       ENDDO
       ENDDO
       
       print *, K, SUM0, SUM1
       
       ENDDO          

       
       DO J=1,128
       
       DO K=0,ND
          PSI_1(0,K) = 0.0D0
	  PSI_1(1,K) = 0.0D0
       DO L=0,ND 
          W0 = K_prop(0,K,L)
	  W1 = K_prop(1,K,L)
	  U0 = PSI(0,L)
	  U1 = PSI(1,L)                   
          PSI_1(0,K) = PSI_1(0,K) + DX*(W0*U0 - W1*U1)
	  PSI_1(1,K) = PSI_1(1,K) + DX*(W0*U1 + W1*U0)
       ENDDO
       ENDDO
       
       DO K=0,ND
         PSI(0,K) = PSI_1(0,K)
	 PSI(1,K) = PSI_1(1,K)
c	 print *, K, PSI(0,K), PSI(1,K)
       ENDDO      
       
       X_av     = 0.0D0 
       PSI_norm = 0.0D0
            
       DO K=0,ND
          PSI_norm = PSI_norm + PSI(0,K)*PSI(0,K) + PSI(1,K)*PSI(1,K)
          X_av     = X_av + (X0+DFLOAT(K)*DX)*( PSI(0,K)*PSI(0,K) + 
     >                                          PSI(1,K)*PSI(1,K) )	
       ENDDO 
       PSI_norm = DX*PSI_norm
       X_av     = DX*X_av
       
       
       print *, J, PSI_norm, X_av 
       
       ENDDO 

       STOP
       END