PROGRAM TDSE 
! Copyright 2007 Matthew Norton
!
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
!
! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run this code
!---------------------------------------------------------------------------
! This program models particles falling in an electric field using quantum
! mechanics. This program outputs the position distribution at t = 0, t = 2,
! and t = 4. Thus producing the movement of the particles as time increases.
! The program also outputs the mean position and standard deviation for each time
! This program advances the wave function in time using the Fourier transform method.
! 
! Programmed by:
!     Matthew Norton
!
!----------------------------------------------------------------------------

IMPLICIT NONE
 REAL*8, PARAMETER :: E0 = 1, Dx0 = 1.0, dX = 0.1, Pie = -1.0, TMax = 51.20
 INTEGER, PARAMETER :: NDIMX = 512
 INTEGER :: I, SIGN, FileUnit = 20
 REAL*8 :: Pi, dw, Xn(NDIMX), PhiK(NDIMX), Psum = 0.0, P2 = 0.0, PStd = 0.0, K(NDIMX)
 REAL*8 :: XSum = 0.0, X2 = 0.0, Ps(NDIMX), Time
 COMPLEX*8 :: PPsi(NDIMX)=(0,0), Psi(NDIMX)=(0,0) , Mom, Sum = (0.0,0.0)
 COMPLEX*8 :: PP(NDIMX)=(0,0)
COMPLEX*8 :: PPX(NDIMX)=(0,0)
 COMPLEX*8, PARAMETER :: Imag = (0,1)
 REAL      T
 INTEGER   ITIME,Q,L
Pi = ACOS(Pie)
DW = (2. * PI) / TMax
 OPEN(UNIT = 20, FILE = "Psi(xt=0).txt")
 OPEN(UNIT = 25, FILE = "Psi(xt=2).txt")
 OPEN(UNIT = 30, FILE = "Psi(xt=4).txt")
 DO I = 1, NDIMX
      Xn(I) = (- NDIMX / 2.0 - 0.5 + I) * dX
      Psi(I) =  DEXP(-(Xn(I)/Dx0) ** 2 / 4.0)
 END DO
        call norm(psi)
!
!           DO L=1,512
!           PPX(L)=PSI(L)
!           ENDDO
 DO I = NDIMX / 2 + 1, NDIMX
      K(I) = (-NDIMX - 1 + I) * Dw
 END DO
 DO I = 1, NDIMX/2
      k(I) = (I-1) * Dw
 END DO
       xsum=0.
       x2=0.
 DO I = 1, NDIMX
      XSum = XSum + Xn(I) *      REAL(Psi(I) * CONJG(Psi(I)))
      X2   = X2   + Xn(I) ** 2 * REAL(Psi(I) * CONJG(Psi(I)))
 END DO
 PRINT *, "XSum", XSum, "X^2", X2, "STDEV", SQRT(X2 -XSum**2)
 SIGN = -1
 CALL FFTBYHAND(SIGN, Psi, PPsi)
            call norm(PPsi)

           DO L=1,512
           PP(L)=PPSI(L)
           ENDDO

 DO ITIME=0,100,50
         T=ITIME*0.04
DO I = 1, NDIMX
      Mom = (0.,1.)/6.*(  (K(I) - T)**3 - K(I)**3  )
      PPsi(I) = PP(I) * EXP(Mom)
 END DO

 SIGN = 1
 CALL FFTBYHAND(SIGN, PPsi, Psi)
 CALL norm(psi) 
DO I = 1, NDIMX
     WRITE(FileUnit,*) Xn(I), REAL(Psi(I) * CONJG(Psi(I)))
END DO
 FileUnit = FileUnit + 5  
 XSUM=0.
 X2=0.
 DO Q = 1, NDIMX
      XSum = XSum + Xn(Q) * ABS(PSI(Q))**2
      X2 = X2 + Xn(Q) ** 2 * ABS(PSI(Q))**2
     END DO
     PRINT *, "XSum", XSum, "X^2", X2, "STDEV", SQRT(ABS(X2 -(XSum * XSum)))
!      WRITE(20,*) T, XSUM,SQRT(ABS(X2 -(XSum * XSum)))
END DO

 ! -------------------------------------------------------------------------------------------------------
 CONTAINS
! -------------------------------------------------------------------------------------------------------
  ! Subroutine for slow fourier transform
  SUBROUTINE FFTBYHAND(SIGN, HIN, HOUT)
  COMPLEX*8 E1, SUM, HIN(NDIMX), HOUT(NDIMX), TPII
  INTEGER I, I1, SIGN
  TPII = 2. * Pi * (0.0,1.0)! 2PI
  DO I = 1, NDIMX
       SUM = (0.0,0.0)
       DO I1 = 1, NDIMX
            E1 = EXP(SIGN * TPII * (I1 - 1) * (I - 1.) / (NDIMX + 0.0))
            SUM = SUM + E1 * HIN(I1)
       END DO
       HOUT(I) = SUM
  END DO
  END SUBROUTINE FFTBYHAND
! --------------------------------------------------------------------------------------------------------
SUBROUTINE norm(hin)
  complex*8 hin(512),sum
  INTEGER L
  sum=0.
  DO l=1,512
       sum=sum+abs( hin(l))**2
  END DO                         
  DO l=1,512
       hin(l)=hin(l)/sqrt(sum)
  END DO
  END SUBROUTINE norm
END PROGRAM