      program metaswe

c**********************************************************************
c**** Program to compute the waves along a string with multiple    ****
c**** finite branches                                              ****
c****                                                              ****
c**** Written 19/02/2021                                           ****
c**********************************************************************

      implicit none

      integer          funitRes
      integer          n
      integer          nk0h0Pts,k0h0Loop
      integer          ntheta0Pts,theta0Loop
      integer          xoh0Loop,yoh0Loop

      double precision pi,arg
      double complex   zi,zarg
      double precision h0oh1p,h0oh1m,h1poh0,h1moh0,h1baroh0
      double precision delta,bigTheta,loh0
      double precision h11hat,h22hat,h12hat
      double precision k0h0,k0h0Start,k0h0Step
      double precision theta0,theta0Start,theta0Step
      double precision k1h1p,k1h1m
      double precision alpha0h0,alpha1h1p,alpha1h1m
      double precision betah0,betah1p,betah1m
      double precision theta1,theta1p,theta1m
      double precision disc,cgx
      double complex   zalpha1h0,zzalpha1h0
      double complex   zEp,zEm,zF
      double complex   zRp,zTp,zRm,zTm,zDen,zRMM,zTMM,zAMM,zBMM
      double complex   zeta
      double precision xoh0,yoh0
      double precision absRp,absTp,absRm,absTm,absRMM,absTMM

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

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

c**********************************************************************
c**** Title and data input                                         ****
c**********************************************************************

      call DATA_INPUT(funitRes,h0oh1p,h0oh1m,delta,bigTheta,loh0,
     +                      theta0Start,theta0Step,ntheta0Pts,
     +                      k0h0Start,k0h0Step,nk0h0Pts)

      theta0Start = theta0Start*pi/180.0d0
      theta0Step = theta0Step*pi/180.0d0

      delta = delta*pi/180.0d0

      h1poh0 = 1.0d0/h0oh1p
      h1moh0 = 1.0d0/h0oh1m

      h1baroh0 = bigTheta*h1poh0+(1.0d0-bigTheta)*h1moh0

c     h1baroh0 = (1.0d0+dsqrt(1.0d0-4.0d0*(dsin(delta)*dcos(delta))**2))
c    +           /(2.0d0*dsin(delta)**2)
c     h1poh0 = (1.0d0-h1baroh0*dsin(delta)**2)/(dcos(delta)**2)
c     print *,h1poh0,h1baroh0

      h11hat = h1poh0*dcos(delta)**2+h1baroh0*dsin(delta)**2
      h22hat = h1poh0*dsin(delta)**2+h1baroh0*dcos(delta)**2
      h12hat = (h1poh0-h1baroh0)*dsin(delta)*dcos(delta)

c     print *,h1poh0,h1baroh0
c     print *,h12hat**2-h11hat*h22hat,h11hat

c**********************************************************************
c**** Loop through values of k0h0 and theta0                       ****
c**********************************************************************

      do 999 k0h0Loop=1,nk0h0Pts

        k0h0 = k0h0Start+dble(k0h0Loop-1)*k0h0Step
        k1h1p = dsqrt(h1poh0)*k0h0
        k1h1m = dsqrt(h1moh0)*k0h0

        do 998 theta0Loop=1,ntheta0Pts

          theta0 = theta0Start+dble(theta0Loop-1)*theta0Step

          alpha0h0 = k0h0*dcos(theta0)
          betah0 = k0h0*dsin(theta0)

c**********************************************************************
c**** Normal bed height h1p extending into x>0                     ****
c**********************************************************************

          if (dabs(dsin(theta0)).lt.dsqrt(h0oh1p)) then

            theta1p = dasin(dsqrt(h1poh0)*dsin(theta0))

            alpha1h1p = k1h1p*dcos(theta1p)
            betah1p = k1h1p*dsin(theta1p)

            if (dabs(betah1p*h0oh1p-betah0).gt.1.0d-12) then
              print *,betah0,betah1p*h0oh1p
            end if

            zRp = (alpha0h0-alpha1h1p)/(alpha0h0+alpha1h1p)
            zTp = 2.0d0*alpha0h0/(alpha0h0+alpha1h1p)

          else

            alpha1h1p = k0h0*h1poh0*dsqrt(dsin(theta0)**2-h0oh1p)

            zRp = (alpha0h0-zi*alpha1h1p)/(alpha0h0+zi*alpha1h1p)
            zTp = 2.0d0*alpha0h0/(alpha0h0+zi*alpha1h1p)

          endif

          absRp = cdabs(zRp)
          absTp = cdabs(zTp)

c**********************************************************************
c**** Normal bed height h1m extending into x>0                     ****
c**********************************************************************

          if (dabs(dsin(theta0)).lt.dsqrt(h0oh1m)) then

            theta1m = dasin(dsqrt(h1moh0)*dsin(theta0))

            alpha1h1m = k1h1m*dcos(theta1m)
            betah1m = k1h1m*dsin(theta1m)

            if (dabs(betah1m*h0oh1m-betah0).gt.1.0d-12) then
              print *,betah0,betah1m*h0oh1m
            end if

            zRm = (alpha0h0-alpha1h1m)/(alpha0h0+alpha1h1m)
            zTm = 2.0d0*alpha0h0/(alpha0h0+alpha1h1m)

          else

            alpha1h1m = k0h0*h1moh0*dsqrt(dsin(theta0)**2-h0oh1m)

            zRm = (alpha0h0-zi*alpha1h1m)/(alpha0h0+zi*alpha1h1m)
            zTm = 2.0d0*alpha0h0/(alpha0h0+zi*alpha1h1m)

          endif

          absRm = cdabs(zRm)
          absTm = cdabs(zTm)

c**********************************************************************
c**** Metamaterial bed transmission into x>0                       ****
c**********************************************************************

          disc = (betah0*h12hat)**2-h11hat*(h22hat*betah0**2-k0h0**2)

          if (disc.gt.0.0d0) then 
          
            zalpha1h0 = (-h12hat*betah0+dsqrt(disc))/h11hat

          else

            zalpha1h0 = (-h12hat*betah0+zi*dsqrt(-disc))/h11hat

          end if
          
c group velocity component in the x-direction, should never be -ve
          cgx = dble(zalpha1h0)*h11hat+betah0*h12hat

          if (cgx.lt.-1.0d-12) then 
            print *,'Oooops'
          end if

          theta1 = 180.0d0*datan2(h12hat*dble(zalpha1h0)+h22hat*betah0,
     +                           cgx)/pi

          zRMM = (alpha0h0-zalpha1h0*h11hat-betah0*h12hat)/
     +           (alpha0h0+zalpha1h0*h11hat+betah0*h12hat)
          zTMM = 2.0d0*alpha0h0/
     +           (alpha0h0+zalpha1h0*h11hat+betah0*h12hat)

          absRMM = cdabs(zRMM)
          absTMM = cdabs(zTMM)

c**********************************************************************
c**** Scattering over a finite shelf                               ****
c**********************************************************************

          if (loh0.gt.1.0d-12) then

c**********************************************************************
c**** Metamaterial bed transmission                                ****
c**********************************************************************

          disc = (betah0*h12hat)**2-h11hat*(h22hat*betah0**2-k0h0**2)

          if (disc.gt.0.0d0) then 
          
            zalpha1h0 = (-h12hat*betah0+dsqrt(disc))/h11hat
            zzalpha1h0 = (-h12hat*betah0-dsqrt(disc))/h11hat

c         print *,dsqrt(disc)*Loh0/(pi*h11hat)
          else

            zalpha1h0 = (-h12hat*betah0+zi*dsqrt(-disc))/h11hat
            zzalpha1h0 = (-h12hat*betah0-zi*dsqrt(-disc))/h11hat

          end if
          
c         print *,zalpha1h0/k0h0
c         print *,zzalpha1h0/k0h0
c         print *,datan2(betah0,alpha0h0)*180.0d0/pi
c         print *,datan2(betah0,dble(zalpha1h0))*180.0d0/pi
c         print *,datan2(betah0,dble(zzalpha1h0))*180.0d0/pi

          zEp = ((zalpha1h0*h11hat+betah0*h12hat)/alpha0h0)-1.0d0
          zEm = ((zzalpha1h0*h11hat+betah0*h12hat)/alpha0h0)-1.0d0
          zF = zexp(zi*(zzalpha1h0-zalpha1h0)*loh0)

          zDen = zEm*zEp*(1.0d0-zF)+2.0d0*zEp-2.0d0*zEm*zF
          zRMM =-zEm*zEp*(1.0d0-zF)/zDen
          zTMM = 2.0d0*zexp(zi*zzalpha1h0*loh0)*(zEp-zEm)/zDen
          zAMM =-2.0d0*zF*zEm/zDen
          zBMM = 2.0d0*zEp/zDen

          absRMM = cdabs(zRMM)
          absTMM = cdabs(zTMM)

          if (dabs(absRMM**2+absTMM**2-1.0d0).gt.1.0d-10) then 
            print *,'Energy'
          end if

          end if

c**********************************************************************
c**** Toggle calculation of the surface elevation                  ****
c**********************************************************************

          if (ntheta0pts.eq.1.and.nk0h0pts.eq.1.and.loh0.gt.0.0d0) then

            do 120 xoh0Loop=0,200

              xoh0 = (-1.0d0+3.0d0*dble(xoh0Loop)/200.0d0)*loh0

              do 110 yoh0Loop=0,200

                yoh0 = (-1.0d0+3.0d0*dble(yoh0Loop)/200.0d0)*loh0

                if (xoh0.lt.0.0d0) then

                  zeta = zexp(zi*betah0*yoh0)*
     +            (zexp(zi*alpha0h0*xoh0)+zRMM*zexp(-zi*alpha0h0*xoh0))

                else if (xoh0.gt.loh0) then

                  zeta = zexp(zi*betah0*yoh0)*
     _            zTMM*zexp(zi*alpha0h0*(xoh0-loh0))

                else 

                  zeta = zexp(zi*betah0*yoh0)*
     +                  (zAMM*zexp(zi*zalpha1h0*xoh0)
     +                   +zBMM*zexp(zi*zzalpha1h0*xoh0))

                end if

                write (funitRes,'(4(1x,g10.4))') 
     +                    xoh0,yoh0,dble(zeta),dimag(zeta)

 110          continue

              write (funitRes,'('''')') 

 120        continue

          else

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

            write (*,'(9(1x,g10.4))') 
     +      k0h0,180.0d0*theta0/pi,absRp,absTp,absRm,absTm,absRMM,absTMM
     +          ,theta1
            write (fUnitRes,'(9(1x,g12.6))')  
     +      k0h0,180.0d0*theta0/pi,absRp,absTp,absRm,absTm,absRMM,absTMM
     +          ,theta1

          end if

 998    continue

 999  continue

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

      end
