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, trajStart, dTraj, lastTraj(trajCount)
	complex(8) wavefunction(sites), gradPsi(sites), F, B, C, D, bar1, bar2
	integer i, j, p, r, thisIndex, barrierStartIndex, barrierEndIndex




	! 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
	trajStart = -deviceLength/c2 + initialOffset - 3.0*stdDeviation
	dtraj = 6.0*stdDeviation / trajCount
	barrierStartIndex = nint((deviceLength/c2 - barrierWidth/c2)/dLength)
	barrierEndIndex = nint((deviceLength/c2 + barrierWidth/c2)/dLength)


	! 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)
	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')
	open (11, file='t.dat')
	do p=1,timeSteps
		t = (p-1)*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
			F = phi(i) * cdexp(-cimag*k*barrierWidth) * c1/(dcosh(kappa*barrierWidth) - cimag/c2*(k**2-kappa**2)		&
					/k/kappa*dsinh(kappa*barrierWidth))
			B = -cimag/c2 * F * (k**2 + kappa**2)/k/kappa * dsinh(kappa*barrierWidth)
			do j=1,barrierStartIndex-1
				x = -deviceLength/c2 + (j-1)*dLength
				bar1 = phi(i) * cdexp(cimag * (k * (x + deviceLength/c2 - initialOffset) - hbar * k**2 *t/c2/m))		&		
								+ B * cdexp(cimag * (-k * (x - deviceLength/c2 + initialOffset) - hbar * k**2 *t/c2/m))
				wavefunction(j) = wavefunction(j) + bar1
			end do
			do j=sites,barrierEndIndex+1,-1
				x = -deviceLength/c2 + (j-1)*dLength
				bar2 = F * cdexp(cimag * (k * (x + deviceLength/c2 - initialOffset) - hbar * k**2 *t/c2/m))
				wavefunction(j) = wavefunction(j) + bar2
			end do
			do j=barrierStartIndex,barrierEndIndex
				x = -deviceLength/c2 + (j-1)*dLength
				wavefunction(j) = wavefunction(j) + bar1 * dexp(kappa * x)										&										
									+ bar2 * dexp(-kappa*x)
			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)



		! Find the derivative of the wavefunction
		gradPsi(1) = wavefunction(2) - wavefunction(1)
		do r=2,sites
			gradPsi(r) = wavefunction(r) - wavefunction(r-1)
		end do

		! Get the trajectories
		do r=1,trajCount
			if (p.eq.1) then
				lastTraj(r) = trajStart + dTraj * (r-1)
			else
				thisIndex = nint((lastTraj(r) + deviceLength/c2)/dLength)
				if ((thisIndex.le.sites).and.(thisIndex.gt.0)) then
					lastTraj(r) = lastTraj(r) + dreal(hbar/c2/m/cimag * (dconjg(wavefunction(thisIndex)) * gradPsi(thisIndex)	&
								- wavefunction(thisIndex) * dconjg(gradPsi(thisIndex)))											&
								/ dconjg(wavefunction(thisIndex)) / wavefunction(thisIndex)) * dTime / dLength
				end if
			end if
			write (10, '(ES13.5)', advance='NO') lastTraj(r)
		end do
		write(10,*)
		write (11,*) t



	end do
	close(10)
	close(11)

97	format(A30,ES13.5)
end