module wei_damp_mod
  use sep3d_struct_mod
  use wei_step_mod
  use sep
  implicit none

  contains


  subroutine init_damping(wave)
    type(wspace) :: wave
    integer :: ncx,ncy,nhx,nhy
#ifdef DEBUG
    write(0,*) "IN INIT PADDING"
#endif

     call from_param("npad_mx",ncx,0)            ! How much padding in x
     call from_param("npad_my",ncy,0)            ! How much padding in y
     call from_param("npad_hx",nhx,0)            ! How much padding in x
     call from_param("npad_hy",nhy,0)            ! How much padding in y
     wave%npad=wave%n;
     wave%opad=wave%o;
     wave%opad(1:4)=wave%o(1:4)-&
         (/ncx,ncy,nhx,nhy/)*wave%d(1:4)
     wave%dpad=wave%d
     wave%npad(1:4)=wave%npad(1:4)+2*(/ncx,ncy,nhx,nhy/)

#ifdef DEBUG
    write(0,*) "OUT INIT PADDING"
#endif

  end subroutine


  !BUILD PADDING/EDGE EFFECT
  subroutine build_padding(step,npad)
    type(westep)       :: step          !step object
    integer            :: npad(4)   !padding in x and y
#ifdef DEBUG
    write(0,*) "IN BUILD PADDING"
#endif
    allocate(step%amp_cx(step%wave%npad(1)))
    allocate(step%amp_cy(step%wave%npad(2)))
    allocate(step%amp_hx(step%wave%npad(3)))
    allocate(step%amp_hy(step%wave%npad(4)))
    call build_amp(step%wave%npad(1),npad(1),step%amp_cx)
    call build_amp(step%wave%npad(2),npad(2),step%amp_cy)
    call build_amp(step%wave%npad(3),npad(3),step%amp_hx)
    call build_amp(step%wave%npad(4),npad(4),step%amp_hy)
#ifdef DEBUG
    write(0,*) "OUT BUILD PADDING"
#endif
  end subroutine


  !BUILD AMPLITUDE FUNCTION FOR DAMPING
  subroutine build_amp(nlen,npad,amp)
    integer :: i,nlen,npad,n2
    real :: amp(:),f

#ifdef DEBUG
    write(0,*) "iN BUILD AMP"
#endif
    amp=1.
    if(nlen/=1 .and. npad/=0) then
      do i=1,npad
        f=real(npad+1-i)/real(npad)
        amp(i)=cos(f*atan(1.)*2)
        amp(nlen-i+1)=cos(f*atan(1.)*2)
        write(0,*) "CHECK",i,f,amp(i),amp(nlen-i+1)
      end do
    end if
#ifdef DEBUG
    write(0,*) "OUT BUILD AMP"
#endif
  end subroutine

  !DAMP THE WAVEFIELD
  subroutine damp_wave(wfld,step)
    complex         :: wfld(:,:,:,:,:)
    type(westep)    :: step 
    integer         :: i1,i2,i3,i4,iws
    integer         :: ibcx,ibcy,ibhx,ibhy

#ifdef DEBUG
    write(0,*) "IN DAMP WAVE"
#endif
    iws=step%iw

    ibcx=(step%wave%npad(1)-step%wave%n(1))/2
    ibcy=(step%wave%npad(2)-step%wave%n(2))/2
    ibhx=(step%wave%npad(3)-step%wave%n(3))/2
    ibhy=(step%wave%npad(4)-step%wave%n(4))/2


    if(ibcy/=0) then
      do i4=1,size(wfld,4)
        do i3=1,size(wfld,3)
          do i2=1,ibcy
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_cy(i2)
            end do 
          end do 
        end do 
      end do 
      do i4=1,size(wfld,4)
        do i3=1,size(wfld,3)
          do i2=size(wfld,2)-ibcy+1,size(wfld,2)
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_cy(i2)
            end do 
          end do 
        end do 
      end do 
    end if

    if (ibcx /=0) then
      do i4=1,size(wfld,4)
        do i3=1,size(wfld,3)
          do i2=1,size(wfld,2)
            do i1=1,ibcx
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_cx(i1)
            end do 
          end do 
        end do 
      end do 
      do i4=1,size(wfld,4)
        do i3=1,size(wfld,3)
          do i2=1,size(wfld,2)
            do i1=size(wfld,1)-ibcx+1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_cx(i1)
            end do 
          end do 
        end do 
      end do 
    end if

    if (ibhy /=0) then
      do i4=1,ibhy
        do i3=1,size(wfld,3)
          do i2=1,size(wfld,2)
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_hy(i1)
            end do 
          end do 
        end do 
      end do 
      do i4=size(wfld,4)-ibhy+1,size(wfld,4)
        do i3=1,size(wfld,3)
          do i2=1,size(wfld,2)
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_hy(i1)
            end do 
          end do 
        end do 
      end do 
    end if
    if (ibhx /=0) then
      do i4=1,size(wfld,4)
        do i3=1,ibhx
          do i2=1,size(wfld,2)
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_hx(i3)
            end do 
          end do 
        end do 
      end do 
      do i4=1,size(wfld,4)
        do i3=size(wfld,3)-ibhx+1,size(wfld,3)
          do i2=1,size(wfld,2)
            do i1=1,size(wfld,1)
              wfld(i1,i2,i3,i4,iws)=wfld(i1,i2,i3,i4,iws)*step%amp_hx(i3)
            end do 
          end do 
        end do 
      end do 
    end if



#ifdef DEBUG
    write(0,*) "OUT DAMP WAVE"
#endif
  end subroutine

end module
