module wei_taper_mod
  use sep
  implicit none
  integer                        , private :: nxpad_beg, nxpad_end
  integer                        , private :: nypad_beg, nypad_end 
  integer                        , private :: taper_type
  real, parameter                , private :: PI = 3.14159265359
contains

  integer function wei_taper_2d_init(taper_type_in,  &
                   nxpad_beg_in, nxpad_end_in, nypad_beg_in, nypad_end_in) result(stat)
    integer       :: taper_type_in
    integer       :: nxpad_beg_in, nxpad_end_in
    integer       :: nypad_beg_in, nypad_end_in 
    nxpad_beg = nxpad_beg_in
    nxpad_end = nxpad_end_in
    nypad_beg = nypad_beg_in
    nypad_end = nypad_end_in
    taper_type = taper_type_in

    stat = 0 
  end function
!===========================================================
  integer function wei_taper_2d_oper(wfld, taper_pct) result(stat)
    complex, dimension(:,:,:), pointer :: wfld
    real                               :: taper_pct

    if (taper_type == 1) then
       call wei_rectangular_taper(wfld, taper_pct)
    else if (taper_type == 2) then
       call wei_circular_taper(wfld, taper_pct)
    else
       call seperr("you have to specify the type of wavefield taper")
    end if
    stat = 0
  end function
!================================================================
  integer function wei_taper_2d_clos() result(stat)
    nxpad_beg = 0
    nxpad_end = 0
    nypad_beg = 0
    nypad_end = 0
    taper_type = 0
    stat = 0
  end function

!--------------------------------------------------------
  subroutine wei_rectangular_taper(wfld, taper_pct) 
    complex, dimension(:,:,:), pointer :: wfld
    real                               :: taper_pct
    real, dimension(:), allocatable    :: taper_x, taper_y
    integer                            :: taplen, err
    integer                            :: ix, iy, iw, nx, ny, nw
    nx = size(wfld, 1) - nxpad_beg - nxpad_end
    ny = size(wfld, 2) - nypad_beg - nypad_end
    nw = size(wfld, 3)
    taplen = taper_pct * nx

    if (nx > 1) then
       allocate(taper_x(nx), stat=err)
       if (err /= 0) call seperr("trouble allocating taper_x")
       call hamming_window_1d(taper_x, taplen)
       do iy=1, ny
          do iw=1, nw
             wfld(:, iy, iw) = wfld(:, iy, iw) * taper_x
          end do
       end do
       deallocate(taper_x)
    end if


    if (ny > 1) then
       allocate(taper_y(ny), stat=err)
       if (err /= 0) call seperr("trouble allocating taper_y")
       call hamming_window_1d(taper_y, taplen)
       do ix=1, nx
          do iw=1, nw
             wfld(ix, :, iw) = wfld(ix, :, iw) * taper_y
          end do
       end do
    end if

  end subroutine
  !----------------------------------------------
  ! 2D Hamming window taper
  !----------------------------------------------
  subroutine wei_circular_taper(wfld, taper_pct)
    complex, dimension(:,:,:), pointer :: wfld
    real                               :: taper_pct
    integer                            :: taplen
    integer                            :: nx, ny, nw, ix, iy, iw
    real                               :: x_orig, y_orig
    real                               :: taper_radius_min
    real                               :: taper_radius_max
    real                               :: taper, dist


    nx = size(wfld, 1) - nxpad_beg - nxpad_end
    ny = size(wfld, 2) - nypad_beg - nypad_end
    nw = size(wfld, 3)
    taplen = taper_pct * min(nx, ny)

    x_orig = nx * 0.5
    y_orig = ny * 0.5

    taper_radius_min = min(nx, ny) * 0.5 - taplen
    taper_radius_max = min(nx, ny)

    do iy=1, ny
       do ix=1, nx
          dist = sqrt((ix-x_orig)**2 + (iy-y_orig)**2)
          if (dist <= taper_radius_min ) then
             taper = 1.
          else if (dist > taper_radius_min .and. dist <= taper_radius_max) then
             taper = 0.54 + 0.46 * cos(PI*(taper_radius_min-dist)/(taper_radius_max-taper_radius_min))
          else
             taper = 0.
          end if
          do iw=1, nw
             wfld(ix+nxpad_beg, iy+nypad_beg, iw) =  wfld(ix+nxpad_beg, iy+nypad_beg, iw) * taper
          end do
       end do
    end do
  end subroutine
  !----------------------------------------------
  ! Hamming window taper
  !----------------------------------------------
  subroutine hamming_window_1d(taper, length)
    real, dimension(:)          :: taper
    integer, intent(in)         :: length
    integer                     :: n1, n2, n3, n4, i
    n1 = 1
    n2 = length
    n3 = size(taper)-length+1
    n4 = size(taper)
    taper = 1.
    do i=n1, n2
       taper(i) = 0.54 + 0.46 * cos(PI*(i-n1)/(n2-n1)-PI)
    end do
    do i=n3, n4
       taper(i) = 0.54 + 0.46 * cos(PI*(n3-i)/(n4-n3))
    end do
  end subroutine


end module
