      program plateeig

c**********************************************************************
c**** Program to find the eigenvalues and eigenvectors for the     ****
c**** rectangular thin elastic plate                               ****
c**** Compiled under gfortran and linked to the NAG library        ****
c**** Input is a via an input file, default is called default.inp  ****
c #nTrunc
c 24
c #nu
c 0.225
c #aob
c 2
c #pTrunc
c 4
c**** Output is a file default.res with the eigenvalues in them    ****
c**** (format: mode #, SS, SA, AS, AA) and files fort.1x where     ****
c**** x=2,3,etc for modes 1,2,etc with data to plot surfaces       ****
c**********************************************************************

      implicit none
      integer          maxNTrunc
      parameter        (maxNTrunc=48)
      integer          maxPTrunc
      parameter        (maxPTrunc=16)

      integer          maxMatSize,LWORK
      parameter        (maxMatSize=(maxNTrunc+1)**2,
     +                  LWORK=5*maxMatSize)
 
      integer          fUnitRes
      integer          nTrunc,pTrunc,matSize,ifail
      integer          n,m,p,q,mnindx,pqindx
      integer          xLoop,yLoop,nXPts,nYPts

      parameter        (nXPts=100,nYPts=100)

      double precision pi
      double precision aob,nu

      double precision kn(0:2*maxNTrunc+1)

      double precision CnE(0:maxNTrunc),CnO(0:maxNTrunc)
      double precision LnE(0:maxNTrunc),LnO(0:maxNTrunc)
      double precision JmpE(0:maxNTrunc,0:maxNTrunc)
      double precision JmpO(0:maxNTrunc,0:maxNTrunc)

      double precision KmnpqSS(maxMatSize,maxMatSize)
      double precision KmnpqSA(maxMatSize,maxMatSize)
      double precision KmnpqAS(maxMatSize,maxMatSize)
      double precision KmnpqAA(maxMatSize,maxMatSize)

      double precision W(maxMatSize)
      double precision lSS(maxPTrunc)
      double precision lSA(maxPTrunc)
      double precision lAS(maxPTrunc)
      double precision lAA(maxPTrunc)
      double precision VSS(maxPTrunc,0:maxNTrunc,0:maxNTrunc)
      double precision VSA(maxPTrunc,0:maxNTrunc,0:maxNTrunc)
      double precision VAS(maxPTrunc,0:maxNTrunc,0:maxNTrunc)
      double precision VAA(maxPTrunc,0:maxNTrunc,0:maxNTrunc)
      double precision WORK(LWORK)

      double precision xoa,yob
      double precision wnEx(0:maxNTrunc),wnEy(0:maxNTrunc)
      double precision wnOx(0:maxNTrunc),wnOy(0:maxNTrunc)
      double precision wSS,wSA,wAS,WAA

      common /const/   pi

c**********************************************************************
c**** Define some constants like pi, i etc.                        ****
c**********************************************************************

      pi = 4.0d0*datan(1.0d0)
      funitRes = 11

c**********************************************************************
c**** Title and data input and the convert theta to radians        ****
c**********************************************************************

      call DATA_INPUT(funitRes,nTrunc,nu,aob,pTrunc)

c**********************************************************************
c**** Calculate roots of plate equation                            ****
c**********************************************************************

      call PROOTS(nTrunc,kn)

c**********************************************************************
c**** Define Cn and Ln for the Even and Odd cases                  ****
c**********************************************************************

      do 10 n=0,nTrunc

        if (n.eq.0) then

          CnE(n) = dsqrt(1.0d0/2.0d0)
          CnO(n) = dsqrt(1.0d0/6.0d0)

          LnE(n) = 0.0d0
          LnO(n) = 0.5d0

        else

          CnE(n) = dsqrt(1.0d0/8.0d0)
          CnO(n) = dsqrt(1.0d0/8.0d0)

          LnE(n) = 0.5d0*kn(2*n)*dtanh(kn(2*n))
          LnO(n) = 0.5d0*kn(2*n+1)/dtanh(kn(2*n+1))

        end if

  10  continue
 
c**********************************************************************
c**** Define Jmp for Even and Odd orders                           ****
c**********************************************************************

      do 30 m=0,nTrunc
        do 20 p=0,nTrunc

          if (m.eq.0) then

            JmpE(m,p) = 0.0d0
            JmpO(m,p) = 0.0d0

          else if (p.eq.0) then

            JmpE(m,p) = LnE(m)-LnE(0)
            JmpO(m,p) = LnO(m)-LnO(0)

          else if (m.eq.p) then

            JmpE(m,m) = 0.25d0*LnE(m)-0.5d0*(LnE(m)**2)

            JmpO(m,m) = 0.25d0*LnO(m)-0.5d0*(LnO(m)**2)

          else

            JmpE(m,p) = (LnE(m)-LnE(p))/(1.0d0-(kn(2*p)/kn(2*m))**4)
            JmpO(m,p) = (LnO(m)-LnO(p))/(1.0d0-(kn(2*p+1)/kn(2*m+1))**4)

          end if

  20    continue
  30  continue
     
c**********************************************************************
c**** Construct the four symmetry matrices                         ****
c**********************************************************************

      do 70 m=0,nTrunc
        do 60 n=0,nTrunc

          mnindx = m*(nTrunc+1)+n+1

          do 50 p=0,nTrunc
            do 40 q=0,nTrunc

              pqindx = p*(nTrunc+1)+q+1

              KmnpqSS(mnindx,pqindx) = (aob**2)*(
     +               nu*(JmpE(m,p)*JmpE(q,n)+JmpE(p,m)*JmpE(n,q))+
     +         (2.0d0-2.0d0*nu)*(LnE(m)-JmpE(m,p))*(LnE(n)-JmpE(n,q)))

              KmnpqSA(mnindx,pqindx) = (aob**2)*(
     +               nu*(JmpE(m,p)*JmpO(q,n)+JmpE(p,m)*JmpO(n,q))+
     +         (2.0d0-2.0d0*nu)*(LnE(m)-JmpE(m,p))*(LnO(n)-JmpO(n,q)))

              KmnpqAS(mnindx,pqindx) = (aob**2)*(
     +               nu*(JmpO(m,p)*JmpE(q,n)+JmpO(p,m)*JmpE(n,q))+
     +         (2.0d0-2.0d0*nu)*(LnO(m)-JmpO(m,p))*(LnE(n)-JmpE(n,q)))

              KmnpqAA(mnindx,pqindx) = (aob**2)*(
     +               nu*(JmpO(m,p)*JmpO(q,n)+JmpO(p,m)*JmpO(n,q))+
     +         (2.0d0-2.0d0*nu)*(LnO(m)-JmpO(m,p))*(LnO(n)-JmpO(n,q)))

              KmnpqSS(mnindx,pqindx) = KmnpqSS(mnindx,pqindx)/
     +                                 (CnE(m)*CnE(n)*CnE(p)*CnE(q))
              KmnpqSA(mnindx,pqindx) = KmnpqSA(mnindx,pqindx)/
     +                                 (CnE(m)*CnO(n)*CnE(p)*CnO(q))
              KmnpqAS(mnindx,pqindx) = KmnpqAS(mnindx,pqindx)/
     +                                 (CnO(m)*CnE(n)*CnO(p)*CnE(q))
              KmnpqAA(mnindx,pqindx) = KmnpqAA(mnindx,pqindx)/
     +                                 (CnO(m)*CnO(n)*CnO(p)*CnO(q))

              if (mnindx.eq.pqindx) then

                KmnpqSS(mnindx,pqindx) = KmnpqSS(mnindx,pqindx)+
     +            (kn(2*m)**4+(aob*kn(2*n))**4)

                KmnpqSA(mnindx,pqindx) = KmnpqSA(mnindx,pqindx)+
     +            (kn(2*m)**4+(aob*kn(2*n+1))**4)

                KmnpqAS(mnindx,pqindx) = KmnpqAS(mnindx,pqindx)+
     +            (kn(2*m+1)**4+(aob*kn(2*n))**4)

                KmnpqAA(mnindx,pqindx) = KmnpqAA(mnindx,pqindx)+
     +            (kn(2*m+1)**4+(aob*kn(2*n+1))**4)

              end if

  40        continue
  50      continue
  60    continue
  70  continue
  
c**********************************************************************
c**** Call routine to find eigenvalues and eigenvectors            ****
c**********************************************************************

      matSize = (nTrunc+1)**2

      ifail = 0

      call F08FAF('V','U',matSize,KmnpqSS,maxMatSize,W,WORK,LWORK,ifail)

      do 90 m=0,nTrunc
        do 80 n=0,nTrunc

          mnindx = m*(nTrunc+1)+n+1

          if (mnindx.le.pTrunc) then

            lSS(mnindx) = W(mnindx)

            do 91 p=0,nTrunc
              do 81 q=0,nTrunc

                pqindx = p*(nTrunc+1)+q+1
                VSS(mnindx,p,q) = KmnpqSS(pqindx,mnindx)

  81          continue
  91        continue

          end if

  80    continue
  90  continue

      ifail = 0

      call F08FAF('V','U',matSize,KmnpqSA,maxMatSize,W,WORK,LWORK,ifail)

      do 110 m=0,nTrunc
        do 100 n=0,nTrunc

          mnindx = m*(nTrunc+1)+n+1

          if (mnindx.le.pTrunc) then

            lSA(mnindx) = W(mnindx)

            do 111 p=0,nTrunc
              do 101 q=0,nTrunc

                pqindx = p*(nTrunc+1)+q+1
                VSA(mnindx,p,q) = KmnpqSA(pqindx,mnindx)

 101          continue
 111        continue

          end if

 100    continue
 110  continue

      ifail = 0

      call F08FAF('V','U',matSize,KmnpqAS,maxMatSize,W,WORK,LWORK,ifail)

      do 130 m=0,nTrunc
        do 120 n=0,nTrunc

          mnindx = m*(nTrunc+1)+n+1

          if (mnindx.le.pTrunc) then

            lAS(mnindx) = W(mnindx)

            do 131 p=0,nTrunc
              do 121 q=0,nTrunc

                pqindx = p*(nTrunc+1)+q+1
                VAS(mnindx,p,q) = KmnpqAS(pqindx,mnindx)

 121          continue
 131        continue

          end if

 120    continue
 130  continue

      ifail = 0

      call F08FAF('V','U',matSize,KmnpqAA,maxMatSize,W,WORK,LWORK,ifail)

      do 150 m=0,nTrunc
        do 140 n=0,nTrunc

          mnindx = m*(nTrunc+1)+n+1

          if (mnindx.le.pTrunc) then

            lAA(mnindx) = W(mnindx)

            do 151 p=0,nTrunc
              do 141 q=0,nTrunc

                pqindx = p*(nTrunc+1)+q+1
                VAA(mnindx,p,q) = KmnpqAA(pqindx,mnindx)

 141          continue
 151        continue

          end if

 140    continue
 150  continue

c**********************************************************************
c**** Print values to screen                                       ****
c**********************************************************************
    
      do 180 p=1,pTrunc
        write (*,'(1x,i4,4(1x,f14.5))') p,lSS(p),lSA(p),lAS(p),lAA(p)
        write (fUnitRes,'(1x,i4,4(1x,f14.5))') 
     +                                  p,lSS(p),lSA(p),lAS(p),lAA(p)
 180  continue

      close(funitRes)

c**********************************************************************
c**** Calculate mode shapes                                        ****
c**********************************************************************

      do 270 xLoop=0,nXpts
        xoa = -1.0d0+2.0d0*dble(xLoop)/dble(nXpts)
        do 260 yLoop=0,nYpts
          yob = -1.0d0+2.0d0*dble(yLoop)/dble(nYpts)

          wnEx(0) = 0.5d0/CnE(0)
          wnOx(0) = 0.5d0*xoa/CnO(0)
          wnEy(0) = 0.5d0/CnE(0)
          wnOy(0) = 0.5d0*yob/CnO(0)

          do 220 n=1,nTrunc

            wnEx(n) = 0.25d0*((dexp(kn(2*n)*(xoa-1.0d0))+
     +                        dexp(-kn(2*n)*(xoa+1.0d0)))/
     +                        (1.0d0+dexp(-2.0d0*kn(2*n)))
     +       +dcos(kn(2*n)*xoa)/dcos(kn(2*n)))/CnE(n)
            
            wnOx(n) = 0.25d0*((dexp(kn(2*n+1)*(xoa-1.0d0))-
     +                         dexp(-kn(2*n+1)*(xoa+1.0d0)))/
     +                        (1.0d0-dexp(-2.0d0*kn(2*n+1)))
     +       +dsin(kn(2*n+1)*xoa)/dsin(kn(2*n+1)))/CnO(n)
            
            wnEy(n) = 0.25d0*((dexp(kn(2*n)*(yob-1.0d0))+
     +                         dexp(-kn(2*n)*(yob+1.0d0)))/
     +                        (1.0d0+dexp(-2.0d0*kn(2*n)))
     +       +dcos(kn(2*n)*yob)/dcos(kn(2*n)))/CnE(n)

            wnOy(n) = 0.25d0*((dexp(kn(2*n+1)*(yob-1.0d0))-
     +                         dexp(-kn(2*n+1)*(yob+1.0d0)))/
     +                         (1.0d0-dexp(-2.0d0*kn(2*n)))
     +       +dsin(kn(2*n+1)*yob)/dsin(kn(2*n+1)))/CnO(n)

 220      continue

          do 250 p=1,pTrunc
 
            wSS = 0.0d0
            wSA = 0.0d0
            wAS = 0.0d0
            wAA = 0.0d0

            do 240 m=0,nTrunc
              do 230 n=0,nTrunc

                wSS = wSS+VSS(p,m,n)*wnEx(m)*wnEy(n)
                wSA = wSA+VSA(p,m,n)*wnEx(m)*wnOy(n)
                wAS = wAS+VAS(p,m,n)*wnOx(m)*wnEy(n)
                wAA = wAA+VAA(p,m,n)*wnOx(m)*wnOy(n)

 230          continue
 240        continue

            write (fUnitRes+p,'(6(1x,f10.4))') xoa,yob/aob,
     +             wSS,wSA,wAS,wAA
          
 250      continue
            
 260    continue

        do 261 p=1,pTrunc
          write (fUnitRes+p,'('''')')
 261    continue

 270  continue

c**********************************************************************
c**** End program                                                  ****
c**********************************************************************

      end


      subroutine DATA_INPUT(fUnitRes,nTrunc,nu,aob,pTrunc)

c**********************************************************************
c**** Reads in data from datafiles                                 ****
c**********************************************************************
c                                                     column 72 here -->
      implicit none
      integer          maxNTrunc
      parameter        (maxNTrunc=48)
      integer          maxPTrunc
      parameter        (maxPTrunc=16)

      character*80     fNameInp,fNameRes,dummyChar

      integer          funitInp,funitRes
      integer          ind,LEN_TO_CHAR
      integer          nTrunc,pTrunc

      double precision aob,nu

c**********************************************************************
c**** Define the file unit numbers to be used                      ****
c**********************************************************************

      funitInp = 10

c**********************************************************************
c**** title and ask for control file name default available        ****
c**********************************************************************

      write (*,'(/18x,''EIGENMODES FOR RECTANGULAR ELASTIC PLATES'')')
      write (*,'(18x,''*****************************************'')')
      write (*,'(/''Enter input file (default <CR>) : ''$)')
      read (*,'(a)') fnameInp
      if (fnameInp.eq.'') then
        fnameInp = 'default.inp'
      endif

c**********************************************************************
c**** Open input file and read in all relevant data                ****
c**********************************************************************

      call FOPEN(funitInp,fnameInp,'r')

      write (*,'('' Reading nTrunc ... '')')
      read (funitInp,*) dummyChar
      read (funitInp,*) nTrunc
      write (*,'(1x,i4)') nTrunc
      if (nTrunc.gt.maxNTrunc) then
        write (*,'(/'' Error in DATA_INPUT: nTrunc > maxNTrunc'')')
        close(funitInp)
        stop
      endif

      write (*,'('' Reading nu ... '')')
      read (funitInp,*) dummyChar
      read (funitInp,*) nu
      write (*,'(1x,f10.6)') nu

      write (*,'('' Reading a/b, beta ... '')')
      read (funitInp,*) dummyChar
      read (funitInp,*) aob
      write (*,'(1x,f10.6)') aob

      write (*,'('' Reading pTrunc ... '')')
      read (funitInp,*) dummyChar
      read (funitInp,*) pTrunc
      write (*,'(1x,i4)') pTrunc
      if (pTrunc.gt.min(maxPTrunc,nTrunc)) then
        write (*,'(/'' Error in DATA_INPUT: pTrunc too big '')')
        close(funitInp)
        stop
      endif

c**********************************************************************
c**** Close control file                                           ****
c**********************************************************************

      close(funitInp)

c**********************************************************************
c**** Define a file for results to go into and open it for writing ****
c**********************************************************************

      ind = LEN_TO_CHAR(fnameInp,'.')
      fnameRes = fnameInp(1:ind)//'res'
      call FOPEN(funitRes,fnameRes,'w')

      return

      end

      subroutine FOPEN(funit,fname,rw_flag)

c**********************************************************************
c**** DESCRIPTION: Open a file for either reading or writing       ****
c**** USAGE:       On input:                                       ****
c****              funit          The file unit number             ****
c****               (integer)                                      ****
c****              fname          The file name                    ****
c****               (char*80)                                      ****
c****              rw_flag        Flag set to 'r' to read and      ****
c****               (char*1)      to 'w' to write                  ****
c****              On output:                                      ****
c**** DEPENDENCIES:LEN_TO_CHAR                                     ****
c**** AUTHORS:     R Porter                                        ****
c**** CHANGE_LOG:  Written 13/08/99 (tested 17/08/99)              ****
c**********************************************************************
c                                                     column 72 here -->
      implicit none

      character*80     fname
      character*1      rw_flag,ans

      logical          fexist

      integer          funit,LEN_TO_CHAR

      if (funit.lt.7) then
        write (*,'(/'' Error in FOPEN: funit less than 7'')')
        stop
      endif

 10   inquire(file=fname,exist=fexist)
      if (rw_flag.eq.'r') then
        if (fexist) then
          open(funit,file=fname,form='formatted',status='old')
          rewind(funit)
        else
          write (*,'(/'' Error in FOPEN: cannot open file '',a)')
     +               fname(1:LEN_TO_CHAR(fname,' '))
          write (*,'('' Do you want to try again? (y/n) : ''$)')
          read (*,'(a)') ans
          if (ans.eq.'y'.or.ans.eq.'Y') goto 10
          stop
        endif
      else if (rw_flag.eq.'w') then
        if (fexist) then
          write(*,'(/'' The file '',a,'' already exists.'')')
     +                   fname(1:LEN_TO_CHAR(fname,' '))
          write(*,'('' Do you want to overwrite it? (y/n) : ''$)')
          read(*,'(a)') ans
          if (ans.ne.'n'.or.ans.ne.'N') then
            open(funit,file=fname,form='formatted',status='old')
            rewind(funit)
          else
            write(*,'('' Do you want to choose a new'',
     +                '' filename? (y/n) :''$)')
            read (*,'(a)') ans
            if (ans.ne.'n'.or.ans.ne.'N') then
              write(*,'('' Enter new filename : ''$)')
              read (*,'(a)') fname
              goto 10
            endif
            stop
          endif
        else
          open(funit,file=fname,form='formatted',status='new')
        endif
      else
        write (*,'('' Error in FOPEN: rw_flag not "r" or "w"'')')
        pause
      endif

      return

      end

      integer function LEN_TO_CHAR(string,ch)

c**********************************************************************
c**** DESCRIPTION: Finds the length of a string to first occurence ****
c****              of a character                                  ****
c**** USAGE:       On input:                                       ****
c****              string         A string                         ****
c****               (char*(*))                                     ****
c****              ch             A character                      ****
c****               (char*1)                                       ****
c****              On output:                                      ****
c****              LEN_TO_CHAR:   Returns length of string to first****
c****               (integer)     occurrence of character ch       ****
c**** DEPENDENCIES:                                                ****
c**** AUTHORS:     R Porter                                        ****
c**** CHANGE_LOG:  Written 13/08/99                                ****
c**********************************************************************
c                                                     column 72 here -->
      implicit none

      character*(*)   string
      character*1     ch

      integer         i

      i = 1
      do 10 while (string(i:i).ne.ch.and.i.lt.len(string))
        i = i+1
 10   continue

      LEN_TO_CHAR = i

      return

      end

      subroutine PROOTS(nTrunc,kn)

c**********************************************************************
c**** Calculate roots of plate equation                            ****
c**********************************************************************

      implicit none
      include 'trunc.h'

      double precision eps,tol
      parameter        (eps=1.0d-02,tol=1.0d-12)

      integer          nTrunc,n

      double precision pi,sgn
      double precision kl,ku,km
      double precision fl,fu,fm
      double precision kn(0:2*maxNTrunc+1)

      common /const/   pi

c**********************************************************************
c**** First two are explicit                                       ****
c**********************************************************************

      kn(0) = 0.0d0
      kn(1) = 0.0d0

      sgn =-1.0d0

c**********************************************************************
c**** Loop over n=2 to n=2*nTrunc+1 roots                          ****
c**********************************************************************

      do 30 n=2,2*nTrunc+1

c**********************************************************************
c**** Set up upper and lower bounds                                ****
c**********************************************************************

        kl = 0.5d0*dble(n)*pi-0.25d0*pi-eps
        ku = 0.5d0*dble(n)*pi-0.25d0*pi+eps

        fl = dtanh(kl)-sgn*dtan(kl)
        fu = dtanh(ku)-sgn*dtan(ku)

        km = 0.5d0*(kl+ku)

        fm = dtanh(km)-sgn*dtan(km)

c**********************************************************************
c**** Bisection method to find roots                               ****
c**********************************************************************

        do while (0.5d0*(ku-kl).gt.tol.and.dabs(fm).gt.tol)

          if (fm*fl.lt.0.0d0) then
            ku = km
            fu = dtanh(ku)-sgn*dtan(ku)
          else
            kl = km
            fl = dtanh(kl)-sgn*dtan(kl)
          end if

          km = 0.5d0*(kl+ku)
          fm = dtanh(km)-sgn*dtan(km)

        end do 

c**********************************************************************
c**** Next time round and the equation for the roots changes       ****
c**********************************************************************

        sgn =-sgn
        kn(n) = km

  30  continue

      return
      end
