c     Program lattice_osc.f    by Gilmore Jan 14, 2009

c     The tridiagonal symmetric matrix describing normal
c     mode oscillations (d(i)=+2, e(i)=-1) is created
c     and diagonalized.  the energies and lowest three
c     eigenvectors are written to an output file.

      implicit   none
      integer    i,j,k,n,np
      parameter(n=99,np=n)
      real*8     xleft,xright,x,delta,ham(np,np)
      real*8     d(np),e(np),z(np,np),a(np,np),norm
      real*8     pythag
      external   pythag


c     begin



      do i=1,n
      do j=1,n
        z(i,j) = 0.0
      end do
      end do

      do i=1,n
        d(i) = 2.0
        e(i) =-1.0
        z(i,i) = 1.0
      end do

      CALL tqli(d,e,n,np,z)
      CALL eigsrt(d,z,n,np)

      write(*,'(2x,11f12.6)')(d(i),i=1,n)

      open(11,file='harm_osc_e.dat',status='unknown')
      do i=1,n
      write(11,'(i6,f12.6)')i,d(i)
      end do
      close(11)

      open(13,file='harm_osc_v.dat',status='unknown')
      do i=1,n
      write(13,'(i6,3f12.6)')i,(z(i,j),j=1,3)
      end do
      close(11)



900   stop
      end



      SUBROUTINE tqli(d,e,n,np,z)
      INTEGER n,np
      REAL*8    d(np),e(np),z(np,np)
CU    USES pythag
      INTEGER    i,iter,k,l,m
      REAL*8   b,c,dd,f,g,p,r,s,pythag,one
      do 11 i=2,n
        e(i-1)=e(i)
11    continue
      e(n)=0.
      do 15 l=1,n
        iter=0
1       do 12 m=l,n-1
          dd=abs(d(m))+abs(d(m+1))
          if (abs(e(m))+dd.eq.dd) goto 2
12      continue
        m=n
2       if(m.ne.l)then
          if(iter.eq.30)pause 'too many iterations in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.*e(l))
          one = 1.0
          r=pythag(g,one)
          g=d(m)-d(l)+e(l)/(g+sign(r,g))
          s=1.
          c=1.
          p=0.
          do 14 i=m-1,l,-1
            f=s*e(i)
            b=c*e(i)
            r=pythag(f,g)
            e(i+1)=r
            if(r.eq.0.)then
              d(i+1)=d(i+1)-p
              e(m)=0.
              goto 1
            endif
            s=f/r
            c=g/r
            g=d(i+1)-p
            r=(d(i)-g)*s+2.*c*b
            p=s*r
            d(i+1)=g+p
            g=c*r-b
C     Omit lines from here ...
            do 13 k=1,n
              f=z(k,i+1)
              z(k,i+1)=s*z(k,i)+c*f
              z(k,i)=c*z(k,i)-s*f
13          continue
C     ... to here when finding only eigenvalues.
14        continue
          d(l)=d(l)-p
          e(l)=g
          e(m)=0.
          goto 1
        endif
15    continue
      return
      END


      FUNCTION pythag(a,b)
      REAL*8    a,b,pythag
      REAL*8    absa,absb
      absa=abs(a)
      absb=abs(b)
      if(absa.gt.absb)then
        pythag=absa*sqrt(1.+(absb/absa)**2)
      else
        if(absb.eq.0.)then
          pythag=0.
        else
          pythag=absb*sqrt(1.+(absa/absb)**2)
        endif
      endif
      return
      END



      SUBROUTINE eigsrt(d,v,n,np)
      INTEGER n,np
      REAL*8     d(np),v(np,np)
      INTEGER    i,j,k
      REAL*8     p
      do 13 i=1,n-1
        k=i
        p=d(i)
        do 11 j=i+1,n
          if(d(j).le.p)then   !!!
            k=j
            p=d(j)
          endif
11      continue
        if(k.ne.i)then
          d(k)=d(i)
          d(i)=p
          do 12 j=1,n
            p=v(j,i)
            v(j,i)=v(j,k)
            v(j,k)=p
12        continue
        endif
13    continue
      return
      END


      SUBROUTINE tred2(a,n,np,d,e)
      INTEGER n,np
      REAL*8     a(np,np),d(np),e(np)
      INTEGER    i,j,k,l
      REAL*8     f,g,h,hh,scale
      do 18 i=n,2,-1
        l=i-1
        h=0.
        scale=0.
        if(l.gt.1)then
          do 11 k=1,l
            scale=scale+abs(a(i,k))
11        continue
          if(scale.eq.0.)then
            e(i)=a(i,l)
          else
            do 12 k=1,l
              a(i,k)=a(i,k)/scale
              h=h+a(i,k)**2
12          continue
            f=a(i,l)
            g=-sign(sqrt(h),f)
            e(i)=scale*g
            h=h-f*g
            a(i,l)=f-g
            f=0.
            do 15 j=1,l
C     Omit following line if finding only eigenvalues
              a(j,i)=a(i,j)/h
              g=0.
              do 13 k=1,j
                g=g+a(j,k)*a(i,k)
13            continue
              do 14 k=j+1,l
                g=g+a(k,j)*a(i,k)
14            continue
              e(j)=g/h
              f=f+e(j)*a(i,j)
15          continue
            hh=f/(h+h)
            do 17 j=1,l
              f=a(i,j)
              g=e(j)-hh*f
              e(j)=g
              do 16 k=1,j
                a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
16            continue
17          continue
          endif
        else
          e(i)=a(i,l)
        endif
        d(i)=h
18    continue
C     Omit following line if finding only eigenvalues.
      d(1)=0.
      e(1)=0.
      do 24 i=1,n
C     Delete lines from here ...
        l=i-1
        if(d(i).ne.0.)then
          do 22 j=1,l
            g=0.
            do 19 k=1,l
              g=g+a(i,k)*a(k,j)
19          continue
            do 21 k=1,l
              a(k,j)=a(k,j)-g*a(k,i)
21          continue
22        continue
        endif
C     ... to here when finding only eigenvalues.
        d(i)=a(i,i)
C     Also delete lines from here ...
        a(i,i)=1.
        do 23 j=1,l
          a(i,j)=0.
          a(j,i)=0.
23      continue
C     ... to here when finding only eigenvalues.
24    continue
      return
      END

