************************************************************************
       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,M,N,KK,LDA,LDB,LDC
       REAL*8 DX, DT, PI, TWOPI, A_HALF, WF_norm, X0, XD
       REAL*8 DX_3, X_start, X_av, PSI_norm, W0, W1, U0, U1
       REAL*8 Kdt_r, Kprop_r, Kw_r, PSI_r, PSI1_r, C0, SUM0, SUM1
       REAL*8 Kdt_i, Kprop_i, Kw_i, PSI_i, PSI1_i, Ktemp_r, Ktemp_i
       REAL*8 ALPHA, BETA
       PARAMETER (M=1001,N=1001,KK=1001)
       DIMENSION Kdt_r(M,M), Kprop_r(M,M), Kw_r(M,M)
       DIMENSION Kdt_i(M,M), Kprop_i(M,M), Kw_i(M,M)
       DIMENSION Ktemp_r(M,M), Ktemp_i(M,M)
       DIMENSION PSI_r(M), PSI_i(M), PSI1_r(M), PSI1_i(M)
       DATA PI/3.141592653589793D0/
       
       LDA = M
       LDB = M
       LDC = M
       
       ND       =  1000
       X0       = -4.0D0
       XD       =  4.0D0
       X_start  =  0.75D0
       A_HALF   =  1.0D0
       DX       =  (XD-X0)/DFLOAT(ND)
       WF_norm  =  (2.0D0*A_HALF/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) )
       
       ALPHA = DX
       BETA  = 0.0D0
 
       DO K = 0,ND 
        PSI_r(K+1) = WF_norm*DEXP(-A_HALF*(X0+DX*DFLOAT(K)-X_start)**2)	 
	PSI_i(K+1) = 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	 
                 
         Kdt_r(K+1,L+1) = C0*(DSIN(W0) + DCOS(W0))
         Kdt_i(K+1,L+1) = C0*(DSIN(W0) - DCOS(W0))
         Kw_r(K+1,L+1)  = Kdt_r(K+1,L+1)
         Kw_i(K+1,L+1)  = Kdt_i(K+1,L+1)	 
       ENDDO
       ENDDO 
       
       X_av     = 0.0D0 
       PSI_norm = 0.0D0
            
       DO K=0,ND
          PSI_norm = PSI_norm + PSI_r(K+1)*PSI_r(K+1) + 
     >               PSI_i(K+1)*PSI_i(K+1)
          X_av     = X_av + (X0+DFLOAT(K)*DX)*( PSI_r(K+1)*PSI_r(K+1) + 
     >                                          PSI_i(K+1)*PSI_i(K+1) )	
       ENDDO 
          PSI_norm = DX*PSI_norm
	  X_av     = DX*X_av
       
       print *, PSI_norm, X_av
              
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kw_r,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kw_i,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kprop_r(K,L) = Ktemp_r(K,L) - Ktemp_i(K,L) 
       ENDDO
       ENDDO	 


       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kw_i,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kw_r,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kprop_i(K,L) = Ktemp_r(K,L) + Ktemp_i(K,L) 
       ENDDO
       ENDDO	


       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kprop_r,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kprop_i,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kw_r(K,L) = Ktemp_r(K,L) - Ktemp_i(K,L) 
       ENDDO
       ENDDO	 


       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kprop_i,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kprop_r,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kw_i(K,L) = Ktemp_r(K,L) + Ktemp_i(K,L) 
       ENDDO
       ENDDO


       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kw_r,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kw_i,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kprop_r(K,L) = Ktemp_r(K,L) - Ktemp_i(K,L) 
       ENDDO
       ENDDO	 


       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_r,LDA,Kw_i,LDB,
     >                                           BETA,Ktemp_r,LDC)       
       CALL DGEMM ('N','N',M,N,KK,ALPHA,Kdt_i,LDA,Kw_r,LDB,
     >                                           BETA,Ktemp_i,LDC) 
     
       DO K = 1, M
       DO L = 1, M       
         Kprop_i(K,L) = Ktemp_r(K,L) + Ktemp_i(K,L) 
       ENDDO
       ENDDO

c      stop
       
c       DO K=0,ND
c
c       SUM0 = 0.0D0
c       SUM1 = 0.0D0
c       
c       DO J=1,M
c       DO I=1,M
c       
c        W0 =  Kprop_r(I,J)
c	W1 = -Kprop_i(I,J)
c       U0 =  Kprop_r(I,K)
c       U1 =  Kprop_i(I,K)
c	
c	SUM0 = SUM0 + DX*DX*(W0*U0 - W1*U1)
c	SUM1 = SUM1 + DX*DX*(W0*U1 + W1*U0)
c
c       ENDDO
c       ENDDO
c      
c       print *, K, SUM0, SUM1
c       
c       ENDDO          

       
       DO J=1,128
       
       DO K=1,M
          PSI1_r(K) = 0.0D0
	  PSI1_i(K) = 0.0D0
       DO L=1,M
          W0 = Kprop_r(K,L)
	  W1 = Kprop_i(K,L)
	  U0 = PSI_r(L)
	  U1 = PSI_i(L)                   
          PSI1_r(K) = PSI1_r(K) + DX*(W0*U0 - W1*U1)
	  PSI1_i(K) = PSI1_i(K) + DX*(W0*U1 + W1*U0)
       ENDDO
       ENDDO
       
       DO K=1,M
         PSI_r(K) = PSI1_r(K)
	 PSI_i(K) = PSI1_i(K)
       ENDDO      
       
       X_av     = 0.0D0 
       PSI_norm = 0.0D0
            
       DO K=1,M
          PSI_norm = PSI_norm + PSI_r(K)*PSI_r(K) + PSI_i(K)*PSI_i(K)
          X_av     = X_av + (X0+DFLOAT(K-1)*DX)*( PSI_r(K)*PSI_r(K) + 
     >                                            PSI_i(K)*PSI_i(K) )	
       ENDDO 
       PSI_norm = DX*PSI_norm
       X_av     = DX*X_av
       
       print *, J, PSI_norm, X_av 
       
       ENDDO 

       STOP
       END