c
c plank.f - Solves for the eigenfrequencies and eigenmodes
c           of a flexible reed with a rectangular cross section
c           suspended from one end at a fixed angle and position,
c           with the other end free.  The solution is obtained
c           for small oscillations about the equilibrium (flat).
c
c author: richard.t.jones at uconn.edu
c version: April 17, 2012
c
c Model parameters include:
c	1) L = length of the reed (m)
c	2) T(x) = thickness of the reed (m) for 0<x<L
c	3) W(x) = width of the reed (m) for 0<x<L
c	4) E(x) = Young's elastic modulus for the reed (Pa)
c	5) rho(x) = density of the reed material (kg/m^3)
c
      subroutine plank()
      include 'plank.inc'
      data L/5e-3/
      data Coef(1,1)/0/
      pi_=datan2(0.d0,-1.d0)
      end

      real function T(x)
      real x
      include 'plank.inc'
c      T = 3e-4
      T = 2e-5
      end

      real function W(x)
      real x
      include 'plank.inc'
c      W = 4e-3
c      W = 4e-3*abs(1-x/L)
      W = 4e-3*(1.01-abs(1-2*x/L))
      w = min(w,2e-3)
c      W = 4e-3*sqrt(1-(2*x/L-1)**2)
      end

      real function E(x)
      real x
      include 'plank.inc'
      E = 1.22e12
      end

      real function rho(x)
      real x
      include 'plank.inc'
      rho = 3.52e3
      end

      double precision function rootkL(ord)
      integer ord
      include 'plank.inc'
      double precision kL,kLlast
      double precision arga
      integer osign,org
      if (ord/2*2.eq.ord) then
        org = ord-1
        osign = 1
      else
        org = ord
        osign = -1
      endif
      kLlast=-1
      kL=0
      do while (abs(kL-kLlast).gt.1e-17*kLlast)
        kLlast=kL
        kL=pi_*org+osign*arga(kL)
      enddo
      rootkL=kL
      print *, 'root',ord,' found to be',rootkL
      end

      subroutine basis_solve(ord)
      integer ord
      include 'plank.inc'
c
c  Uniform plank solution is of the general form, where x = 0..L : 
c
c               ikx     -ikx     kx      -kx
c     f(x) = f e   + f e    + f e   + f e
c             1       2        3       4
c
c  Normalization chosen below is f(L) = 1.  The fact that f(x) is real
c  constrains complex coefficient f_2 = conjg(f_1) and f_3 to pure real.
c  Below I represent f_1 = f_1R + i f_1I and f_2 = f_1R - i f_1I, where
c  f_1R and f_1I are pure real.  The boundary conditions to be enforced
c  are written as follows, in terms of the coefficients with f_4 = 1.
c
c     2 f_1R + f_3 = -1
c    -2 f_1I + f_3 = +1
c    -2 f_1R aR + 2 f_1I aI + f_3 b = -1/b
c     2 f_1I aR + 2 f_1R aI + f_3 b = +1/b
c
c  Only three of these equations are independent, the fourth being
c  guaranteed by the solution to the eigenvalue equation provided in
c  the rootkL() function.
c
      double precision C(3,3),v(3)
      integer work(3)
      double precision aR,aI,b
      double precision kL,fL
      double precision rootkL
      kL=rootkL(ord)
      ksol=kL/L
      aR=dcos(kL)
      aI=dsin(kL)
      b=dexp(kL)
      C(1,1)=2
      C(1,2)=0
      C(1,3)=1
      v(1)=-1
      C(2,1)=0
      C(2,2)=-2
      C(2,3)=1
      v(2)=1
      C(3,1)=-2*aR
      C(3,2)=+2*aI
      C(3,3)=b
      v(3)=-1/b
      call DEQN(3,C,3,work,ifail,1,v)
      if (ifail.ne.0) then
        print *, 'REQN returned',ifail,', solution failed!'
        f1R=0
        f1I=0
        f3=0
        f4=0
        return
      endif
      f1R=v(1)
      f1I=v(2)
      f3=v(3)
      f4=1
c
c  Renormalize the result to give f(L)=1
c
      fL=2*f1R*aR-2*f1I*aI+f3*b+1/b
      f1R=f1R/fL
      f1I=f1I/fL
      f3=f3/fL
      f4=1/fL
      print *, 'Solution for basis mode',ord,' follows:'
      print *, '  f1R = ',f1R
      print *, '  f1I = ',f1I
      print *, '  f3 = ',f3
      print *, '  f4 = ',f4
      print *, 'Use function fsol(x) to plot it.'
      print *, 'Test boundary condition 1:',2*f1R+f3+f4
      print *, 'Test boundary condition 2:',-2*f1I+f3-f4
      print *, 'Test boundary condition 3:',-2*f1R*aR+2*f1I*aI+f3*b+f4/b
      print *, 'Test boundary condition 4:',2*f1I*aR+2*f1R*aI+f3*b-f4/b
      end

      real function fsol(x)
      real x
      include 'plank.inc'
      double precision kx
      if (x.ge.0 .and. x.le.L) then
         kx=ksol*x
         fsol=2*f1R*cos(kx)-2*f1I*sin(kx)+f3*exp(kx)+f4*exp(-kx)
      else
         fsol=0
      endif
      end

      real function fp(x)
      real x
      include 'plank.inc'
      double precision kx
      if (x.ge.0 .and. x.le.L) then
         kx=ksol*x
         fp=-2*f1R*sin(kx)-2*f1I*cos(kx)+f3*exp(kx)-f4*exp(-kx)
      else
         fp=0
      endif
      end

      real function fpp(x)
      real x
      include 'plank.inc'
      double precision kx
      if (x.ge.0 .and. x.le.L) then
         kx=ksol*x
         fpp=-2*f1R*cos(kx)+2*f1I*sin(kx)+f3*exp(kx)+f4*exp(-kx)
      else
         fpp=0
      endif
      end

      real function fppp(x)
      real x
      include 'plank.inc'
      double precision kx
      if (x.ge.0 .and. x.le.L) then
         kx=ksol*x
         fppp=+2*f1R*sin(kx)+2*f1I*cos(kx)+f3*exp(kx)-f4*exp(-kx)
      else
         fppp=0
      endif
      end

      double precision function arga(kL)
      double precision kL
      include 'plank.inc'
      double precision b
      b=dexp(kL)
      arga=datan2(b**2-1,2*b)
      end

      subroutine eigensolve
      include 'plank.inc'
c
c  Solution is of the general form, where x = 0..L : 
c
c             _M_  
c     y (x) = \   C   f (x)
c      n      /__  mn  m
c             m=1
c
c  where f_m is basis function m from the set of M basis functions.
c  The basis functions form an orthogonal set, but are not normalized
c  to unity.  The full model is given by the following diff. equation.
c
c                     2   /       2  \
c      /    E    \   d   |    3  d    |                2
c     |  -------  | ---- | W T  ----  | y (x)  =  omega   y (x)
c      \ rho W T /     2 |         2  |  n             n   n
c                    dx   \      dx  /
c
c  where n is the mode index.  The basis functions are the solutions
c  to this differential equation for the case that W,T,rho,E are all
c  independent of x.  In the more general problem, where all of these
c  parameters can depend on x, is solved below.
c
      integer m,n,ord
      double precision sum1,sumH0,sumH1,sumH2
      double precision fn(Xpoints)
      double precision WT3(3)
      real x,dx

      parameter (lwork = 500*Modes)
      double precision work(lwork)
      integer info

      real fsol,fp,fpp,fppp
      external fsol,fp,fpp,fppp
      real rho,W,T,E
      external rho,W,T,E

      do ord=1,Modes
        call basis_solve(ord)
        kmode(ord) = ksol
        f1Rmode(ord) = f1R
        f1Imode(ord) = f1I
        f3mode(ord) = f3
        f4mode(ord) = f4
      enddo

c  Use the fact the H is symmetric to reduce the total number
c  of matrix elements that must be separately computed.
      dx = L/Xpoints
      do n=1,Modes
        ksol = kmode(n)
        f1R = f1Rmode(n)
        f1I = f1Imode(n)
        f3 = f3mode(n)
        f4 = f4mode(n)
        Norm(n) = 0
        do i=1,Xpoints-1
          x = i*dx
          fn(i) = fsol(x)
          Norm(n) = Norm(n)+fn(i)**2
        enddo
        Norm(n) = (Norm(n)+0.5)*dx
        do m=1,Modes
          ksol = kmode(m)
          f1R = f1Rmode(m)
          f1I = f1Imode(m)
          f3 = f3mode(m)
          f4 = f4mode(m)
          WT3(2) = W(0)*T(0)**3
          WT3(3) = W(dx)*T(dx)**3
          sumH0 = 0
          sumH1 = 0
          sumH2 = 0
          do i=1,Xpoints-1
            x = i*dx
            WT3(1) = WT3(2)
            WT3(2) = WT3(3)
            WT3(3) = W(x+dx)*T(x+dx)**3
            dWT3dx = (WT3(3)-WT3(1))/(2*dx)
            d2WT3dx2 = (WT3(3)+WT3(1)-2*WT3(2))/dx**2
            Efact = E(x)/(rho(x)*W(x)*T(x))
            sumH0 = sumH0+fn(i)*fsol(x)*Efact*WT3(2)*ksol**4
            sumH1 = sumH1+fn(i)*fppp(x)*Efact*dWT3dx*ksol**3
            sumH2 = sumH2+fn(i)*fpp(x)*Efact*d2WT3dx2*ksol**2
          enddo
          x = L
          sumH0 = sumH0+0.5*(E(x)/rho(x))*(T(x)*ksol**2)**2
          sumH0 = sumH0*dx;
          sumH1 = sumH1*dx;
          sumH2 = sumH2*dx;
c         print 1000, n,m,sumH0,sumH1,sumH2,Norm(n)
c1000     format(i2,',',i2,' : ',4(e12.6,' '))
          Coef(n,m) = (sumH0+2*sumH1+sumH2)/Norm(n)
          if (n.eq.m) then
            Unit(n,m) = 1
          else
            Unit(n,m) = 0
          endif
        enddo
      enddo

      call dggev('V','V',Modes,Coef,Modes,Unit,Modes,
     +           alphaR,alphaI,beta,Vl,Modes,Vr,Modes,
     +           work,lwork,info)
      if (info.lt.0) then
        print *, 'Error from DSYEV: problem with argument',-info
        return
      elseif (info.gt.0) then
        print *, 'Error from DSYEV: convergence failure, info =',info
        return
      elseif (work(1).gt.lwork) then
        print *, 'Warning from DSYEV: algorithm wants lwork =',work(1),
     +           ', but only received',lwork
        return
      endif

      do ord=1,Modes
        if (alphaI(ord).ne.0) then
          print *, 'Error - complex eigenvalue found in position',ord,
     +             ' -- cannot continue.'
          return
       elseif (alphaR(ord).lt.0) then
         print *, 'Error - negative eigenvalue found in position',ord,
     +             ' -- cannot continue.'
         return
        else
          evalue(ord) = alphaR(ord)/beta(ord)
        endif
      enddo
      call sortzv(evalue,iorder,Modes,1,0,0)

      print *, 'Eigenfrequencies are as follows:'
      do ord=1,Modes
        print 2000, ord,sqrt(E(0)/rho(0))*kmode(ord)**2*T(0)/(2*pi_),
     +              sqrt(evalue(iorder(ord)))/(2*pi_)
 2000   format(' mode',i3,': ',e15.8,'  ',e15.8)
      enddo
      end

      real function ysol(x,n)
      real x
      integer n
      include 'plank.inc'
      double precision ysum
      integer ord
      if (Coef(1,1).eq.0) then
        call eigensolve
      endif
      ysol = 0
      ysum = 0
      do ord=1,Modes
        ksol = kmode(ord)
        f1R = f1Rmode(ord)
        f1I = f1Imode(ord)
        f3 = f3mode(ord)
        f4 = f4mode(ord)
        ysol = ysol+Vr(ord,iorder(n))*fsol(x)
        ysum = ysum+Vr(ord,iorder(n))
      enddo
      ysol = ysol/ysum
      end