Program HomeworkThreeThree use fund_const use MSIMSL implicit none real(8) dLength, dK, k0, variance, kStart, phi(sites), k, kappa, E, x, t, pos real(8) barrierStart, trajStart, dTraj, lastTraj(trajCount), wignerR complex(8) wavefunction(sites), F, B, C, D, bar1, bar2, Wigner integer i, j, r, thisIndex, barrierStartIndex, barrierEndIndex, wigner1, wigner2 ! Initialize some convenience variables dLength = deviceLength / sites dK = kLength / sites 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 t = totalTime 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 (10, file='psi.dat') do i=1,sites wavefunction(i) = wavefunction(i) * dK / dsqrt(c2*pi) write (10,*) dreal(wavefunction(i)), dimag(wavefunction(i)) end do close(10) ! Find the expectation of position wignerR = c0 do i=1,sites x = -deviceLength/c2 + (i-1)*dLength wignerR = wignerR + dconjg(wavefunction(i)) * wavefunction(i) * x end do wignerR = wignerR * dLength ! Find the Wigner function open (10, file='wigner.dat') do i=1,sites k = kStart + (i-1)*dk Wigner = c0 do j=1,sites pos = -deviceLength/c2 + (j-1)*dLength wigner1 = nint((wignerR - pos/c2 + deviceLength/c2)/dLength) wigner2 = nint((wignerR + pos/c2 + deviceLength/c2)/dLength) if (((wigner1.ge.1).and.(wigner1.le.sites)).and.((wigner2.ge.1).and.(wigner2.le.sites))) then Wigner = Wigner + dconjg(wavefunction(wigner2)) * wavefunction(wigner1) * cdexp(cimag * k * pos) end if end do Wigner = Wigner * dLength / pi / hbar write (10,*) dreal(Wigner), dimag(Wigner) end do close(10) 97 format(A30,ES13.5) end