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