!! wavefield extrapolation via Split-Step Fourier for one frequency and one depth step
!! should initialize in advance for each frequency and depth step 
module wei_source_boundary_mod
  use  wei_myfft_mod
  use  wei_type_mod
  use  wei_vel_mod
  use  sep
  implicit none
  integer, external                      :: sreed_window
  real, dimension(:)  , pointer, private :: kx2
  real, dimension(:)  , pointer, private :: ky2 
  real, dimension(:,:), pointer, private :: vmod_sb
  real, dimension(:)  , pointer, private :: vref_sb
  real                         , private :: kz2_max
contains
!=========================================================================================
  integer function wei_source_boundary_init(vgrid, pgrid, kx2_in, ky2_in, maxang_in) result(stat)
    type(vel_grid)              :: vgrid
    type(pro_grid)              :: pgrid
    real, dimension(:), pointer :: kx2_in
    real, dimension(:), pointer :: ky2_in
    real                        :: maxang_in
    integer, dimension(3)       :: n_cube, j_wind, n_wind, f_wind
    integer                     :: n_dim_cube, err

    n_dim_cube = 3
    n_cube = (/vgrid%ax%n, vgrid%ay%n, vgrid%az%n /)
    j_wind = (/1, 1, 1/)
    n_wind = (/pgrid%ax%n, pgrid%ay%n, 1/)
    f_wind = (/vgrid%index_provel_xmin-1, vgrid%index_provel_ymin-1, 0/)

    kx2     => kx2_in
    ky2     => ky2_in


    kz2_max =  cos(maxang_in*4.*atan(1.)/180.)**2

    allocate(vmod_sb(pgrid%ax%n, pgrid%ay%n))
    allocate(vref_sb(1))

    vmod_sb =  0.
    vref_sb =  0.

    err = sreed_window("vel",n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,vmod_sb)
    if (err .lt. 0) call seperr("wei_source_boundary_init: trouble reading velocity file")    

    stat = wei_vel_oper(vmod_sb, vref_sb)
    write(0,*)"vmod_sb", maxval(vmod_sb), minval(vmod_sb)
    write(0,*)"vref_sb", vref_sb(1)
  end function
!=========================================================================================
  integer function wei_source_boundary_clos() result(stat)
    deallocate(vmod_sb)
    deallocate(vref_sb)
    nullify(kx2)
    nullify(ky2) 
    stat = 0
  end function
!=========================================================================================
  integer function wei_source_boundary_oper(aws, wfld_d) result(stat)
    type(axis)                         :: aws
    complex, dimension(:,:,:) :: wfld_d

    integer                            :: nkx, nky
    integer                            :: ikx, iky, iws
    real                               :: lamda, kz, kz2, ws, wov, wov2
    complex                            :: shaping

    nkx = size(wfld_d, 1) 
    nky = size(wfld_d, 2)
    do iws=1, aws%n
       ws   = aws%o + (iws-1)*aws%d
       wov  = ws / vref_sb(1)
       wov2 = wov * wov
       call FFT2DX(wfld_d(:,:,iws),  1)
       do iky=1, nky
          do ikx=1, nkx
             kz2 = 1.0 - kx2(ikx)/wov2 - ky2(iky)/wov2
             if (kz2 >= kz2_max) then
                kz = sqrt(kz2)
                lamda = wov * kz
                if (abs(lamda) > 0.000001) then
                   wfld_d(ikx, iky, iws) = wfld_d(ikx, iky, iws) * cmplx(0., 0.5 / lamda )
                end if
             end if
          end do
       end do
       call FFT2DX(wfld_d(:,:,iws), -1)
    end do      
    stat = 0
  end function
!==================================================================================
!  subroutine wei_source_boundary_freq(ws, vref_sb, wfld_slice_d)
!    real                             :: ws
!    real                             :: vref_sb
!    complex, dimension(:,:), pointer :: wfld_slice_d
!    real                             :: ws, wov, wov2, kz2, kz, lamda
!    wov  = ws / vref_sb
!    wov2 = wov * wov
!    do iky=1, size(ky2)
!       do ikx=1, size(kx2)
!          kz2 = 1.0 - kx2(ikx)/wov2 - ky2(iky)/wov2
!          if (kz2 >= 0. ) then
!             kz = sqrt(kz2)
!          else
!             kz = sqrt(-kz2)
!          end if
!          lamda   = wov * kz
!          wfld_slice_d(ikx, iky) = wfld_slice_d(ikx, iky) * cmplx(0., 0.5 / lamda)
!       end do
!    end do
!  end subroutine
!====================================================================================
!  subroutine wei_source_boundary_spac(ws, vref_sb, vmod_sb, wfld_slice_d)
!    do iy=1, size(ky2)
!       do ix=1, size(kx2)
!          lamda = 0.5 * ws * (1./vmod_sb(ix,iy) - 1./vref_sb)
!          wfld_slice_d(ix, iy) = 
!       end do
!    end do
!  end subroutine
end module














