c bands.f
c computes band structure of 2-dimension kronig-penny model.
c r. b. laughlin - 24 jan 07
      parameter(n=20)
      dimension f(n),e(n)
c establish desired numerical accuracy.
      eps=0.0001
      s=sqrt(2.0)
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 open output file.
      open(7,file=\"out1.dat\",status=\"new\")
      open(8,file=\"out2.dat\",status=\"new\")
c loop q from gamma to x.
      do 2 j=1,40
      qx=0.0125*(j-0.5)
      qy=0.0
      call bands(qx,qy,e,f)
      x=qx
      write(7,20) x,(e(k), k=1,n)
    2 continue
      write(7,21)
      x=0.5
      write(8,20) x
c loop q from x to m.
      do 3 j=1,40
      qy=0.0125*(j-0.5)
      call bands(qx,qy,e,f)
      x=0.5+qy
      write(7,20) x,(e(k), k=1,n)
    3 continue
      write(7,21)
      x=1.0
      write(8,20) x
c loop q from m to gamma.
      do 4 j=1,40
      qx=0.5-0.0125*(j-0.5)
      qy=qx
      call bands(qx,qy,e,f)
      x=1.0+s*(0.5-qx)
      write(7,20) x,(e(k), k=1,n)
    4 continue
      x=1.0+0.5*s
      write(8,20) x
   20 format(30f10.4)
   21 format(1h )
      stop
      end

      subroutine bands(qx,qy,e,f)
      parameter(n=20)
      dimension f(n),e(n)
      complex z,g
c establish desired numerical accuracy.
      eps=0.0001
c establish value for potential. warning: must not be zero.
      v=-0.05
      vi=1.0/v
c initialize energy band counter.
      l=1
c initialize old green's function to a negative value.
      grl=-eps
c now hunt through energies to find eigenvalues.
      de=0.002
      do 4 i=1,4000
      e0=(i-1)*de-1.0
c for this energy, initialize new green's function to zero
      z=e0+cmplx(0.0,eps)
      g=0.0
c and then update by summing over reciprocal lattice.
      nn=2*n-1
      do 1 jj=1,nn
      do 1 kk=1,nn
      j=jj-n
      k=kk-n
      x=qx+j
      y=qy+k
      if(j.lt.0) j=-j
      if(k.lt.0) k=-k
      g=g+f(j+1)*f(k+1)/(z-(x**2+y**2)/2.0)
    1 continue
c energy is found if real part of g crossed the axis downward.
      gr=real(g)
      gi=aimag(g)
      if(gr.ge.vi) go to 3
      if(grl.le.vi) go to 2
      e(l)=e0-(gr-vi)/(gr-grl)*de
      l=l+1
c break out if we found enough bands.
      if(l.gt.n) go to 5
    2 continue
    3 continue
c remember most recent value of g.
      grl=gr
    4 continue
    5 continue
      write(*,6) qx,qy,l
    6 format(2f10.4,i5)
      return
      end