      program resslit

c**********************************************************************
c**** Program to compute the complex eigenvalues for a             ****
c**** opening of a cavity with a divider                           ****
c**********************************************************************

      implicit none
      include 'trunc.h'

      integer          funitRes
      integer          nkdPts,kdLoop
      integer          pTrunc,nTrunc,np,nh,nz,ifail,info
      integer          p,q,n,i,j,k
      integer          niter
      integer          pivotA(2*maxPTrunc+2)

      double complex   zi,zarg,zt1,zt2,ztanh
      double precision tol
      double precision pi,c1o6,c7o6,c1o3,piby2
      double precision kdStart,kdStep,beta0d
      double complex   kd,ka,kdOld,kdEps
      double precision aod,bod,theta0
      double precision gp(10),gw(10)
      double precision dnpi,dnpi2
      double complex   alphand(maxNTrunc),kappand(maxNTrunc)

      double complex   zBesVec(2*maxPTrunc+2)
      double precision Fpn1(0:maxPTrunc,maxNTrunc)
      double precision Fpn2(0:maxPTrunc,maxNTrunc)
      double precision Gpn1(0:maxPTrunc,maxNTrunc)
      double precision Gpn2(0:maxPTrunc,maxNTrunc)

      double precision tlim,h,t
      double complex   tint
      double complex   KpqA(2*maxPTrunc+2,2*maxPTrunc+2)
      double complex   RhsA(2*maxPTrunc+2)
      double complex   f,f1
      double precision aa0A
      double complex   a0A,a1A,qsuma

      double precision S14AAF
      external         S14AAF

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

      pi = 4.0d0*datan(1.0d0)
      piby2 = 0.5d0*pi
      c1o6 = 1.0d0/6.0d0
      c7o6 = 7.0d0/6.0d0
      c1o3 = 1.0d0/3.0d0
      funitRes = 11
      zi = dcmplx(0.0d0,1.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**** Title and data input and the convert theta to radians        ****
c**********************************************************************

      call DATA_INPUT(funitRes,pTrunc,nTrunc,np,aod,bod,theta0,
     +                kdStart,kdStep,nkdPts)

c**********************************************************************
c**** Loop through complex values of kd via Newton method with     ****
c**** initial guess                                                ****
c**********************************************************************

c++++ MANUAL INSERT

      do 998 k=0,1

      aod = 0.34d0-dble(k)/dble(100)

c++++ END MANUAL INSERT

      tol = 2.0d-13
      niter = 0

      kd = dcmplx(kdStart,-0.1d0)
      kdEps = tol*dcmplx(1.0d0,1.0d0)

      kdOld = dcmplx(kdStart,10.0d0)

      do 999 while (cdabs(kd-kdOld).gt.1.0d-13.and.niter.lt.nkdPts)

        kdOld = kd

        ka = kd*aod

c**********************************************************************
c**** Find \alpha_n, \kappa_n and F_{pn}^{(1)}                     ****
c**********************************************************************

        do 50 n=1,nTrunc   

          dnpi = dble(n)*pi
          dnpi2 = (dble(n)-0.5d0)*pi

          alphand(n) = zsqrt((dnpi/aod)**2-kd**2)
          kappand(n) = zsqrt((dnpi2/aod)**2-kd**2)

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(c1o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 10 p=0,pTrunc
            Fpn1(p,n) = dble(zBesVec(2*p+1))/(dnpi**c1o6)
 10       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(c7o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 20 p=0,pTrunc
            Fpn2(p,n) = dble(zBesVec(2*p+1))/(dnpi2**c1o6)
 20       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 30 p=0,pTrunc
            Gpn1(p,n) =-((-1.0d0)**n)*dble(zBesVec(2*p+1))
 30       continue

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 40 p=0,pTrunc
            Gpn2(p,n) = ((-1.0d0)**n)*dble(zBesVec(2*p+1))
 40       continue

 50     continue
     
c**********************************************************************
c**** Initialise arrays KS(p,q) and KA(p,q) and then add any       ****
c**** infinite sums                                                ****
c**********************************************************************

        do 80 p=0,pTrunc
          do 70 q=0,pTrunc

            KpqA(2*p+1,2*q+1) = 0.0d0
            KpqA(2*p+1,2*q+2) = 0.0d0
            KpqA(2*p+2,2*q+1) = 0.0d0
            KpqA(2*p+2,2*q+2) = 0.0d0

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

               zt1 = zsin(kd*(1.0d0-bod))/zcos(kd*(1.0d0-bod))
               KpqA(2*p+1,2*q+1) = 1.0d0/(ka*zt1)

            end if

            do 60 n=1,nTrunc

              zt1 = ztanh(kappand(n)*bod)
              zt2 = ztanh(alphand(n)*(1.0d0-bod))

              KpqA(2*p+1,2*q+1) = KpqA(2*p+1,2*q+1)
     +                -2.0d0*Gpn1(p,n)*Gpn1(q,n)/(kappand(n)*aod*zt1)
     +                -2.0d0*Gpn2(p,n)*Gpn2(q,n)/(alphand(n)*aod*zt2)

              KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)
     +                +2.0d0*Fpn2(p,n)*Fpn2(q,n)/(kappand(n)*aod*zt1)

              if (cdabs(kappand(n)*bod).lt.100d0) then 
                zt1 = 1.0d0/(-zi*zsin(zi*kappand(n)*bod))
              else
                zt1 = 0.0d0
              end if 

              KpqA(2*p+1,2*q+2) = KpqA(2*p+1,2*q+2)
     +                -2.0d0*Gpn1(p,n)*Fpn2(q,n)*zt1/(kappand(n)*aod)

              KpqA(2*p+2,2*q+1) = KpqA(2*p+2,2*q+1)
     +                -2.0d0*Fpn2(p,n)*Gpn1(q,n)*zt1/(kappand(n)*aod)

 60         continue

 70       continue
 80     continue

c**********************************************************************
c**** The next bit calculates the integrals in 3 parts             ****
c**** 1. we start with the finite integral over (0, kd) maps to    ****
c****    0 < t < \pi/2                                             ****
c**********************************************************************

        tlim = pi/2.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 150 j=0,nh-1

          do 140 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dsin(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dcos(t)/
     +             (zarg*(dble(ka)*dsin(t))**c1o3)

            do 130 p=0,pTrunc
              do 120 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 120          continue
 130        continue

 140      continue
 150    continue

c**********************************************************************
c**** 2. next, do (kd,2kd)                                         ****
c**********************************************************************

        tlim = dlog(2.0d0+dsqrt(3.0d0))
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 190 j=0,nh-1

          do 180 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dcosh(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dsinh(t)/
     +             (zarg*(dble(ka)*dcosh(t))**c1o3)

            do 170 p=0,pTrunc
              do 160 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 160          continue
 170        continue

 180      continue
 190    continue

c**********************************************************************
c**** 3. finally, do (2kd,infty)                                   ****
c**********************************************************************

        tlim = 1000.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 230 j=0,nh-1

          do 220 i=1,10

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

            ifail = 0
            zarg = dcmplx(t,0.0d0)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(t**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if

            do 210 p=0,pTrunc
              do 200 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)
                tint = 2.0d0*tint/(pi*zarg*(t**c1o3))

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 200          continue
 210        continue

 220      continue
 230    continue

c**********************************************************************
c**** Invert matrices                                              ****
c**********************************************************************

        call CLUDCMP(KpqA,2*pTrunc+2,2*maxPTrunc+2,pivotA,info)

c**********************************************************************
c**** Calculate the monopolar/dipolar coefficients                 ****
c**********************************************************************

        f = dcmplx(dble(info),0.0d0)
        do 240 p=1,2*pTrunc+2

          f = f*KpqA(p,p)

 240    continue

        kd = kdOld+kdEps

        ka = kd*aod

c**********************************************************************
c**** Find \alpha_n, \kappa_n and F_{pn}^{(1)}                     ****
c**********************************************************************

        do 51 n=1,nTrunc   

          dnpi = dble(n)*pi
          dnpi2 = (dble(n)-0.5d0)*pi

          alphand(n) = zsqrt((dnpi/aod)**2-kd**2)
          kappand(n) = zsqrt((dnpi2/aod)**2-kd**2)

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(c1o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 11 p=0,pTrunc
            Fpn1(p,n) = dble(zBesVec(2*p+1))/(dnpi**c1o6)
 11       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(c7o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 21 p=0,pTrunc
            Fpn2(p,n) = dble(zBesVec(2*p+1))/(dnpi2**c1o6)
 21       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 31 p=0,pTrunc
            Gpn1(p,n) =-((-1.0d0)**n)*dble(zBesVec(2*p+1))
 31       continue

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 41 p=0,pTrunc
            Gpn2(p,n) = ((-1.0d0)**n)*dble(zBesVec(2*p+1))
 41       continue

 51     continue
     
c**********************************************************************
c**** Initialise arrays KS(p,q) and KA(p,q) and then add any       ****
c**** infinite sums                                                ****
c**********************************************************************

        do 81 p=0,pTrunc
          do 71 q=0,pTrunc

            KpqA(2*p+1,2*q+1) = 0.0d0
            KpqA(2*p+1,2*q+2) = 0.0d0
            KpqA(2*p+2,2*q+1) = 0.0d0
            KpqA(2*p+2,2*q+2) = 0.0d0

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

               zt1 = zsin(kd*(1.0d0-bod))/zcos(kd*(1.0d0-bod))
               KpqA(2*p+1,2*q+1) = 1.0d0/(ka*zt1)

            end if

            do 61 n=1,nTrunc

              zt1 = ztanh(kappand(n)*bod)
              zt2 = ztanh(alphand(n)*(1.0d0-bod))

              KpqA(2*p+1,2*q+1) = KpqA(2*p+1,2*q+1)
     +                -2.0d0*Gpn1(p,n)*Gpn1(q,n)/(kappand(n)*aod*zt1)
     +                -2.0d0*Gpn2(p,n)*Gpn2(q,n)/(alphand(n)*aod*zt2)

              KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)
     +                +2.0d0*Fpn2(p,n)*Fpn2(q,n)/(kappand(n)*aod*zt1)

              if (cdabs(kappand(n)*bod).lt.100d0) then 
                zt1 = 1.0d0/(-zi*zsin(zi*kappand(n)*bod))
              else
                zt1 = 0.0d0
              end if 

              KpqA(2*p+1,2*q+2) = KpqA(2*p+1,2*q+2)
     +                -2.0d0*Gpn1(p,n)*Fpn2(q,n)*zt1/(kappand(n)*aod)

              KpqA(2*p+2,2*q+1) = KpqA(2*p+2,2*q+1)
     +                -2.0d0*Fpn2(p,n)*Gpn1(q,n)*zt1/(kappand(n)*aod)

 61         continue

 71       continue
 81     continue

c**********************************************************************
c**** The next bit calculates the integrals in 3 parts             ****
c**** 1. we start with the finite integral over (0, kd) maps to    ****
c****    0 < t < \pi/2                                             ****
c**********************************************************************

        tlim = pi/2.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 151 j=0,nh-1

          do 141 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dsin(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dcos(t)/
     +             (zarg*(dble(ka)*dsin(t))**c1o3)

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

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 121          continue
 131        continue

 141      continue
 151    continue

c**********************************************************************
c**** 2. next, do (kd,2kd)                                         ****
c**********************************************************************

        tlim = dlog(2.0d0+dsqrt(3.0d0))
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 191 j=0,nh-1

          do 181 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dcosh(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dsinh(t)/
     +             (zarg*(dble(ka)*dcosh(t))**c1o3)

            do 171 p=0,pTrunc
              do 161 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 161          continue
 171        continue

 181      continue
 191    continue

c**********************************************************************
c**** 3. finally, do (2kd,infty)                                   ****
c**********************************************************************

        tlim = 1000.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 231 j=0,nh-1

          do 221 i=1,10

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

            ifail = 0
            zarg = dcmplx(t,0.0d0)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(t**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if

            do 211 p=0,pTrunc
              do 201 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)
                tint = 2.0d0*tint/(pi*zarg*(t**c1o3))

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 201          continue
 211        continue

 221      continue
 231    continue

c**********************************************************************
c**** Invert matrices                                              ****
c**********************************************************************

        call CLUDCMP(KpqA,2*pTrunc+2,2*maxPTrunc+2,pivotA,info)

c**********************************************************************
c**** Calculate the monopolar/dipolar coefficients                 ****
c**********************************************************************

        f1 = dcmplx(dble(info),0.0d0)
        do 241 p=1,2*pTrunc+2

          f1 = f1*KpqA(p,p)

 241    continue

c++++ Newton step

        kd = kdOld-f*kdEps/(f1-f)

        niter = niter+1

 999  continue

c++++ NOW RUN SCATTERING PROBLEM AT REAL VALUE OF kd

      kd = dble(kd)

        ka = kd*aod
      beta0d = dble(kd)*dsin(theta0)

c**********************************************************************
c**** Find \alpha_n, \kappa_n and F_{pn}^{(1)}                     ****
c**********************************************************************

        do 52 n=1,nTrunc   

          dnpi = dble(n)*pi
          dnpi2 = (dble(n)-0.5d0)*pi

          alphand(n) = zsqrt((dnpi/aod)**2-kd**2)
          kappand(n) = zsqrt((dnpi2/aod)**2-kd**2)

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(c1o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 12 p=0,pTrunc
            Fpn1(p,n) = dble(zBesVec(2*p+1))/(dnpi**c1o6)
 12       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(c7o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 22 p=0,pTrunc
            Fpn2(p,n) = dble(zBesVec(2*p+1))/(dnpi2**c1o6)
 22       continue

          ifail = 0
          zarg = dcmplx(dnpi2,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 32 p=0,pTrunc
            Gpn1(p,n) =-((-1.0d0)**n)*dble(zBesVec(2*p+1))
 32       continue

          ifail = 0
          zarg = dcmplx(dnpi,0.0d0)
          call S17DEF(0.0d0,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 42 p=0,pTrunc
            Gpn2(p,n) = ((-1.0d0)**n)*dble(zBesVec(2*p+1))
 42       continue

 52     continue
     
c**********************************************************************
c**** Initialise arrays KS(p,q) and KA(p,q) and then add any       ****
c**** infinite sums                                                ****
c**********************************************************************

        do 82 p=0,pTrunc
          do 72 q=0,pTrunc

            KpqA(2*p+1,2*q+1) = 0.0d0
            KpqA(2*p+1,2*q+2) = 0.0d0
            KpqA(2*p+2,2*q+1) = 0.0d0
            KpqA(2*p+2,2*q+2) = 0.0d0

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

               zt1 = zsin(kd*(1.0d0-bod))/zcos(kd*(1.0d0-bod))
               KpqA(2*p+1,2*q+1) = 1.0d0/(ka*zt1)

            end if

            do 62 n=1,nTrunc

              zt1 = ztanh(kappand(n)*bod)
              zt2 = ztanh(alphand(n)*(1.0d0-bod))

              KpqA(2*p+1,2*q+1) = KpqA(2*p+1,2*q+1)
     +                -2.0d0*Gpn1(p,n)*Gpn1(q,n)/(kappand(n)*aod*zt1)
     +                -2.0d0*Gpn2(p,n)*Gpn2(q,n)/(alphand(n)*aod*zt2)

              KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)
     +                +2.0d0*Fpn2(p,n)*Fpn2(q,n)/(kappand(n)*aod*zt1)

              if (cdabs(kappand(n)*bod).lt.100d0) then 
                zt1 = 1.0d0/(-zi*zsin(zi*kappand(n)*bod))
              else
                zt1 = 0.0d0
              end if 

              KpqA(2*p+1,2*q+2) = KpqA(2*p+1,2*q+2)
     +                -2.0d0*Gpn1(p,n)*Fpn2(q,n)*zt1/(kappand(n)*aod)

              KpqA(2*p+2,2*q+1) = KpqA(2*p+2,2*q+1)
     +                -2.0d0*Fpn2(p,n)*Gpn1(q,n)*zt1/(kappand(n)*aod)

 62         continue

 72       continue
 82     continue

c**********************************************************************
c**** Find RHS F_{p0}^{(2)} depends on \theta_0 = 0                ****
c**********************************************************************

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

          ifail = 0
          zarg = dcmplx(beta0d*aod,0.0d0)
          call S17DEF(c7o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 102 p=0,pTrunc

            RhsA(2*p+1) = dcmplx(0.0d0,0.0d0)
            RhsA(2*p+2) =-2.0d0*zi*dble(zBesVec(2*p+1))/
     +                              ((beta0d*aod)**c1o6)

 102      continue

        else

          do 112 p=0,pTrunc

            RhsA(2*p+1) = dcmplx(0.0d0,0.0d0)
            RhsA(2*p+2) = dcmplx(0.0d0,0.0d0)

 112      continue

        end if

c**********************************************************************
c**** The next bit calculates the integrals in 3 parts             ****
c**** 1. we start with the finite integral over (0, kd) maps to    ****
c****    0 < t < \pi/2                                             ****
c**********************************************************************

        tlim = pi/2.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 152 j=0,nh-1

          do 142 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dsin(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dcos(t)/
     +             (zarg*(dble(ka)*dsin(t))**c1o3)

            do 132 p=0,pTrunc
              do 122 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 122          continue
 132        continue

 142      continue
 152    continue

c**********************************************************************
c**** 2. next, do (kd,2kd)                                         ****
c**********************************************************************

        tlim = dlog(2.0d0+dsqrt(3.0d0))
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 192 j=0,nh-1

          do 182 i=1,10

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

            ifail = 0
            zarg = dble(ka)*dcosh(t)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(zarg**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if
            zarg = (2.0d0/pi)*dble(ka)*dsinh(t)/
     +             (zarg*(dble(ka)*dcosh(t))**c1o3)

            do 172 p=0,pTrunc
              do 162 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)*zarg

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 162          continue
 172        continue

 182      continue
 192    continue

c**********************************************************************
c**** 3. finally, do (2kd,infty)                                   ****
c**********************************************************************

        tlim = 1000.0d0
        h = 2.0d0*pi/dble(np)
        nh = int(tlim/(2.0d0*h))
        h = tlim/(2.0d0*dble(nh))

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

        do 232 j=0,nh-1

          do 222 i=1,10

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

            ifail = 0
            zarg = dcmplx(t,0.0d0)
            call S17DEF(c1o6,zarg,2*pTrunc+2,'U',zBesVec,nz,ifail)

            zarg = zsqrt(t**2-ka**2)
            if (dimag(zarg).lt.0.0d0) then
              zarg =-zarg
            end if

            do 212 p=0,pTrunc
              do 202 q=0,pTrunc

                tint = zBesVec(2*p+2)*zBesVec(2*q+2)
                tint = 2.0d0*tint/(pi*zarg*(t**c1o3))

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+h*gw(i)*tint

 202          continue
 212        continue

 222      continue
 232    continue

c**********************************************************************
c**** Invert matrices                                              ****
c**********************************************************************

        call CLUDCMP(KpqA,2*pTrunc+2,2*maxPTrunc+2,pivotA,info)
        call CLUBKSB(KpqA,2*pTrunc+2,2*maxPTrunc+2,pivotA,RhsA)

c**********************************************************************
c**** Calculate the monopolar/dipolar coefficients                 ****
c**********************************************************************

c Leading order
        ifail = 0
        a1A =-RhsA(2)*ka/
     +        (dsin(theta0)*(2.0d0**c7o6)*S14AAF(13.0d0/6.0d0,ifail))
c divided by sin theta0
        a0A =-RhsA(1)/(ka*dsin(dble(kd)*(1.0d0-bod)))
        aa0A =-pi*dexp(1.0d0)*dsin(theta0)*dexp(0.5d0*pi*bod/aod)/
     +         ((aod*0.5d0*pi/(1.0d0-bod))**2)

        ifail = 0
        a1A = dcmplx(0.0d0,0.0d0)

c Exact

        do 382 p=0,pTrunc

          qsuma = 0.0d0
          do 381 q=0,pTrunc

            qsuma = qsuma+((-0.25d0)**q*((0.25d0*ka)**(2*q+2*p))*
     +                    S14AAF(2.0d0*(q+p)+3.0d0,ifail))/
     +      (S14AAF(dble(q)+1.0d0,ifail)*
     +       S14AAF(dble(q)+2.0d0*dble(p)+13.0d0/6.0d0,ifail)*
     +       (S14AAF(dble(q)+dble(p)+2.0d0,ifail))**2)

 381      continue
          a1A = a1A+0.25d0*zi*ka*qsuma*RhsA(2*p+2)/2.0d0**c1o6

 382    continue

c Divided by sin(theta0)
        a1A = a1A/dsin(theta0)

c**********************************************************************
c**** Print to screen/file                                         ****
c**********************************************************************

        write (*,'(6(1x,g16.10))') 
     +           aod,dble(kd),cdabs(a1a),cdabs(a0a),dabs(aa0A)

        write (fUnitRes,'(5(1x,g16.10))') 
     +           aod,dble(kd),cdabs(a1a),cdabs(a0a),dabs(aa0A)

 998  continue

c++++ END MANUAL INSERT

c**********************************************************************
c**** End program                                                  ****
c**********************************************************************
    
      close(funitRes)

      end

      double complex function ztanh(z)

      double complex z
      double precision x,y,ex

      x = dble(z)
      y = dimag(z)

      ex = dexp(-2.0d0*dabs(x))

      ztanh = dcmplx(
     +     dtanh(2.0d0*x)/(1.0d0+2.0d0*ex*dcos(2.0d0*y)/(1.0d0+ex**2)),
     +     dsin(2.0d0*y)*ex/(0.5d0+0.5d0*ex**2+dcos(2.0d0*y)*ex))


      return
      end

