c contour.f
c makes contour plot for 2-dimension kronig-penny model.
c r. b. laughlin - 24 jan 07
parameter(n=20)
parameter(np=81)
dimension f(n),p(np,np)
complex psi(np,np)
nn=2*n-1
c establish value for pi.
pi=3.1415927
pi2=2.0*pi
c establish numerical accuracy.
eps=0.0001
c establish value of e.
e0=-0.5658
c establish values of qx and qy.
q0x=0.0
q0y=0.0
c open output file
open(7,file='output.dat',status='new')
c fill reciprocal lattice regularization array.
x=n**2
beta=-2.0*alog(eps)/x
do 1 j=1,n
x=j-1
f(j)=exp(-beta*x**2/2.0)
1 continue
c normalize wavefunction.
b=0.0
do 2 jj=1,nn
do 2 kk=1,nn
j=jj-n
k=kk-n
qx=q0x+j
qy=q0y+k
a=e0-(qx**2+qy**2)/2.0
a=1.0/(a**2+eps**2)
if(j.lt.0) j=-j
if(k.lt.0) k=-k
b=b+f(j+1)*f(k+1)*a
2 continue
b=1.0/sqrt(b)
c calculate space increment.
x=np-1
dx=1.0/x
c zero out charge density accumulator.
do 3 l=1,np
do 3 m=1,np
psi(l,m)=0.0
3 continue
c loop on momenta.
do 5 jj=1,nn
do 5 kk=1,nn
j=jj-n
k=kk-n
qx=q0x+j
qy=q0y+k
a=e0-(qx**2+qy**2)/2.0
a=a/(a**2+eps**2)*b
if(j.lt.0) j=-j
if(k.lt.0) k=-k
c=f(j+1)*f(k+1)*a
c for this momentum, loop over sites.
do 4 l=1,np
x=dx*(l-1)-0.5
do 4 m=1,np
y=dx*(m-1)-0.5
z=(x*qx+y*qy)*pi2
psi(l,m)=psi(l,m)+c*cexp(cmplx(0.0,z))
4 continue
5 continue
c square wavefunction to make probability.
do 6 l=1,np
do 6 m=1,np
p(l,m)=real(psi(l,m))**2+aimag(psi(l,m))**2
6 continue
s=0.5
c pick contour.
do 24 ip=1,20
p0=2.0*ip
c loop over little boxes.
do 23 l=2,np
do 23 m=2,np
pp=(p(l,m)+p(l-1,m)+p(l,m-1)+p(l-1,m-1))/4.0
c first little triangle.
i=1
x=(p0-p(l,m))*(p0-p(l-1,m))
if (x.gt.0.0) go to 11
x=(p0-p(l,m))/(p(l-1,m)-p(l,m))
x=(l-1-x)*dx
y=(m-1)*dx
write(7,30) x,y
i=i+1
11 continue
x=(p0-p(l,m))*(p0-pp)
if(x.gt.0.0) go to 12
x=s*(p0-p(l,m))/(pp-p(l,m))
y=(m-1-x)*dx
x=(l-1-x)*dx
write(7,30) x,y
i=i+1
12 continue
x=(p0-p(l-1,m))*(p0-pp)
if(x.gt.0.0) go to 13
x=s*(p0-p(l-1,m))/(pp-p(l-1,m))
y=(m-1-x)*dx
x=(l-2+x)*dx
write(7,30) x,y
i=i+1
13 continue
if(i.gt.1) write(7,31)
c second little triangle.
i=1
x=(p0-p(l,m))*(p0-p(l,m-1))
if (x.gt.0.0) go to 14
x=(p0-p(l,m))/(p(l,m-1)-p(l,m))
y=(m-1-x)*dx
x=(l-1)*dx
write(7,30) x,y
i=i+1
14 continue
x=(p0-p(l,m))*(p0-pp)
if(x.gt.0.0) go to 15
x=s*(p0-p(l,m))/(pp-p(l,m))
y=(m-1-x)*dx
x=(l-1-x)*dx
write(7,30) x,y
i=i+1
15 continue
x=(p0-p(l,m-1))*(p0-pp)
if(x.gt.0.0) go to 16
x=s*(p0-p(l,m-1))/(pp-p(l,m-1))
y=(m-2+x)*dx
x=(l-1-x)*dx
write(7,30) x,y
i=i+1
16 continue
if(i.gt.1) write(7,31)
c third little triangle.
i=1
x=(p0-p(l,m-1))*(p0-p(l-1,m-1))
if (x.gt.0.0) go to 17
x=(p0-p(l,m-1))/(p(l-1,m-1)-p(l,m-1))
x=(l-1-x)*dx
y=(m-2)*dx
write(7,30) x,y
i=i+1
17 continue
x=(p0-p(l,m-1))*(p0-pp)
if(x.gt.0.0) go to 18
x=s*(p0-p(l,m-1))/(pp-p(l,m-1))
y=(m-2+x)*dx
x=(l-1-x)*dx
write(7,30) x,y
i=i+1
18 continue
x=(p0-p(l-1,m-1))*(p0-pp)
if(x.gt.0.0) go to 19
x=s*(p0-p(l-1,m-1))/(pp-p(l-1,m-1))
y=(m-2+x)*dx
x=(l-2+x)*dx
write(7,30) x,y
i=i+1
19 continue
if(i.gt.1) write(7,31)
c fourth little triangle.
i=1
x=(p0-p(l-1,m))*(p0-p(l-1,m-1))
if (x.gt.0.0) go to 20
x=(p0-p(l-1,m))/(p(l-1,m-1)-p(l-1,m))
y=(m-1-x)*dx
x=(l-2)*dx
write(7,30) x,y
i=i+1
20 continue
x=(p0-p(l-1,m))*(p0-pp)
if(x.gt.0.0) go to 21
x=s*(p0-p(l-1,m))/(pp-p(l-1,m))
y=(m-1-x)*dx
x=(l-2+x)*dx
write(7,30) x,y
i=i+1
21 continue
x=(p0-p(l-1,m-1))*(p0-pp)
if(x.gt.0.0) go to 22
x=s*(p0-p(l-1,m-1))/(pp-p(l-1,m-1))
y=(m-2+x)*dx
x=(l-2+x)*dx
write(7,30) x,y
i=i+1
22 continue
if(i.gt.1) write(7,31)
23 continue
24 continue
30 format(2f10.4)
31 format(1h )
stop
end