      program flttrunccyl 

c**********************************************************************
c**** Program to compute the scattering by waves by a metamaterial ****
c**** truncated cylinder                                           ****
c**********************************************************************
      
      implicit none
      integer          maxNTrunc,maxQTrunc,maxMatSize
      parameter        (maxNTrunc=16,maxQTrunc=128,
     +                  maxMatSize=(maxQTrunc+1)*(2*maxNTrunc+1))

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

      double complex   zi,zarg
      double precision pi,umn
      double precision beta,aoh,doh,hod,aod,ddoh
      double precision ka,kh,kd,kr,khStart,khStep
      double precision bigKh,bigKd
      double precision kmh(maxQTrunc),kmd(maxQTrunc)
      double precision N0sqrt,NmSqrt(maxQTrunc)
      double precision t,h,gw(10),gp(10)
      double precision mud,muqd(maxQTrunc)
      double complex   zB1Vec(maxNTrunc+3),zB2Vec(maxNTrunc+3)
      double complex   CWRK(maxNTrunc+3)
      double precision bJka(-maxNTrunc:maxNTrunc) 
      double precision bJdka(-maxNTrunc:maxNTrunc)
      double complex   bHka(-maxNTrunc:maxNTrunc)
      double complex   bHdka(-maxNTrunc:maxNTrunc)
      double precision bKkma(-maxNTrunc:maxNTrunc,maxQTrunc)
      double precision bKdkma(-maxNTrunc:maxNTrunc,maxQTrunc)
      double complex   bJmuqa(-maxNTrunc-1:maxNTrunc+1,0:maxQTrunc)
      double complex   bJdmuqa(-maxNTrunc-1:maxNTrunc+1,0:maxQTrunc)
      double complex   FmqP(0:maxQTrunc,0:maxQTrunc)
      double complex   FmqM(0:maxQTrunc,0:maxQTrunc)
      double complex   Mpqnm
      double complex   Kpqnm
      double complex   Mat(maxMatSize,maxMatSize)
      double complex   Rhs(maxMatSize)
      double complex   bpq(-maxNTrunc:maxNTrunc,0:maxQTrunc)
      double complex   anm(-maxNTrunc:maxNTrunc,0:maxQTrunc)
      double precision xoaMin,xoaMax,yoaMin,yoaMax
      double precision xoa,yoa,roa,rod,theta,psi0,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,qTrunc,nTrunc,iTrunc,aoh,doh,beta,
     +                ddoh,khStart,khStep,nkhPts,
     +                xoaMin,xoaMax,nXPts,yoaMin,yoaMax,nYPts)

      matSize = (qTrunc+1)*(2*nTrunc+1)

      beta = beta*pi/180.0d0
      hod = 1.0d0/doh
      aod = aoh*hod

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**********************************************************************

c Either loop through kh or 
      do 999 khLoop=1,nkhPts

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

c comment out 2 lines above and uncomment three below to run over angles
c     do 999 s=0,90
c       beta = dble(s)*pi/180.0d0
c       kh = khStart

        ka = kh*aoh
        kd = kh*doh

        bigKh = kh*dtanh(kh)
        bigKd = bigKh*doh

        N0sqrt = dsqrt(0.5d0*(1.0d0+dsinh(2.0d0*kh)/(2.0d0*kh)))

c Auto control of truncation 
c       nTrunc = int(ka*2)+2
c       print *,nTrunc
c       qTrunc = int(12.0/(ka+0.1d0))+4
c       print *,qTrunc
c       matSize = (qTrunc+1)*(2*nTrunc+1)
c       iTrunc = 30*qTrunc
c       print *,iTrunc

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

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

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

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

        do 10 n=1,nTrunc

          umn = (-1.0d0)**n

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

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

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

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

  10    continue

c**********************************************************************
c**** Calculate the eigenvalues and normalising factors            ****
c**********************************************************************

        call EIGVAL(bigKh,kmh,qTrunc,maxQTrunc)

        do 30 m=1,qTrunc

          kmd(m) = kmh(m)*doh

          NmSqrt(m) = 
     +        dsqrt(0.5d0*(1.0d0+dsin(2.0d0*kmh(m))/(2.0d0*kmh(m))))

c**********************************************************************
c**** Calculate and store K_n(k_m a), K_n'(k_m a) SCALED           ****
c**********************************************************************

          zarg = dcmplx(kmh(m)*aoh,0.0d0)
          ifail = 0
          call S18DCF(0.0d0,zarg,nTrunc+2,'S',zB2Vec,nz,ifail)

          bKkma(0,m) = dble(zB2Vec(1))
          bKdkma(0,m) =-dble(zB2Vec(2))

          do 20 n=1,nTrunc

            bKkma(n,m) = dble(zB2Vec(n+1))
            bKkma(-n,m) = bKkma(n,m)

            bKdkma(n,m) =-0.5d0*(dble(zB2Vec(n))+dble(zB2Vec(n+2)))
            bKdkma(-n,m) = bKdkma(n,m)

  20      continue

  30    continue

c**********************************************************************
c**** Set up the matrix                                            ****
c**********************************************************************

        do 50 inm=1,matSize
          do 40 ipq=1,matSize

            Mat(inm,ipq) = dcmplx(0.0d0,0.0d0)

  40      continue
  50    continue

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

        h = 0.5d0*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**** Determine the real root of the \mu-dispersion relation       ****
c**** and calculate the Bessel functions J_n(\mu_0 a)              ****
c**********************************************************************

            call MUDISPER(mud,t,hod,bigKd)

            zarg = dcmplx(mud*aod,0.0d0)
            ifail = 0
            call S17DEF(0.0d0,zarg,nTrunc+3,'U',zB1Vec,nz,ifail)

            bJmuqa(0,0) = dble(zB1Vec(1))
            bJdmuqa(0,0) =-dble(zB1Vec(2))

            do 60 n=1,nTrunc+1

              umn = (-1.0d0)**n

              bJmuqa(n,0) = dble(zB1Vec(n+1))
              bJmuqa(-n,0) = umn*bJmuqa(n,0)

              bJdmuqa(n,0) = 0.5d0*(dble(zB1Vec(n))-dble(zB1Vec(n+2)))
              bJdmuqa(-n,0) = umn*bJdmuqa(n,0)

  60        continue

c**********************************************************************
c**** Determine the imag roots of the \mu-dispersion relation      ****
c**** and calculate the Bessel functions I_n(\mu_q a)              ****
c**** Either use BesselI calls and fiddle with prefactors or       ****
c**** use BesselJ functions. In both cases, SCALED                 ****
c**********************************************************************

            call MUEIGVAL(muqd,qTrunc,t,hod,bigKd)
c if you want to output variation of mu_q(t)
c           write (fUnitRes,'(11(1x,f12.6))') t,(muqd(q),q=1,qTrunc)

            do 56 q=1,qTrunc
 
              zarg = dcmplx(muqd(q)*aod,0.0d0)
              ifail = 0
              call S18DEF(0.0d0,zarg,nTrunc+3,'S',zB1Vec,nz,ifail)

              bJmuqa(0,q) = dble(zB1Vec(1))
              bJdmuqa(0,q) =-zi*dble(zB1Vec(2))

              do 57 n=1,nTrunc+1

                umn = (-1.0d0)**n

                bJmuqa(n,q) = (zi**n)*dble(zB1Vec(n+1))
                bJmuqa(-n,q) = umn*bJmuqa(n,q) 

                bJdmuqa(n,q) = (zi**(n-1))* 
     +                         0.5d0*(dble(zB1Vec(n))+dble(zB1Vec(n+2)))
                bJdmuqa(-n,q) = umn*bJdmuqa(n,q) 

  57          continue

  56        continue

c Alternative computation of J_n(\mu_q(t)a)
c           do 65 q=1,qTrunc
c
c             zarg = dcmplx(0.0d0,muqd(q)*aod)
c             ifail = 0
c             call S17DEF(0.0d0,zarg,nTrunc+3,'S',zB1Vec,nz,ifail)

c             bJmuqa(0,q) = zB1Vec(1)
c             bJdmuqa(0,q) =-zB1Vec(2)

c             do 70 n=1,nTrunc+1

c               umn = (-1.0d0)**n

c               bJmuqa(n,q) = zB1Vec(n+1)
c               bJmuqa(-n,q) = umn*bJmuqa(n,q)

c               bJdmuqa(n,q) = 0.5d0*(zB1Vec(n)-zB1Vec(n+2))
c               bJdmuqa(-n,q) = umn*bJdmuqa(n,q)

c 70          continue

c 65        continue

c**********************************************************************
c**** Now compute F_{m0}^+, F_{m0}^-                               ****
c**********************************************************************

            FmqP(0,0) = (1.0d0/(N0sqrt*(kd**2-mud**2)))*
     +                (mud*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (mud*dcosh(mud)-bigKd*dsinh(mud)))

            FmqM(0,0) =-(1.0d0/(N0sqrt*(kd**2-(mud*dcos(t))**2)))*
     +                (mud*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (mud*dcosh(mud)-bigKd*dsinh(mud)))

            do 80 m=1,QTrunc
     
              FmqP(m,0) =-(1.0d0/(Nmsqrt(m)*(kmd(m)**2+mud**2)))*
     +                  (mud*dcos(kmd(m)*(hod-1.0d0))*
     +                   (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                   +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                   (mud*dcosh(mud)-bigKd*dsinh(mud)))

              FmqM(m,0) = 
     +              (1.0d0/(Nmsqrt(m)*(kmd(m)**2+(mud*dcos(t))**2)))
     +                 *(mud*dcos(kmd(m)*(hod-1.0d0))*
     +                   (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                   +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                   (mud*dcosh(mud)-bigKd*dsinh(mud)))

  80        continue

c**********************************************************************
c**** Now compute F_{mq}^+, F_{mq}^-                               ****
c**********************************************************************

            do 100 q=1,qTrunc

              FmqP(0,q) = (zi/(N0sqrt*(kd**2+muqd(q)**2)))*
     +                (muqd(q)*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

              FmqM(0,q) =-(zi/(N0sqrt*(kd**2+(muqd(q)*dcos(t))**2)))*
     +                (muqd(q)*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

              do 90 m=1,qTrunc

                FmqP(m,q) =-(zi/(Nmsqrt(m)*(kmd(m)**2-muqd(q)**2)))*
     +                 (muqd(q)*dcos(kmd(m)*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

                FmqM(m,q) = 
     +                (zi/(Nmsqrt(m)*(kmd(m)**2-(muqd(q)*dcos(t))**2)))*
     +                  (muqd(q)*dcos(kmd(m)*(hod-1.0d0))*
     +                  (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                  +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                  (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

  90          continue

 100        continue

c**********************************************************************
c**** Matrix elements put together - eeek                          ****
c**********************************************************************

            do 140 m=0,qTrunc
              do 130 n=-nTrunc,nTrunc

                inm = nTrunc+n+1+(2*nTrunc+1)*m

                do 120 q=0,qTrunc
                  do 110 p=-nTrunc,nTrunc

                    ipq = nTrunc+p+1+(2*nTrunc+1)*q

                    umn = 2.0d0*(1.0d0+(-1.0d0)**(n+p))

                    if (dabs(umn).gt.1.0d-10) then

                    if (m.eq.0.and.q.eq.0) then

                      Mpqnm = bJmuqa(n,0)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (mud*bHka(n)/(kd*BHdka(n)))
     +                     *(bJdmuqa(n,0)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                 +0.5d0*dcos(t)*(bJmuqa(n-1,0)*dcos(dble(n-p-1)*t)
     +                   -bJmuqa(n+1,0)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else if (m.eq.0.and.q.gt.0) then

                      Mpqnm = bJmuqa(n,q)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (zi*muqd(q)*bHka(n)/(kd*BHdka(n)))*(
     +                      bJdmuqa(n,q)*FmqP(m,q)*dcos(dble(n-p)*t)
     +              +0.5d0*dcos(t)*(bJmuqa(n-1,q)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,q)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else if (m.gt.0.and.q.eq.0) then

                      Mpqnm = bJmuqa(n,0)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (mud*bKkma(n,m)/(kmd(m)*BKdkma(n,m)))*(
     +                        bJdmuqa(n,0)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                 +0.5d0*dcos(t)*(bJmuqa(n-1,0)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,0)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else 

                      Mpqnm = bJmuqa(n,q)*
     +                        (FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (zi*muqd(q)*bKkma(n,m)/
     +                        (kmd(m)*BKdkma(n,m)))*(
     +                    bJdmuqa(n,q)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                +0.5d0*dcos(t)*(bJmuqa(n-1,q)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,q)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    end if

                    Mpqnm = doh*umn*Mpqnm
                    Kpqnm = doh*umn*Kpqnm

                    else

                      Mpqnm = dcmplx(0.0d0,0.0d0)
                      Kpqnm = dcmplx(0.0d0,0.0d0)

                    end if

c**********************************************************************
c**** Integrate                                                    ****
c**********************************************************************

                    Mat(inm,ipq) = Mat(inm,ipq)+gw(i)*h*(Mpqnm-Kpqnm)

 110              continue
 120            continue

 130          continue
 140        continue

 150      continue
 160    continue

c**********************************************************************
c**** Set up the RHS                                               ****
c**********************************************************************

        do 180 m=0,qTrunc
          do 170 n=-nTrunc,nTrunc

            inm = nTrunc+n+1+(2*nTrunc+1)*m

            if (m.eq.0) then

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

            else
 
              Rhs(inm) = dcmplx(0.0d0,0.0d0)

            end if 

 170      continue
 180    continue

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

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

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

        do 200 q=0,qTrunc
          do 190 p=-nTrunc,nTrunc

            ipq = nTrunc+p+1+(2*nTrunc+1)*q

            bpq(p,q) = Rhs(ipq)

c           print *,p,q,bpq(p,q)

 190      continue
 200    continue

c**********************************************************************
c**** Calculate the outer coefficients                             ****
c**** There are two ways to do it                                  ****
c**********************************************************************

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

c two ways of getting the amn coefficients
            if (m.eq.0) then 
              anm(n,m) =-zexp(-zi*dble(n)*beta)*bJka(n)/bHka(n)
c             anm(n,m) =-zexp(-zi*dble(n)*beta)*bJdka(n)/bHdka(n)
            else
              anm(n,m) = dcmplx(0.0d0,0.0d0)
            end if

 230      continue
 240    continue

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

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

        do 161 j=0,iTrunc-1

          do 151 i=1,10

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

c**********************************************************************
c**** Determine the real root of the \mu-dispersion relation       ****
c**** and calculate the Bessel functions J_n(\mu_0 a)              ****
c**********************************************************************

            call MUDISPER(mud,t,hod,bigKd)

            zarg = dcmplx(mud*aod,0.0d0)
            ifail = 0
            call S17DEF(0.0d0,zarg,nTrunc+3,'U',zB1Vec,nz,ifail)

            bJmuqa(0,0) = dble(zB1Vec(1))
            bJdmuqa(0,0) =-dble(zB1Vec(2))

            do 61 n=1,nTrunc+1

              umn = (-1.0d0)**n

              bJmuqa(n,0) = dble(zB1Vec(n+1))
              bJmuqa(-n,0) = umn*bJmuqa(n,0)

              bJdmuqa(n,0) = 0.5d0*(dble(zB1Vec(n))-dble(zB1Vec(n+2)))
              bJdmuqa(-n,0) = umn*bJdmuqa(n,0)

  61        continue

c**********************************************************************
c**** Determine the imag roots of the \mu-dispersion relation      ****
c**** and calculate the Bessel functions I_n(\mu_q a)              ****
c**********************************************************************

            call MUEIGVAL(muqd,qTrunc,t,hod,bigKd)

            do 58 q=1,qTrunc
 
              zarg = dcmplx(muqd(q)*aod,0.0d0)
              ifail = 0
              call S18DEF(0.0d0,zarg,nTrunc+3,'S',zB1Vec,nz,ifail)

              bJmuqa(0,q) = dble(zB1Vec(1))
              bJdmuqa(0,q) =-zi*dble(zB1Vec(2))

              do 59 n=1,nTrunc+1

                umn = (-1.0d0)**n

                bJmuqa(n,q) = (zi**n)*dble(zB1Vec(n+1))
                bJmuqa(-n,q) = umn*bJmuqa(n,q) 

                bJdmuqa(n,q) = (zi**(n-1))* 
     +                         0.5d0*(dble(zB1Vec(n))+dble(zB1Vec(n+2)))
                bJdmuqa(-n,q) = umn*bJdmuqa(n,q) 

  59          continue

  58        continue

c Alternative computation of J_n(\mu_q(t)a)
c           do 66 q=1,qTrunc
 
c             zarg = dcmplx(0.0d0,muqd(q)*aod)
c             ifail = 0
c             call S17DEF(0.0d0,zarg,nTrunc+3,'S',zB1Vec,nz,ifail)

c             bJmuqa(0,q) = zB1Vec(1)
c             bJdmuqa(0,q) =-zB1Vec(2)

c             do 71 n=1,nTrunc+1

c               umn = (-1.0d0)**n

c               bJmuqa(n,q) = zB1Vec(n+1)
c               bJmuqa(-n,q) = umn*bJmuqa(n,q)

c               bJdmuqa(n,q) = 0.5d0*(zB1Vec(n)-zB1Vec(n+2))
c               bJdmuqa(-n,q) = umn*bJdmuqa(n,q)

c 71          continue

c 66        continue

c**********************************************************************
c**** Now compute F_{m0}^+, F_{m0}^-                               ****
c**********************************************************************

            FmqP(0,0) = (1.0d0/(N0sqrt*(kd**2-mud**2)))*
     +                (mud*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (mud*dcosh(mud)-bigKd*dsinh(mud)))

            FmqM(0,0) =-(1.0d0/(N0sqrt*(kd**2-(mud*dcos(t))**2)))*
     +                (mud*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (mud*dcosh(mud)-bigKd*dsinh(mud)))

            do 81 m=1,qTrunc
     
              FmqP(m,0) =-(1.0d0/(Nmsqrt(m)*(kmd(m)**2+mud**2)))*
     +                  (mud*dcos(kmd(m)*(hod-1.0d0))*
     +                   (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                   +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                   (mud*dcosh(mud)-bigKd*dsinh(mud)))

              FmqM(m,0) = 
     +              (1.0d0/(Nmsqrt(m)*(kmd(m)**2+(mud*dcos(t))**2)))
     +                 *(mud*dcos(kmd(m)*(hod-1.0d0))*
     +                   (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                   +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                   (mud*dcosh(mud)-bigKd*dsinh(mud)))

  81        continue

c**********************************************************************
c**** Now compute F_{mq}^+, F_{mq}^-                               ****
c**********************************************************************

            do 101 q=1,qTrunc

              FmqP(0,q) = (zi/(N0sqrt*(kd**2+muqd(q)**2)))*
     +                (muqd(q)*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

              FmqM(0,q) =-(zi/(N0sqrt*(kd**2+(muqd(q)*dcos(t))**2)))*
     +                (muqd(q)*dcosh(kd*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 -kd*dsinh(kd*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

              do 91 m=1,qTrunc

                FmqP(m,q) =-(zi/(Nmsqrt(m)*(kmd(m)**2-muqd(q)**2)))*
     +                 (muqd(q)*dcos(kmd(m)*(hod-1.0d0))*
     +                 (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                 +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                 (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

                FmqM(m,q) = 
     +                (zi/(Nmsqrt(m)*(kmd(m)**2-(muqd(q)*dcos(t))**2)))*
     +                  (muqd(q)*dcos(kmd(m)*(hod-1.0d0))*
     +                  (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                  +kmd(m)*dsin(kmd(m)*(hod-1.0d0))*
     +                  (muqd(q)*dcos(muqd(q))-bigKd*dsin(muqd(q))))

  91          continue

 101        continue

c**********************************************************************
c**** Matrix elements put together - eeek                          ****
c**********************************************************************

            do 141 m=0,qTrunc
              do 131 n=-nTrunc,nTrunc

                inm = nTrunc+n+1+(2*nTrunc+1)*m

                do 121 q=0,qTrunc
                  do 111 p=-nTrunc,nTrunc

                    ipq = nTrunc+p+1+(2*nTrunc+1)*q

                    umn = 2.0d0*(1.0d0+(-1.0d0)**(n+p))

                    if (dabs(umn).gt.1.0d-10) then

                    if (m.eq.0.and.q.eq.0) then

                      Mpqnm = bJmuqa(n,0)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (mud*bHka(n)/(kd*BHdka(n)))
     +                     *(bJdmuqa(n,0)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                 +0.5d0*dcos(t)*(bJmuqa(n-1,0)*dcos(dble(n-p-1)*t)
     +                   -bJmuqa(n+1,0)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else if (m.eq.0.and.q.gt.0) then

                      Mpqnm = bJmuqa(n,q)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (zi*muqd(q)*bHka(n)/(kd*BHdka(n)))*(
     +                      bJdmuqa(n,q)*FmqP(m,q)*dcos(dble(n-p)*t)
     +              +0.5d0*dcos(t)*(bJmuqa(n-1,q)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,q)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else if (m.gt.0.and.q.eq.0) then

                      Mpqnm = bJmuqa(n,0)
     +                        *(FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (mud*bKkma(n,m)/(kmd(m)*BKdkma(n,m)))*(
     +                        bJdmuqa(n,0)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                 +0.5d0*dcos(t)*(bJmuqa(n-1,0)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,0)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    else 

                      Mpqnm = bJmuqa(n,q)*
     +                        (FmqP(m,q)+FmqM(m,q))*dcos(dble(n-p)*t)

                      Kpqnm = (zi*muqd(q)*bKkma(n,m)/
     +                        (kmd(m)*BKdkma(n,m)))*(
     +                    bJdmuqa(n,q)*FmqP(m,q)*dcos(dble(n-p)*t)
     +                +0.5d0*dcos(t)*(bJmuqa(n-1,q)*dcos(dble(n-p-1)*t)
     +                    -bJmuqa(n+1,q)*dcos(dble(n-p+1)*t))*FmqM(m,q))

                    end if

                    Mpqnm = doh*umn*Mpqnm
                    Kpqnm = doh*umn*Kpqnm

                    else

                      Mpqnm = dcmplx(0.0d0,0.0d0)
                      Kpqnm = dcmplx(0.0d0,0.0d0)

                    end if

c**********************************************************************
c**** Integrate                                                    ****
c**********************************************************************

c two ways of getting the amn coefficients
                    if (m.eq.0) then
                      anm(n,m) = anm(n,m)+bpq(p,q)*gw(i)*h*Mpqnm/bHka(n)
c                     anm(n,m) = anm(n,m)+bpq(p,q)*gw(i)*h*Kpqnm/bHka(n)
                    else
                      anm(n,m) = anm(n,m)+bpq(p,q)*gw(i)*h*Mpqnm/
     +                                    (bKkma(n,m)*dexp(-kmd(m)*aod))
c                     anm(n,m) = anm(n,m)+bpq(p,q)*gw(i)*h*Kpqnm/
c    +                                    (bKkma(n,m)*dexp(-kmd(m)*aod))
                    end if

 111              continue
 121            continue

 131          continue
 141        continue

 151      continue
 161    continue

        do 241 m=0,qTrunc
          do 231 n=-nTrunc,nTrunc

c print the scattering coefficients
c           print *,n,m,anm(n,m)

 231      continue
 241    continue

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

        if (nkhPts.eq.1) then 

          do 350 xLoop=0,nXPts

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

            do 340 yLoop=0,nYPts

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

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

              kr = kd*rod

              if (roa.ge.1.0d0) then 

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

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

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

                bHka(0) = dcmplx(dble(zB1Vec(1)),dble(zB2Vec(1)))

                do 250 n=1,nTrunc

                  umn = (-1.0d0)**n

                  bHka(n) = dcmplx(dble(zB1Vec(n+1)),dble(zB2Vec(n+1)))
                  bHka(-n) = umn*bHka(n)

 250            continue

                do 260 n=-nTrunc,nTrunc

                  phi = phi+(zi**n)*anm(n,0)*bHka(n)*
     +                      zexp(zi*dble(n)*theta)

 260            continue

                do 290 m=1,qTrunc

                  psi0 = N0sqrt*dcos(kmh(m))/(NmSqrt(m)*dcosh(kh))

                  zarg = dcmplx(kmd(m)*rod,0.0d0)
                  ifail = 0
                  call S18DCF(0.0d0,zarg,nTrunc+1,'U',zB2Vec,nz,ifail)

                  bKkma(0,m) = dble(zB2Vec(1))

                  do 270 n=1,nTrunc

                    bKkma(n,m) = dble(zB2Vec(n+1))
                    bKkma(-n,m) = bKkma(n,m)

 270              continue

                  do 280 n=-nTrunc,nTrunc

                    phi = phi+(zi**n)*anm(n,m)*bKkma(n,m)*
     +                        zexp(zi*dble(n)*theta)*psi0

 280              continue

 290            continue

              else 

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

                phi = dcmplx(0.0d0,0.0d0)

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

                do 330 j=0,iTrunc-1

                  do 320 i=1,10

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

c**********************************************************************
c**** Determine the roots of the \mu-dispersion relation           ****
c**********************************************************************

                    call MUDISPER(mud,t,hod,bigKd)
                    call MUEIGVAL(muqd,qTrunc,t,hod,bigKd)

                    do 310 p=-nTrunc,nTrunc

                      umn = (-1.0d0)**p

                      phi = phi+gw(i)*h*bpq(p,0)*mud*(
     +                      zexp(zi*(mud*rod*dcos(theta-t)+dble(p)*t))+
     +                      zexp(zi*(mud*rod*dcos(theta+t)-dble(p)*t))+
     +                  umn*zexp(-zi*(mud*rod*dcos(theta+t)+dble(p)*t))+
     +                  umn*zexp(-zi*(mud*rod*dcos(theta-t)-dble(p)*t)))

                      do 300 q=1,qTrunc

                        phi = phi+gw(i)*h*bpq(p,q)*zi*muqd(q)
     +                               *dexp(-muqd(q)*aod)*(
     +                   zexp(-muqd(q)*rod*dcos(theta-t)+zi*dble(p)*t)+
     +                   zexp(-muqd(q)*rod*dcos(theta+t)-zi*dble(p)*t)+
     +               umn*zexp(muqd(q)*rod*dcos(theta+t)-zi*dble(p)*t)+
     +               umn*zexp(muqd(q)*rod*dcos(theta-t)+zi*dble(p)*t))

 300                  continue
 310                continue

 320              continue
 330            continue

                psi0 = dcosh(kh)/N0Sqrt
                phi = phi/psi0


              end if

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

 340        continue

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

 350      continue

        else

c uncomment to run over angles
c         goto 997

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

          do 360 n=-nTrunc,nTrunc

            sigma1 = sigma1+cdabs(anm(n,0))**2
            sigma2 = sigma2-dble(anm(n,0)*zexp(zi*dble(n)*beta))

 360      continue

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

          h = 0.5d0*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)

c**********************************************************************
c**** Determine the roots of the \mu-dispersion relation           ****
c**********************************************************************

              call MUDISPER(mud,t,hod,bigKd)
              call MUEIGVAL(muqd,qTrunc,t,hod,bigKd)

c**********************************************************************
c**** We need J_1(\mu_q(t)a) for q=0 and q >=1                     ****
c**********************************************************************

              zarg = dcmplx(mud*aod,0.0d0)
              ifail = 0
              call S17DEF(0.0d0,zarg,2,'U',zB1Vec,nz,ifail)

              bJmuqa(1,0) = dble(zB1Vec(2))

              do 370 q=1,qTrunc
 
                zarg = dcmplx(0.0d0,muqd(q)*aod)
                ifail = 0
                call S17DEF(0.0d0,zarg,2,'S',zB1Vec,nz,ifail)

                bJmuqa(1,q) = zB1Vec(2)

 370          continue

              do 410 p=-nTrunc,nTrunc

                umn = 2.0d0*(1.0d0-(-1.0d0)**p)

                if (dabs(umn).gt.1.0d-10) then

                  zarg = (bigKd*dcosh(mud)-mud*dsinh(mud))
     +                   /(mud*(dcos(t)**2))

                  Fy = Fy+gw(i)*h*bpq(p,0)*dsin(t)*dsin(dble(p)*t)
     +                               *zi*umn*zarg*bJmuqa(1,0)

                  do 380 q=1,qTrunc

                    zarg = (bigKd*dcos(muqd(q))+muqd(q)*dsin(muqd(q)))
     +                     /(zi*muqd(q)*(dcos(t)**2))

                    Fy = Fy+gw(i)*h*bpq(p,q)*dsin(t)*dsin(dble(p)*t)
     +                                 *zi*umn*zarg*bJmuqa(1,q)

 380              continue

                end if

 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*kd*bHdka(1)*N0Sqrt
     +           /(dsinh(kh)*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

c**********************************************************************
c**** or compute the diffraction coefficient                       ****
c**********************************************************************

c uncomment to run over angles
c997      do 450 p=-180,180

c           theta = pi*dble(p)/180.0d0

c           zarg = dcmplx(0.0d0,0.0d0)

c           do 440 n=-nTrunc,nTrunc

c             zarg = zarg+anm(n,0)*zexp(zi*dble(n)*theta)

c440        continue

c           write (*,'(4(1x,f12.6))') dble(s),dble(p),cdabs(zarg)
c           write (fUnitRes,'(4(1x,f12.6))') dble(s),dble(p),cdabs(zarg)

c450      continue

c         write (fUnitRes,'('''')') 
c uncomment to run over angles -- down to here

        end if

 999  continue

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

      close(funitRes)

      end


      subroutine MUDISPER(mud,t,hod,bigKd)

c**********************************************************************
c**** Real dispersion relation                                     ****
c**********************************************************************

      implicit none

      integer          n

      double precision mud,t,hod,bigKd
      double precision mudlow,mudupp,mudnew,fupp,flow,fnew
      double precision cost

      cost = dcos(t)

c**** Upper limit on real root

      call DISPER(bigKd,mudupp)
      mudupp = mudupp+1.0d-10

c**** Lower limit on real root depends on Kd

      if (bigKd.le.1.0d0) then

        mudlow = 1.0d-10

      else

c**** if Kd > 1 then determine asymptote of R(\mu) by iteration

        mudlow = bigKd

        do 10 n=1,20
          mudlow = bigKd*dtanh(mudlow)
  10    continue

        mudlow = mudlow+1.0d-10

      end if

c**** Now use a bisection method

      flow = cost*dtanh(mudlow*(hod-1.0d0)*cost)
     +      -(bigKd-mudlow*dtanh(mudlow))/(mudlow-bigKd*dtanh(mudlow))
      fupp = cost*dtanh(mudupp*(hod-1.0d0)*cost)
     +      -(bigKd-mudupp*dtanh(mudupp))/(mudupp-bigKd*dtanh(mudupp))

      mudnew = 0.5d0*(mudupp+mudlow)

      do 20 while (mudupp-mudlow.gt.1.0d-13)

        fnew = cost*dtanh(mudnew*(hod-1.0d0)*cost)
     +      -(bigKd-mudnew*dtanh(mudnew))/(mudnew-bigKd*dtanh(mudnew))

        if (fnew*fupp.gt.0.0d0) then

          fupp = fnew
          mudupp = mudnew

        else

          flow = fnew
          mudlow = mudnew

        end if

        mudnew = 0.5d0*(mudupp+mudlow)

  20  continue

      mud = mudnew

      return

      end 

      subroutine MUEIGVAL(muqd,qTrunc,t,hod,bigKd)

c**********************************************************************
c**** Real dispersion relation                                     ****
c**********************************************************************

      implicit none

      integer          maxQTrunc
      parameter        (maxQTrunc=128)

      integer          n,q,qTrunc,qR,qL

      double precision pi
      double precision mudR(maxQTrunc+1),mudL(maxQTrunc+1)
      double precision mudRL(0:maxQTrunc)
      double precision muqd(maxQTrunc),t,hod,bigKd
      double precision mudlow,mudupp,mudnew,fupp,flow,fnew
      double precision cost

c**** Housekeeping

      pi = 4.0d0*datan(1.0d0)

      cost = dcos(t)

c**** Determine the location of the first zero of R at zero
c**** or by iterating to a fixed point

      if (bigKd.ge.1.0d0) then

        mudR(1) = 1.0d-10

      else

        mudR(1) = bigKd

        do 10 n=1,20
          mudR(1) = datan(mudR(1)/bigKd)
  10    continue
          
      end if 

c**** Locate the first qTrunc asymptotes of R and L
c**** L is explicit

      mudL(1) = 0.5d0*pi/((hod-1.0d0)*cost)

      do 30 q=2,qTrunc+1

        mudL(q) = (dble(q)-0.5d0)*pi/((hod-1.0d0)*cost)

c**** Determine the location of the other zeros of R 
c**** by iterating to a fixed point

        mudR(q) = (dble(q-1)+0.25d0)*pi

        do 20 n=1,20
          mudR(q) = datan(mudR(q)/bigKd)+dble(q-1)*pi
  20    continue
          
  30  continue
    
c**** We need to order the R zeros and L asymptotes

      qL = 1
      qR = 1

      do 40 q=0,qTrunc

        if (mudR(qR).lt.mudL(qL)) then

          mudRL(q) = mudR(qR)
          qR = qR+1
 
        else

          mudRL(q) = mudL(qL)
          qL = qL+1

        end if
 
  40  continue

c**** NOW we can find the roots of F 

      do 60 q=1,qTrunc

c**** The roots lie in the between zeros of R and asymptotes of L

        mudlow = mudRL(q-1)+1.0d-10
        mudupp = mudRL(q)-1.0d-10

c**** Now use a secant/bisection method

        flow = cost*dsin(mudlow*(hod-1.0d0)*cost)*
     +              (bigKd*dsin(mudlow)-mudlow*dcos(mudlow))-
     +             ((bigKd*dcos(mudlow)+mudlow*dsin(mudlow))*
     +              dcos(mudlow*(hod-1.0d0)*cost))
        fupp = cost*dsin(mudupp*(hod-1.0d0)*cost)*
     +              (bigKd*dsin(mudupp)-mudupp*dcos(mudupp))-
     +             ((bigKd*dcos(mudupp)+mudupp*dsin(mudupp))*
     +              dcos(mudupp*(hod-1.0d0)*cost))

        mudnew = 0.5d0*(mudupp+mudlow)

        do 50 while (mudupp-mudlow.gt.1.0d-13)

          fnew = cost*dsin(mudnew*(hod-1.0d0)*cost)*
     +              (bigKd*dsin(mudnew)-mudnew*dcos(mudnew))-
     +             ((bigKd*dcos(mudnew)+mudnew*dsin(mudnew))*
     +              dcos(mudnew*(hod-1.0d0)*cost))

          if (fnew*fupp.gt.0.0d0) then

            fupp = fnew
            mudupp = mudnew

          else

            flow = fnew
            mudlow = mudnew

          end if

          mudnew = 0.5d0*(mudupp+mudlow)

  50    continue

        muqd(q) = mudnew

  60  continue

c**** All done now

      return

      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
