Program GaussianDispersion
use fund_const
use MSIMSL

	implicit none
	integer i, j, xMode, yMode, tIndex
	real(8) dLength, dTime, sum, startX, startY, sigma, gausNorm, k0, kx, ky, x, y, t, E
	complex(8) weight(modeCount, modeCount), wavefunction(sites, sites)




	! Initialize some convenience variables
	dLength = deviceLength / sites
	dTime = totalTime / timeSteps
	startX = deviceLength / 2.0
	startY = deviceLength / 2.0
	sigma = nominalWidth / 6.0
	k0 = mGaAs * initialVelocity / hbar
	kx = 1.0/sqrt(velocitySlope**2 + 1.0) * k0
	ky = velocitySlope * kx
			



	! Write out our coordinate system for convenience
	open (10, file='x.dat')
	open (11, file='y.dat')
	do i=1,sites
		x = i*dLength
		y = i*dLength
		write (10,*) x
		write (11,*) y
	end	do
	close(10)
	close(11)





	! Normalize our Gaussian
!	k0 = 0; kx = 0; ky = 0
	sum = 0
	do i=1,sites
		do j=1,sites
			sum = sum + exp(-2.0 * ( (i*dLength - startX)**2 + (j*dLength - startY)**2 ) / sigma**2)
		end do
	end do
	sum = sum * dLength**2
	gausNorm = 1.0/sqrt(sum)




	! Compute the probabilities at t=0
	weight = c0
	open (10, file='weight.dat')
	do xMode = 1,modeCount
		do yMode = 1,modeCount
			do i=1,sites
				do j=1,sites
					x = i*dLength
					y = j*dLength
					weight(xMode, yMode) =	weight(xMode, yMode) + 												&
											exp(- ((x - startX)**2 + (y - startY)**2) / sigma**2) *				&
											sin(xMode * x * pi/deviceLength) * sin(yMode * y * pi/deviceLength)
				end do
			end do
			weight(xMode, yMode) = weight(xMode, yMode) * 2.0/deviceLength * gausNorm * dLength**2
			write (10,*) real(weight(xMode, yMode)), imag(weight(xMode, yMode))
		end do
	end do
	close(10)





	! Compute the superposition wavefunction at various times
	do tIndex = 0,timeSteps
		t = tIndex * dTime
		wavefunction = c0
		write(*,97) 'time:',t
		do xMode = 1,modeCount
			do yMode = 1,modeCount
				do i=1,sites
					do j=1,sites
						x = i*dLength
						y = j*dLength
						E = (xMode + yMode)**2 * hbar/deviceLength*hbar/deviceLength * pi**2 / 2.0 / mGaAs
						wavefunction(i,j) = wavefunction(i,j) + weight(xMode, yMode) *								&
											2.0/deviceLength * exp(cimag * (kx*x + ky*y - E * t /hbar)) *			&
											dsin(xMode * x * pi/deviceLength) * dsin(yMode * y * pi/deviceLength)
					end do
				end do
			end do
		end do
		open (100+tIndex)
		do i=1,sites
			do j=1,sites
				write (100+tIndex,*) real(wavefunction(i,j)), imag(wavefunction(i,j))
			end do
		end	do
		close(100+tIndex)
	end do


97	format(A30,ES13.5)
end