!!-----------------------------------------------------------------
!! Author: Yaxun Tang, tang@sep.stanford.edu
!!-----------------------------------------------------------------
module wemva_glocal_mod
  use  base_struct_mod
  use  base_progrd_mod
  use  base_outgrd_mod
  use  base_velgrd_mod
  implicit none

contains

  subroutine wemva_glocal(pgrid, ogrid, vgrid, dgrid, igrid, ish)
    type(pro_grid) :: pgrid
    type(out_grid) :: ogrid
    type(vel_grid) :: vgrid
    type(dat_grid) :: dgrid
    type(img_grid) :: igrid
    type(hes_grid) :: hgrid
    integer        :: ish

    integer        :: ishx, ishy
    integer        :: ogrid_extr_step
    real           :: ogrid_hz_max

    if (dgrid%hdrsou) then
       dgrid%sx = dgrid%hsou(1,ish)
       dgrid%sy = dgrid%hsou(2,ish)
    else
       ishx = mod(ish-1, dgrid%asx%n)+1
       ishy = int(1. + real(ish-1)/real(dgrid%asx%n))
       dgrid%sx = dgrid%asx%o + (ishx-1)*dgrid%asx%d
       dgrid%sy = dgrid%asy%o + (ishy-1)*dgrid%asy%d
    end if

    if (dgrid%offset) then
       dgrid%agx%o = dgrid%gxmin_save + dgrid%sx
       dgrid%agy%o = dgrid%gymin_save + dgrid%sy
    end if

    call base_progrd_gen3d(pgrid, dgrid, vgrid)
    call base_outgrd_gen3d(ogrid, igrid, pgrid)
    call base_velgrd_gen3d(vgrid, pgrid)

    ogrid_hz_max     = ogrid%ahz%o + (ogrid%ahz%n-1)*ogrid%ahz%d
    ogrid_extr_step  = 1.5 + (ogrid_hz_max-0.)/ogrid%ahz%d

    pgrid%nzstep = pgrid%nzstep + ogrid_extr_step
    pgrid%nzstep = min(pgrid%az%n, pgrid%nzstep)

    return
  end subroutine

  subroutine wemva_wfs_init(wfds, csou, dgrid, pgrid)
    complex, dimension(:,:,:,:) :: wfds
    complex, dimension(:,:)     :: csou
    type(dat_grid)              :: dgrid 
    type(pro_grid)              :: pgrid

    integer                     :: ix, iy, iws, itr
    integer                     :: indx, indy
    real                        :: gx, gy
    real                        :: wtx, wty
    real                        :: PI

    PI = 3.14159265359

    wfds = cmplx(0., 0.)

    if (dgrid%comsou) then
       do iy=1, dgrid%agy%n
          gy   = dgrid%agy%o + (iy-1)*dgrid%agy%d
          indy = 1.5 + (gy - pgrid%ay%o)/pgrid%ay%d

          do ix=1, dgrid%agx%n
             gx   = dgrid%agx%o + (ix-1)*dgrid%agx%d
             indx = 1.5 + (gx - pgrid%ax%o)/pgrid%ax%d
             itr  = ix + (iy-1)*dgrid%agx%n
             wfds(indx,indy,1,:) = csou(itr, :)
          end do
       end do 

    else 
       !!.. potential error, should use isotropic smoothing
       do iy=max(pgrid%index_datpro_sy-2,1), min(pgrid%index_datpro_sy+2,pgrid%ay%n)
          wty = 0.5 + 0.5 * cos(PI * real(iy-pgrid%index_datpro_sy)/3)
          do ix=max(pgrid%index_datpro_sx-2,1), min(pgrid%index_datpro_sx+2,pgrid%ax%n)
             wtx = 0.5 + 0.5 * cos(PI * real(ix-pgrid%index_datpro_sx)/3)
             wfds(ix,iy,1,:) = wtx*wty*csou(1, :)
          end do
       end do

    end if
    return
  end subroutine

  subroutine wemva_wfr_init(wfdr, crec, dgrid, pgrid)
    complex, dimension(:,:,:,:) :: wfdr
    complex, dimension(:,:)     :: crec
    type(dat_grid)              :: dgrid
    type(pro_grid)              :: pgrid

    integer                     :: ix, iy, itr
    integer                     :: indx, indy 
    real                        :: gx, gy

    wfdr = cmplx(0., 0.)

    do iy=1, dgrid%agy%n
       gy = dgrid%agy%o + (iy-1)*dgrid%agy%d
       indy = 1.5 + (gy - pgrid%ay%o)/pgrid%ay%d
       do ix=1, dgrid%agx%n
           gx = dgrid%agx%o + (ix-1)*dgrid%agx%d
           indx = 1.5 + (gx - pgrid%ax%o)/pgrid%ax%d
           itr  = ix + (iy-1)*dgrid%agx%n
           wfdr(indx,indy,1,:) = crec(itr, :)
       end do
    end do

    return
  end subroutine
end module














