      program dslit

c**********************************************************************
c**** Program to compute the scattering of acoustic waves by       ****
c**** the 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
      integer          pivotS(maxPTrunc+1)
      integer          pivotA(2*maxPTrunc+2)

      double complex   zi,zarg
      double precision pi,c1o6,c7o6,c1o3,piby2
      double precision kd,ka,kdStart,kdStep
      double precision aod,bod,theta0
      double precision gp(10),gw(10)
      double precision beta0d
      double precision dnpi,dnpi2,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,tint
      double complex   KpqS(maxPTrunc+1,maxPTrunc+1)
      double complex   RhsS(maxPTrunc+1)
      double complex   KpqA(2*maxPTrunc+2,2*maxPTrunc+2)
      double complex   RhsA(2*maxPTrunc+2)
      double complex   TMat(maxPTrunc+1,maxPTrunc+1)
      double complex   a0S,a0A,a1A,det,qsums,qsuma
      double precision chi,elu,aa0A

      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)

      theta0 = pi*theta0/180.0d0

c**********************************************************************
c**** Loop through values of kd                                    ****
c**********************************************************************

      do 999 kdLoop=1,nkdPts

        kd = kdStart+dble(kdLoop-1)*kdStep
        ka = kd*aod

        beta0d = kd*dsin(theta0)

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) = dsqrt((dnpi/aod)**2-kd**2)
          kappand(n) = dsqrt((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

            KpqS(p+1,q+1) = 0.0d0

            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
            TMat(p+1,q+1) = 0.0d0

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

               KpqS(p+1,q+1) =-(1.0d0/(ka*dtan(kd)))
     +                   *(6.0d0/((2.0d0**c1o6)*S14AAF(c1o6,ifail)))**2

               KpqA(2*p+1,2*q+1) = 1.0d0/(ka*dtan(kd*(1.0d0-bod)))
               TMat(p+1,q+1) = 1.0d0/(ka*dtan(kd*(1.0d0-bod)))

            end if

            do 60 n=1,nTrunc

              KpqS(p+1,q+1) = KpqS(p+1,q+1)+2.0d0*Fpn1(p,n)*Fpn1(q,n)/
     +                        (alphand(n)*aod*dtanh(alphand(n)))

              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*dtanh(kappand(n)*bod))
     +                    -2.0d0*Gpn2(p,n)*Gpn2(q,n)/
     +                 (alphand(n)*aod*dtanh(alphand(n)*(1.0d0-bod)))

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

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

              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*dtanh(kappand(n)*bod))

              TMat(p+1,q+1) = TMat(p+1,q+1)
     +                    -2.0d0*Gpn1(p,n)*Gpn1(q,n)/(kappand(n)*aod)
     +                    -2.0d0*Gpn2(p,n)*Gpn2(q,n)/
     +                 (alphand(n)*aod*dtanh(alphand(n)*(1.0d0-bod)))

 60         continue

 70       continue
 80     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(c1o6,zarg,2*pTrunc+1,'U',zBesVec,nz,ifail)

          do 90 p=0,pTrunc
            RhsS(p+1) = 2.0d0*dble(zBesVec(2*p+1))/
     +                        ((beta0d*aod)**c1o6)
 90       continue

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

          do 100 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)

 100      continue

        else

          do 110 p=0,pTrunc

            RhsS(p+1) = dcmplx(0.0d0,0.0d0)

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

 110      continue

          ifail = 0
          RhsS(1) = 12.0d0/((2.0d0**c1o6)*S14AAF(c1o6,ifail))

        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 150 j=0,nh-1

          do 140 i=1,10

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

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

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

                tint = dble(zBesVec(2*p+1))*dble(zBesVec(2*q+1))
                tint = (2.0d0/pi)*tint/((ka*dsin(t))**c1o3)

                KpqS(p+1,q+1) = KpqS(p+1,q+1)+zi*h*gw(i)*tint

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

                KpqA(2*p+2,2*q+2) = KpqA(2*p+2,2*q+2)+zi*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 151 j=0,nh-1

          do 141 i=1,10

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

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

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

                tint = dble(zBesVec(2*p+1))*dble(zBesVec(2*q+1))
                tint = (2.0d0/pi)*tint/((ka*dcosh(t))**c1o3)

                KpqS(p+1,q+1) = KpqS(p+1,q+1)+h*gw(i)*tint

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

                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**** 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 152 j=0,nh-1

          do 142 i=1,10

            t = 2.0d0*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)

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

                tint = dble(zBesVec(2*p+1))*dble(zBesVec(2*q+1))
                tint = 2.0d0*tint/(pi*dsqrt(t**2-ka**2)*(t**c1o3))

                KpqS(p+1,q+1) = KpqS(p+1,q+1)+h*gw(i)*tint

                tint = dble(zBesVec(2*p+2))*dble(zBesVec(2*q+2))
                tint = 2.0d0*tint/(pi*dsqrt(t**2-ka**2)*(t**c1o3))

                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**** Special cases                                                ****
c**********************************************************************


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

c           UMat(p+1,q+1) = KpqA(2*p+2,2*q+2)

 170      continue

c         RhsU(p+1) = RhsA(2*p+2)

 180    continue

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

        call CLUDCMP(KpqS,pTrunc+1,maxPTrunc+1,pivotS,info)
        call CLUBKSB(KpqS,pTrunc+1,maxPTrunc+1,pivotS,RhsS)

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

        zarg = dcmplx(dble(info),0.0d0)

c       call CLUDCMP(UMat,pTrunc+1,maxPTrunc+1,pivotS,info)
c       call CLUBKSB(UMat,pTrunc+1,maxPTrunc+1,pivotS,RhsU)

        call CLUDCMP(TMat,pTrunc+1,maxPTrunc+1,pivotS,info)

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

        ifail = 0
        a0S = 6.0d0*RhsS(1)*zi/((2.0d0**c1o6)*S14AAF(c1o6,ifail))
        a1A =-RhsA(2)*ka/
     +        ((2.0d0**c7o6)*S14AAF(13.0d0/6.0d0,ifail))
        a0A =-RhsA(1)/(ka*dsin(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)

        print *,'A0S = ',a0S
        print *,'A1A = ',a1A

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

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

        do 182 p=0,pTrunc

          qsums = 0.0d0
          qsuma = 0.0d0
          do 181 q=0,pTrunc

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

            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)

 181      continue
          a0S = a0S+zi*qsums*RhsS(p+1)/2.0d0**c1o6
          a1A = a1A+0.25d0*zi*ka*qsuma*RhsA(2*p+2)/2.0d0**c1o6

 182    continue
        print *,'A0S = ',a0S
        print *,'A1A = ',a1A

c++++ to compare with -4i
        a1A = a1A/dsin(theta0)

c**********************************************************************
c**** Calculate the trapped mode condition for an infinite channel ****
c**** And zarg is the condition for a pure antisymmetric modes in  ****
c**** and open channel... small Im part implies resonance          ****
c**********************************************************************

        det = dble(info)
        do 190 p=0,pTrunc

          det = det*TMat(p+1,p+1)
          zarg = zarg*KpqA(p+1,p+1)*KpqA(pTrunc+p+2,pTrunc+p+2)

 190    continue

c       print *,dble(det)
c       print *,zarg

c**********************************************************************
c**** Calculate the Evans, Linton and Ursell condition             ****
c**** based on WH and wide spacing                                 ****
c**********************************************************************

        chi = 0.0d0
        do 200 n=1,nTrunc

          dnpi = dble(n)*pi
          chi = 2.0d0*dasin(ka/dnpi)-dasin(2.0d0*ka/dnpi)

 200    continue

        elu = kd*(1.0d0-bod)-pi/2.0d0+ka*dlog(4.0d0)/pi

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

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

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

 999  continue

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

      end
