PROGRAM Quantum
! for parts (1) and (2) set Xmax = 4000 and Xmin to 0
! for extra part, set xmax = 5000 and Xmin to Xmin to -1000
IMPLICIT NONE
INTEGER, PARAMETER :: M = 100, Xmax = 4000, NumTimes = 6, Xmin = 0
COMPLEX*8, PARAMETER :: Imag = 1.0d0 * (0.0d0,1.0d0)
REAL*8, PARAMETER :: HBar = 1.0d0, Mass = 1.0d0, Two = 2.0d0, H = 0.01d0 
REAL*8 :: K0, Lambda0, Pi, A, B, DeltaK = 0.0d0, Time = 0.0d0
REAL*8 :: Vp = 0.0d0, Vg = 0.0d0, Omega, TwoPi, X0, sink = 0.0d0
INTEGER :: I, J, T
REAL*8, DIMENSION(-M:M) :: Kj, Aj, Bj
REAL*8, DIMENSION(Xmin:XMax) :: X
REAL*8, DIMENSION(0:NumTimes,Xmin:XMax) :: AbsPsi = 0.0d0
COMPLEX*8, DIMENSION(-M:M) :: Sum1, Sum2
COMPLEX*8, DIMENSION(Xmin:XMax) :: Psi = (0.0d0,0.0d0)

OPEN (UNIT = 10, FILE = "Quantum.txt", STATUS = "UNKNOWN")
Pi = DACOS(-1.0d0)
WRITE(*,100, ADVANCE = "NO"), "Please input the central wavelength: "
READ *, Lambda0
WRITE(*,100, ADVANCE = "NO"), "Please input the full width at half maximum: "
READ *, B
TwoPi = 2.0d0 * DACOS(-1.0d0) ! calculate 2pi
K0 = (2.0d0 * Pi) / Lambda0 ! calculate k0
A = B / (Two * DSQRT(DLOG(Two))) ! calculate A

X0 = 2.0d0 ! comment out this line for parts (1) and (2)
Vp = (Hbar * K0) / (2.0d0 * Mass) ! phase velocity
Vg = (HBar * K0) / Mass ! group velocity
Omega = (Hbar * K0 ** 2) / (2.0d0 * Mass) ! calculate omega
PRINT *, "Phase velocity is: ", Vp
PRINT *, "Group velocity is: ", Vg
PRINT *, "The period is: ", TwoPi / Omega
DeltaK = 4.0d0 / (A * M) ! calculate delta k
DO J = -M, M
     Kj(J) = K0 + J * DeltaK 
     ! to do parts (1) and (2), comment out lines 40 through 45, to do extra part
     ! comment out line 39 and uncomment lines 40 through 45
      Aj(J) = A / (2.0d0 * DEXP((A ** 2 * (Kj(J) - K0) ** 2) / 4.0d0) * DSQRT(Pi))
     ! Sink = 1.0d0 / 2.0d0 * (Kj(J) - K0) * X0
     ! IF (J == 0) THEN
     !      Aj(J) = X0
     !      CYCLE
     ! END IF
     ! Aj(J) = X0 * Sinc(Sink)
END DO
DO I = Xmin, XMax
     X(I) = I * H ! calculate position
END DO
OPEN (UNIT = 25, file = "vel.txt")

DO T = 0, NumTimes
     Time = DBLE(T) ! convert t to double precision
     DO J = -M, M
          Bj(J) = Kj(J) ** 2 * Time / 2.0d0 ! calcualte bj
     END DO
     DO I = Xmin, XMax
          Psi(I) = (0.0d0,0.0d0) ! reset Psi
     END DO
     DO I = Xmin, XMax
          DO J = -M, M
               Psi(I) = Psi(I) + Aj(J) * EXP(IMag * (Kj(J) * X(I)&
               &- Bj(J))) * DeltaK ! calculate Psi
          END DO
     END DO
     DO I = Xmin, XMax
          AbsPsi(T,I) = Psi(I) * CONJG(Psi(I)) ! find |Psi|^2
     END DO
END DO
DO I = Xmin, XMax
     WRITE(10,200, ADVANCE = 'NO') X(I) ! output position
     DO T = 0, NumTimes
          ! output |Psi(x,t)|^2 for each t to file, all on same line
          WRITE(10,210, ADVANCE = 'NO') AbsPsi(T,I) 
     END DO
     WRITE(10,*) ! new ine
END DO

100 FORMAT(1X, A:)
200 FORMAT(1X, F10.6)
210 FORMAT(1X, ES12.5)

! --------------------------------------------------------------------
 CONTAINS
! --------------------------------------------------------------------

FUNCTION Sinc(X)
REAL*8 :: X, Sinc
! Sinc function
! sinc(x) = sin(x) / x
Sinc = DSIN(X) / X
END FUNCTION

! --------------------------------------------------------------------

END PROGRAM