Appendix A – Source Code for Numerical Experiments
!***********************************************************************
! Program WFSolver
!
! Uses an explicit scheme (based on
delta_t increments) to numerically
! solve the 2D Schroedinger equation for
an arbitrary potential.
!
! The potential is hard-coded (in this
version... that's easy to fix!)
! as an array V(x,y).
!
! The initial wavefunciton is a
Gaussian (this could easily be changed)
! with user specified initial position
and momentum.
!
! The size of the 2D space, the number
of grid points and the timestep are
! all user configurable. The frequency
of output file generation is also
! able to be specified by the user.
!
! Output files are named 'slitXXXX.da'
for backwards compatibility (where
! XXXX is the frame number, padded with
zeros). The total probability
! at each frame is written out to the
file 'slitproby.dat'. The potential
! is saved to 'slitp.dat'. Each of the
'.da' files (and 'slitp.dat') is in
! an appropriate format for use by
GNUplot with the command:
! splot
[0:45] [0:30] [0:0.5] 'slit0001.da' w l , 'slitp.dat' w l
!
! (Note that the specifying the size of
the X, Y and Z dimensions speeds the
! plotting process, but these values
may vary depending on the size of the
! user-defined space.)
!
!
! May be built under un*x by:
!
f90 wfsolver.f90
!
! or using the makefile, do a compile
(if necessary) and load input from
! the file wf.in by:
! make
run
!
!
! Copyright (c) 1998 Stuart Prescott
(daedelus@vislab.usyd.edu.au)
! Source code freely distributable
under the GNU GPL
!
! The results were generated using the
Sydney Regional Visualisation
! Laboratory which is supported by the
Australian Research Council.
!***********************************************************************!
Program wfsolver
Implicit
None
!Set
up arrays for the real and imaginary parts of the WF:
Real,
Dimension(:,:,:), Allocatable:: psr, psi
!Also
need a potential (time indep, so don't need extra dimension)
!and
for the proby (mod(psi))
Real,
Dimension(:,:), Allocatable:: v, ps2
!The
total probability
Real
proby
!Dimensions
of the 2D space
Real
xsize, ysize
!Some
parameters: max=number of mesh elements in each coord
!framestep=steps
between output writes
Integer
maxx, maxy, framestep, time
!Generate
GNUPlot animation file
character*1
dognuplot
Integer
dognuflag
!Also
set up some useful dummy vars
Integer
i, j, n, framenum, writestepi, writestepj
!Used
for generating the initial wavepacket
Complex
exc, zi
!Parameters
for setting up the potential
Real
va, vb, vbig
!Other
parameters in the wf solving routine
Real
dx, dy, dt, dxdy, x0, y0, k0x, k0y, x, y, a1, a2, alpha
!The
filename to write things to
Character*11
filename
!Get
the number of time steps and the number of output frames
Write(*,*)'Enter
the number of time steps to iterate through'
Write(*,*)'(must
be a positive integer)'
Read(*,*)time
Write(*,*)'Enter
the number of time steps to each output frame'
Write(*,*)'(must
be a positive integer)'
Read(*,*)framestep
!SIZE
OF GRID
Write(*,*)'Enter
the number of grid points'
Read(*,*)maxx
maxy
= maxx
Allocate
(psr(maxx+1,maxy+1,2), psi(maxx+1,maxy+1,2))
Allocate
(v(maxx+1,maxy+1), ps2(maxx+1,maxy+1))
Write(*,*)'Enter
the size of the space'
Read(*,*)xsize
ysize
= xsize
dx
= xsize/real(maxx)
dy
= ysize/real(maxy)
!Parameters
for how the info will be written out
writestepi=CEILING(real(maxx)/45.0)
writestepj=CEILING(real(maxy)/30.0)
write(*,*)
'writesteps',writestepi, writestepj
Write(*,*)'Enter
the timestep'
Read(*,*)dt
dxdy = dx*dy
alpha
= 0.5*dt/(dx*dx)
!initial
momentum, position
Write(*,*)'Enter
the initial x momentum'
Read(*,*)k0x
Write(*,*)'Enter
the initial y momentum'
Read(*,*)k0y
Write(*,*)'Enter
the initial x position'
Read(*,*)x0
Write(*,*)'Enter
the initial y position'
Read(*,*)y0
zi = cmplx(0.0D0,1.D0)
!initial
wave function
y
= -ysize/2.0
Do
j=1, maxy+1
x = -xsize/2.0
Do i=1, maxx+1
exc
= exp(zi*(k0x*x+k0y*y))
a1 =
exp(-0.5*(((x-x0))**2.0+((y-y0))**2.0))
!real part of the initial wave function
psr(i,j,1) = REAL(a1*exc)
!imaginay part of the initial wave
function
psi(i,j,1) = AIMAG(a1*exc)
x
= x + dx
EndDo
y = y + dx
EndDo
!set
up the potential slit width: 50-40=10 units
!Do j=1,maxy+1
! Do i=1,maxx+1
!
If((j.eq.35).and.((i.le.40).or.(i.ge.51)))Then
! v(i,j) = 50.0
! Else
! v(i,j) = 0.0
! EndIf
! EndDo
!EndDo
!Set
up a ramp V
!Do
j=1,maxy+1
! v(:,j)=50.0*(1.0-real(j)/max)
!EndDo
!Set
up a completely flat potential
!v(:,:)=0.3
!Set
up a paraboloid potential
!y=-ysize/2
!Do
j=1,maxy+1
! x=-xsize/2
! Do i=1,maxx+1
! v(i,j)=1.0*(x**2+y**2)
! x=x+dx
! EndDo
! y=y+dy
!EndDo
!Set
up a parabolic trough potential
!y=-ysize/2
!Do
j=1,maxy+1
! v(:,j)=0.8*(y**2)
! y=y+dy
!EndDo
!OK,
how about a stadium problem
Write(*,*)'
'
Write(*,*)'Stadium
Setup'
Write(*,*)'Enter
the cap radius'
Read(*,*)va
Write(*,*)'Enter
the side length'
Read(*,*)vb
Write(*,*)'Enter
the outer potential'
Read(*,*)vbig
y
= -ysize/2.0
Do
j=1, maxy+1
If (abs(y) .gt. vb) Then
v(:,j)=vbig
Else
x=-xsize/2.0
Do i=1,maxx+1
If ( (abs(x) .gt. va) .and. &
((x-va)**2.0 + (y)**2.0 .gt. (vb)**2.0)
.and. &
((x+va)**2.0 + (y)**2.0 .gt.
(vb)**2.0) ) Then
v(i,j)=vbig
!*2/3
Else
v(i,j)=0.0
!vbig/3 !=0.0
EndIf
x=x+dx
EndDo
EndIf
y=y+dy
EndDo
!Save
the potential
open(9,FILE='slitp.dat',STATUS='UNKNOWN')
Do
j=2, maxy, writestepj
Do i=2, maxx, writestepi
Write(9,11)v(i,j)*0.005
EndDo
Write(9,11)
EndDo
Close(9)
write
(*,*) 'Potential saved in slitp.dat'
open(10,
FILE='slitproby.dat', STATUS='UNKNOWN')
!
propagate solution through time
Do n=1,time
!calculate the real part of the wavefunction
Do j=2, maxy
Do i=2, maxx
psr(i,j,2)=psr(i,j,1)+2*psi(i,j,1)*(4*alpha+dt*v(i,j)) &
-2*alpha*(psi(i+1,j,1)+psi(i-1,j,1)+psi(i,j+1,1)+psi(i,j-1,1))
EndDo
!force derivative at x-edges to be
zero
psr(1, j, 2) = psr(2, j,
2)
psr(maxx+1, j, 2) = psr(maxx, j, 2)
EndDo
!now do the imaginary part of the
wavefunction
Do j=2, maxy
Do i=2, maxx
psi(i,j,2)=psi(i,j,1)-2*psr(i,j,2)*(4*alpha+dt*v(i,j)) &
+2*alpha*(psr(i+1,j,2)+psr(i-1,j,2)+psr(i,j+1,2)+psr(i,j-1,2))
EndDo
!force derivative at x-edges to be
zero
psi(1, j, 2) = psi(2, j,
2)
psi(maxx+1, j, 2) = psi(maxx, j, 2)
EndDo
psi(:,:,1) = psi(:,:,2)
psr(:,:,1) = psr(:,:,2)
if (mod(n,framestep) .eq. 0) then
framenum=floor(real(n)/real(framestep))
write(*,*)'Writing for frame ',framenum
filename='slit'//&
&
CHAR(mod(framenum,10000)/1000+ICHAR('0'))// &
&
CHAR(mod(framenum,1000 )/100 +ICHAR('0'))// &
&
CHAR(mod(framenum,100 )/10 +ICHAR('0'))// &
&
CHAR(mod(framenum,10 ) +ICHAR('0'))// &
&
'.dat'
Open(9,FILE=filename, STATUS='UNKNOWN')
Do j=2,maxy,writestepj
Do i=2,maxx,writestepi
ps2(i,j)
= psr(i,j,1)*psr(i,j,1)+ &
psi(i,j,1)*psi(i,j,1)
Write(9,11)ps2(i,j)
EndDo
!Extra line for GNUplot format
Write(9,11)
EndDo
Close(9)
write(*,*)'Frames saved at ',n
!Also need to keep track of
probability....
proby=0
Do j=2, maxy
Do i=2, maxx
proby
= proby + ps2(i,j)*dxdy
EndDo
EndDo
write(10,11)proby
EndIf
11 Format (E12.6)
EndDo
close(10) !probability file
End