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