Program HomeworkThree
use fund_const
use MSIMSL

	implicit none
	real(8) dLength, dK, dTime, k0, variance, kStart, phi(sites), k, kappa, E, x, t
	real(8) barrierStart
	complex(8) wavefunction(sites), trans, refl
	integer i, j, p, r




	! Initialize some convenience variables
	dLength = deviceLength / sites
	dK = kLength / sites
	dTime = totalTime / timeSteps
	k0 = dsqrt(c2*m*initialEnergy)/hbar
	variance = stdDeviation**2
	kStart = k0 - kLength/c2
	barrierStart = deviceLength/c2




	! Create our k-space Gaussian at t=0
	open (10, file='k.dat')
	open (11, file='x.dat')
	open (12, file='phi.dat')
	do i=1,sites
		k = kStart + (i-1)*dk
		x = -deviceLength/c2 + (i-1)*dLength
		phi(i) = (variance/c2/pi)**0.25 * dexp(-variance * (k - k0)**2/4)
		write (10,*) k
		write (11,*) x
		write (12,*) phi(i), phi(i), phi(i), phi(i), phi(i), phi(i)
	end do
	close(10)
	close(11)
	close(12)




	! Create our real-space Gaussian from our k-space Gaussian at each time step
	open (10, file='traj.dat')
	do p=0,timeSteps
		t = p*dTime
		wavefunction = c0
		do i=1,sites
			k = kStart + (i-1)*dk
			E = hbar**2 * k**2 / c2 / m
			kappa = sqrt(c2*m*(barrierHeight - E))/hbar
			trans = phi(i) * cdexp(-cimag*k*barrierWidth) * c1/(dcosh(kappa*barrierWidth) + (k**2+kappa**2)/c2		&
					/k/kappa*dsinh(kappa*barrierWidth))
			refl = phi(i) * cdsqrt(1 - (trans/phi(i))**2)
			do j=1,sites
				x = -deviceLength/c2 + (j-1)*dLength
				if (x.lt.c0 - barrierWidth/c2) then
					wavefunction(j) = wavefunction(j)																&
								+ phi(i) * cdexp(cimag * (k * (x + deviceLength/c2 - initialOffset) - hbar * k**2 *t/c2/m))		&		
								+ refl * cdexp(cimag * (-k * (x - deviceLength/c2 + initialOffset) - hbar * k**2 *t/c2/m))
				else if (x.gt.c0 + barrierWidth/c2) then
					wavefunction(j) = wavefunction(j)																&
								+ trans * cdexp(cimag * (k * (x + deviceLength/c2 - initialOffset) - hbar * k**2 *t/c2/m))
				end if
			end do
		end do
		open (100+p)
		do i=1,sites
			wavefunction(i) = wavefunction(i) * dK / dsqrt(c2*pi)
			write (100+p,*) dreal(wavefunction(i)), dimag(wavefunction(i))
		end do
		close(100+p)
	end do
	close(10)




97	format(A30,ES13.5)
end