      program swtrunccyl 

c**********************************************************************
c**** Program to compute the scattering by waves by a metamaterial ****
c**** truncated cylinder                                           ****
c**********************************************************************

      implicit none

      integer          maxNTrunc,maxMatSize
      parameter        (maxNTrunc=64,maxMatSize=2*maxNTrunc+1)

      integer          funitRes,ifail,nz,info
      integer          xLoop,yLoop,khLoop
      integer          nXPts,nYPts,nkhPts
      integer          n,m,nTrunc,matSize,iOpt,i,j,iTrunc
      integer          pivot(maxMatSize)

      double complex   zi,zarg
      double precision pi,umn
      double precision ka,kh,kr,khStart,khStep
      double precision beta,aoh,doh,ddoh,mutma,mutmr
      double complex   zBesJVec(maxNTrunc+3)
      double complex   zBesYVec(maxNTrunc+3)
      double complex   CWRK(maxNTrunc+3)
      double precision t,h,gw(10),gp(10)
      double complex   bn(-maxNTrunc:maxNTrunc)
      double complex   an(-maxNTrunc:maxNTrunc)
      double precision bJka(-maxNTrunc:maxNTrunc)
      double precision bJdka(-maxNTrunc:maxNTrunc)
      double complex   bHka(-maxNTrunc:maxNTrunc)
      double complex   bHdka(-maxNTrunc:maxNTrunc)
      double precision bJmua(-maxNTrunc-1:maxNTrunc+1)
      double precision bJdmua(-maxNTrunc-1:maxNTrunc+1)
      double complex   Mntm,Kntm
      double complex   Mat(maxMatSize,maxMatSize)
      double complex   Rhs(maxMatSize)
      double precision xoaMin,xoaMax,yoaMin,yoaMax
      double precision xoa,yoa,roa,theta,sigma1,sigma2
      double complex   phi,Fy

c**********************************************************************
c**** Set constants                                                ****
c**********************************************************************

      zi = dcmplx(0.0d0,1.0d0)
      funitRes = 11
      pi = 4.0d0*datan(1.0d0)

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

      call DATA_INPUT(funitRes,iOpt,nTrunc,iTrunc,aoh,doh,beta,
     +                ddoh,khStart,khStep,nkhPts,
     +                xoaMin,xoaMax,nXPts,yoaMin,yoaMax,nYPts)

      beta = beta*pi/180.0d0

c**********************************************************************
c**** Gauss points and weights for integration                     ****
c**********************************************************************

      gp(1) =-9.7390652851717174d-01
      gp(2) =-8.6506336668898431d-01
      gp(3) =-6.7940956829902444d-01
      gp(4) =-4.3339539412924716d-01
      gp(5) =-1.4887433898163122d-01
      gp(6) = 1.4887433898163122d-01
      gp(7) = 4.3339539412924716d-01
      gp(8) = 6.7940956829902444d-01
      gp(9) = 8.6506336668898431d-01
      gp(10)= 9.7390652851717174d-01

      gw(1) = 6.6671344308689484d-02
      gw(2) = 1.4945134915057695d-01
      gw(3) = 2.1908636251598679d-01
      gw(4) = 2.6926671930999235d-01
      gw(5) = 2.9552422471475448d-01
      gw(6) = 2.9552422471475448d-01
      gw(7) = 2.6926671930999235d-01
      gw(8) = 2.1908636251598679d-01
      gw(9) = 1.4945134915057695d-01
      gw(10)= 6.6671344308689484d-02

c**********************************************************************
c**** Do loop over kh                                              ****
c**********************************************************************

      do 999 khLoop=1,nkhPts

        kh = khStart+dble(khLoop-1)*khStep
        ka = kh*aoh

c**********************************************************************
c**** Calculate and store J_n(ka), J_n'(ka), H_n(ka), H_n'(ka)     ****
c**********************************************************************

        zarg = dcmplx(ka,0.0d0)
        ifail = 0
        call S17DEF(0.0d0,zarg,nTrunc+2,'U',zBesJVec,nz,ifail)
        ifail = 0
        call S17DCF(0.0d0,zarg,nTrunc+2,'U',zBesYVec,nz,CWRK,ifail)

        bJka(0) = dble(zBesJVec(1))
        bJdka(0) =-dble(zBesJVec(2))

        bHka(0) = dcmplx(bJka(0),dble(zBesYVec(1)))
        bHdka(0) = dcmplx(bJdka(0),-dble(zBesYVec(2)))

        do 20 n=1,nTrunc

          umn = (-1.0d0)**n

          bJka(n) = dble(zBesJVec(n+1))
          bJka(-n) = umn*bJka(n)

          bJdka(n) = 0.5d0*(dble(zBesJVec(n))-dble(zBesJVec(n+2)))
          bJdka(-n) = umn*bJdka(n)

          bHka(n) = dcmplx(bJka(n),dble(zBesYVec(n+1)))
          bHka(-n) = umn*bHka(n)

          bHdka(n) = dcmplx(bJdka(n),
     +              0.5d0*(dble(zBesYVec(n))-dble(zBesYVec(n+2))))
          bHdka(-n) = umn*bHdka(n)

  20    continue

c**********************************************************************
c**** Initialise the matrix                                        ****
c**********************************************************************

        do 22 n=-nTrunc,nTrunc
          do 21 m=-nTrunc,nTrunc

            Mat(n+nTrunc+1,m+nTrunc+1) = dcmplx(0.0d0,0.0d0)

  21      continue
  22    continue

c**********************************************************************
c**** Loop through all the subintervals                            ****
c**** an in each subinterval loop through the Gauss points         ****
c**********************************************************************

        h = 2.0d0*pi/dble(2*iTrunc)

        do 60 j=0,iTrunc-1

          do 50 i=1,10

            t = h*(2.0d0*dble(j)+gp(i)+1.0d0)

c**********************************************************************
c**** Calculate and store M_n(t_m) and K_n(t_m)                    ****
c**********************************************************************

            mutma = ka/dsqrt(ddoh+(doh-ddoh)*(dsin(t)**2))

            zarg = dcmplx(mutma,0.0d0)
            ifail = 0
            call S17DEF(0.0d0,zarg,nTrunc+3,'U',zBesJVec,nz,ifail)

            bJmua(0) = dble(zBesJVec(1))
            bJdmua(0) =-dble(zBesJVec(2))

            do 30 n=1,nTrunc+1

              umn = (-1.0d0)**n

              bJmua(n) = dble(zBesJVec(n+1))
              bJmua(-n) = umn*bJmua(n)

              bJdmua(n) = 0.5d0*(dble(zBesJVec(n))-dble(zBesJVec(n+2)))
              bJdmua(-n) = umn*bJdmua(n)

  30        continue

            do 40 n=-nTrunc,nTrunc

              Mntm = bJmua(n)/bHka(n)

              Kntm = ((mutma/(ka*bHdka(n)))*(ddoh*bJdmua(n)
     +                    -0.5d0*(doh-ddoh)*zi*dsin(t)*
     +           (bJmua(n-1)*zexp(zi*t)+bJmua(n+1)*zexp(-zi*t))))

              do 35 m=-nTrunc,nTrunc

                Mat(n+nTrunc+1,m+nTrunc+1) = Mat(n+nTrunc+1,m+nTrunc+1)
     +                        +h*gw(i)*(Mntm-Kntm)*zexp(zi*dble(m-n)*t)

  35          continue      

  40        continue      

  50      continue

  60    continue

c**********************************************************************
c**** Right hand side                                              ****
c**********************************************************************

        do 80 n=-nTrunc,nTrunc

          Rhs(n+nTrunc+1) = 2.0d0*zi*zexp(-zi*dble(n)*beta)
     +                      /(pi*ka*bHka(n)*bHdka(n))

  80    continue

c**********************************************************************
c**** Invert system of equations                                   ****
c**********************************************************************

        matSize = 2*nTrunc+1

        call CLUDCMP(Mat,matSize,maxMatSize,pivot,info)
        call CLUBKSB(Mat,matSize,maxMatSize,pivot,Rhs)

c**********************************************************************
c**** Hive off the coefficients                                    ****
c**********************************************************************

        do 90 n=-nTrunc,nTrunc

          bn(n) = Rhs(n+nTrunc+1)
c         print *,bn(n)

  90    continue

c**********************************************************************
c**** And calculate the a_n's                                      ****
c**********************************************************************

        do 100 n=-nTrunc,nTrunc

          an(n) =-bJka(n)*zexp(-zi*dble(n)*beta)/bHka(n)
c         an(n) =-bJdka(n)*zexp(-zi*dble(n)*beta)/bHdka(n)

 100    continue

c**********************************************************************
c**** Loop through all the subintervals                            ****
c**** an in each subinterval loop through the Gauss points         ****
c**********************************************************************

        h = 2.0d0*pi/dble(2*iTrunc)

        do 160 j=0,iTrunc-1

          do 150 i=1,10

            t = h*(2.0d0*dble(j)+gp(i)+1.0d0)

c**********************************************************************
c**** Calculate and store M_n(t_m) and K_n(t_m)                    ****
c**********************************************************************

            mutma = ka/dsqrt(ddoh+(doh-ddoh)*(dsin(t)**2))

            zarg = dcmplx(mutma,0.0d0)
            ifail = 0
            call S17DEF(0.0d0,zarg,nTrunc+3,'U',zBesJVec,nz,ifail)

            bJmua(0) = dble(zBesJVec(1))
            bJdmua(0) =-dble(zBesJVec(2))

            do 130 n=1,nTrunc+1

              umn = (-1.0d0)**n

              bJmua(n) = dble(zBesJVec(n+1))
              bJmua(-n) = umn*bJmua(n)

              bJdmua(n) = 0.5d0*(dble(zBesJVec(n))-dble(zBesJVec(n+2)))
              bJdmua(-n) = umn*bJdmua(n)

 130        continue

            do 140 n=-nTrunc,nTrunc

              Mntm = (bJmua(n)/bHka(n))

              Kntm = ((mutma/(ka*bHdka(n)))*(ddoh*bJdmua(n)
     +                    -0.5d0*(doh-ddoh)*zi*dsin(t)*
     +           (bJmua(n-1)*zexp(zi*t)+bJmua(n+1)*zexp(-zi*t))))

              do 135 m=-nTrunc,nTrunc

                an(n) = an(n)+bn(m)*h*gw(i)*Mntm*zexp(zi*dble(m-n)*t)
c               an(n) = an(n)+bn(m)*h*gw(i)*Kntm*zexp(zi*dble(m-n)*t)

 135          continue      

 140        continue      

 150      continue

 160    continue

        do 161 n=-nTrunc,nTrunc

c         print *,n,an(n)

 161    continue

c**********************************************************************
c**** Plot the surface elevation (if nkhPts = 1)                   ****
c**********************************************************************

        if (nkhPts.eq.1) then 

          do 220 xLoop=0,nXPts

            xoa = xoaMin+(xoaMax-xoaMin)*xLoop/nXPts

            do 210 yLoop=0,nYPts

              yoa = yoaMin+(yoaMax-yoaMin)*yLoop/nYPts

              roa = dsqrt(xoa**2+yoa**2)
              theta = datan2(yoa,xoa)

              kr = ka*roa

              if (roa.ge.1.0d0) then 

c**********************************************************************
c**** Outside the metamaterial cylinder                            ****
c**********************************************************************

                zarg = dcmplx(kr,0.0d0)
                ifail = 0
                call S17DEF(0.0d0,zarg,nTrunc+1,'U',zBesJVec,nz,ifail)
                ifail = 0
                call S17DCF(0.0d0,zarg,nTrunc+1,'U',zBesYVec,nz,CWRK,
     +                      ifail)

                bJka(0) = dble(zBesJVec(1))
                bHka(0) = dcmplx(bJka(0),dble(zBesYVec(1)))

                do 170 n=1,nTrunc

                  umn = (-1.0d0)**n
                  bJka(n) = dble(zBesJVec(n+1))
                  bJka(-n) = umn*bJka(n)

                  bHka(n) = dcmplx(bJka(n),dble(zBesYVec(n+1)))
                  bHka(-n) = umn*bHka(n)

 170            continue

                phi = zexp(zi*kr*dcos(theta-beta))

                do 180 n=-nTrunc,nTrunc

                  phi = phi+(zi**n)*an(n)*bHka(n)*zexp(zi*dble(n)*theta)

 180            continue

              else 

c**********************************************************************
c**** Inside the metamaterial cylinder                             ****
c**********************************************************************

                phi = dcmplx(0.0d0,0.0d0)

                h = 2.0d0*pi/dble(2*iTrunc)

                do 260 j=0,iTrunc-1

                  do 250 i=1,10

                    t = h*(2.0d0*dble(j)+gp(i)+1.0d0)

                    mutmr = kr/dsqrt(ddoh+(doh-ddoh)*(dsin(t)**2))

                    do 240 n=-nTrunc,nTrunc

                        phi = phi+bn(n)*h*gw(i)*
     +                        zexp(zi*(mutmr*dcos(theta-t)+dble(n)*t))

  240               continue

  250             continue
                
  260           continue

              end if

c**********************************************************************
c**** Print x,y,phi to file                                        ****
c**********************************************************************
        
              write (fUnitRes,'(4(1x,f12.6))') xoa,yoa,
     +                                         dble(phi),dimag(phi)

  210       continue

            write (fUnitRes,'('''')')

  220     continue

        else

c**********************************************************************
c**** Compute the scattering cross section                         ****
c**********************************************************************
        
          sigma1 = 0.0d0
          sigma2 = 0.0d0

          do 360 n=-nTrunc,nTrunc

            sigma1 = sigma1+cdabs(an(n))**2
            sigma2 = sigma2-dble(an(n)*zexp(zi*dble(n)*beta))

 360      continue

c**********************************************************************
c**** Compute the vertical force                                   ****
c**********************************************************************
        
          Fy = dcmplx(0.0d0,0.0d0)

          h = 2.0d0*pi/dble(2*iTrunc)

          do 430 j=0,iTrunc-1

            do 420 i=1,10

              t = h*(2.0d0*dble(j)+gp(i)+1.0d0)

              mutma = ka/dsqrt(ddoh+(doh-ddoh)*(dsin(t)**2))

c**********************************************************************
c**** We need J_1(\mu(t)a)                                         ****
c**********************************************************************

              zarg = dcmplx(mutma,0.0d0)
              ifail = 0
              call S17DEF(0.0d0,zarg,2,'U',zBesJVec,nz,ifail)

              bJmua(1) = dble(zBesJVec(2))

              do 410 n=-nTrunc,nTrunc

                Fy = Fy+gw(i)*h*bn(n)*dsin(t)*zexp(zi*dble(n)*t)
     +                               *bJmua(1)

 410          continue

 420        continue
 430      continue

c**********************************************************************
c**** Normalise the force                                          ****
c**********************************************************************

          if (dabs(beta).gt.1.0d-10) then
            Fy = 0.5d0*Fy*zi*pi*ka*(ddoh-doh)*bHdka(1)/dsin(beta)
          else
            Fy = dcmplx(0.0d0,0.0d0)
          end if
        
c**********************************************************************
c**** Print Fy to file                                             ****
c**********************************************************************
        
          write (*,'(4(1x,f12.6))') kh,cdabs(Fy),sigma1,sigma2
          write (fUnitRes,'(3(1x,f12.6))') kh,cdabs(Fy),sigma1

        end if

 999  continue

c**********************************************************************
c**** End                                                          ****
c**********************************************************************

      close(funitRes)

      end

      subroutine cludcmp(a,n,np,indx,d)
c**********************************************************************
c**** Perform LU decomposition on complex matrix A                 ****
c**** Input:  a(np,np)                                             ****
c****         np is physical dimension of A                        ****
c****         n actual dimension                                   ****
c**** Output: indx(n) is a vector recording row permutations       ****
c****         d = +1/-1  even/odd interchanges                     ****
c**********************************************************************
      implicit none

      integer          d,i,j,k,imax,n,np,indx(np)
      double complex   a(np,np)
      double precision vv(np)
      double complex   sum
      double precision aamax,tiny,dum
      intrinsic        cdabs

      tiny=1.0d-20
      d = 1
      do 12 i=1,n
        aamax=0.0d0
        do 11 j=1,n
          if (cdabs(a(i,j)).gt.aamax) aamax = cdabs(a(i,j))
  11    continue
        if (aamax.eq.0.0d0) pause 'Singular matrix'
        vv(i) = 1.0d0/aamax
  12  continue
      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum = a(i,j)
            if (i.gt.1) then
              do 13 k=1,i-1
                sum = sum-a(i,k)*a(k,j)
  13          continue
              a(i,j) = sum
             end if
  14       continue
         end if
         aamax = 0.0d0
         do 16 i=j,n
           sum = a(i,j)
           if (j.gt.1) then
             do 15 k=1,j-1
               sum = sum-a(i,k)*a(k,j)
  15         continue
             a(i,j) = sum
           end if
           dum = vv(i)*cdabs(sum)
           if (dum.ge.aamax) then
             imax = i
             aamax = dum
           endif
  16     continue
         if (j.ne.imax) then
           do 17 k=1,n
             sum = a(imax,k)
             a(imax,k) = a(j,k)
             a(j,k) = sum
  17       continue
           d = -d
           vv(imax) = vv(j)
         endif
         indx(j) = imax
         if (j.ne.n) then
           if (cdabs(a(j,j)).eq.0.0d0) a(j,j) = tiny
           sum = 1.0d0/a(j,j)
           do 18 i=j+1,n
             a(i,j) = a(i,j)*sum
  18       continue
         endif
  19   continue
       if (cdabs(a(n,n)).eq.0.0d0) a(n,n) = tiny
       return
       end

      subroutine clubksb(a,n,np,indx,b)
c**********************************************************************
c**** Perform forward and backward substitution for LU             ****
c**** decomposition of a where a is complex                        ****
c**** Input:  a(np,np) in LU decomposition form from CLUDCMP       ****
c****         np phyiscal dimension                                ****
c****         n actual dimension                                   ****
c****         indx(n) contains all permutations                    ****
c****         b(n) contains RHS of a.x=b                           ****
c**** output: b(n) overwrite RHS of a.x=b with x                   ****
c**********************************************************************
      implicit none

      integer        n,np,indx(np)
      double complex a(np,np),b(np)
      integer        i,ii,j,ll
      double complex sum

      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0) then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
  11      continue
        else if (cdabs(sum).ne.0.0d0) then
          ii=i
        end if
        b(i)=sum
  12  continue
      do 14 i=n,1,-1
        sum=b(i)
        if (i.lt.n) then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
  13      continue
        end if
        b(i)=sum/a(i,i)
  14  continue
      return

      end
