      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 theta00,weight
      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 complex   phi(0:200,0:200)
      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
      theta00 = theta0Start+dble(ntheta0Pts/2)*theta0Step

      delta = delta*pi/180.0d0

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

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

      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**********************************************************************
c**** Don't loop through k0h0                                      ****
c**********************************************************************

      k0h0 = k0h0Start
      k1h1p = dsqrt(h1poh0)*k0h0
      k1h1m = dsqrt(h1moh0)*k0h0

c**********************************************************************
c**** Set field to zero everywhere                                 ****
c**********************************************************************

      do 20 xoh0Loop=0,200
        do 10 yoh0Loop=0,200
          phi(xoh0Loop,yoh0Loop) = dcmplx(0.0d0,0.0d0)
 10     continue
 20   continue

c**********************************************************************
c**** Loop through values of theta0                                ****
c**********************************************************************

      do 998 theta0Loop=1,ntheta0Pts

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

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

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

        do 120 xoh0Loop=0,200

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

          do 110 yoh0Loop=0,200

            yoh0 = (-1.5d0+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

            weight = 2.0d0*dsqrt(pi)*theta0Step*dcos(theta00-theta0)*
     +                    dexp(-(2.0d0*pi*dsin(theta00-theta0))**2)

            phi(xoh0Loop,yoh0Loop) = phi(xoh0Loop,yoh0Loop)+zeta*weight

 110      continue

 120    continue

 998  continue

c**********************************************************************
c**** Write to file                                                ****
c**********************************************************************

      do 160 xoh0Loop=0,200
        xoh0 = (-1.0d0+3.0d0*dble(xoh0Loop)/200.0d0)*loh0
        do 150 yoh0Loop=0,200
          yoh0 = (-2.0d0+3.0d0*dble(yoh0Loop)/200.0d0)*loh0

          write (fUnitRes,'(4(1x,f12.6))') xoh0,yoh0,
     +       dble(phi(xoh0Loop,yoh0Loop)),dimag(phi(xoh0Loop,yoh0Loop))

  150   continue
        write (fUnitRes,'('''')')
  160 continue

      close(funitRes)



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

      end
