  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !! . . Fourier-domain phase shift
  !! . . Set up for 2D
  !! . .
  !! . . Inputs:
  !!          dax  - 1D data (in Fourier domain)
  !!          w    - current angular frequency
  !!          vmin - minimum velocity of extrapolation step
  !!          darg - causal/acausal option (source/receiver wavefield)
  !!          nx   - # of samples in array
  !!          dz   - extrapolation step
  !!          okx  - Spatial Nyquist
  !!          dkx  - sampling in spatial wavenumber
  !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine wemig_phs_correction(dax,w,vmin,darg,nx,dz,okx,dkx)
  
  integer :: ikx,nx
  real    :: w,kr,darg,vmin,kx,w,okx,dkx
  complex :: dax(:) 
  real    :: w_v0,sr2,phsft,dz,a1,a2,b1,b2,wv2

  !! . . Coefficients for filtering
  !! . . Depends on what your accuracy in propagation
  a1=-0.040315157; b1=-0.873981642
  a2=-0.457289566; b2=-0.222691983
  
  w_v0=w/vmin
  wv2=w_v0*w_v0

  do ikx=1,nx

     kx  = okx+(ikx-1.)*dkx
     sr2 = kx*kx/wv2

     !! . . Filter out evanescent energy
     if ( sr2 .gt. 1) then
        dax(ikx) = 0.
     else
        phsft=-dz*w_v0*( &
        sqrt(1.-sr2) - &
        (1. + a1*kx**2/wv2/(1.+b1*((kx)/w_v0)**2) + &
              a2*kx**2/wv2/(1.+b2*((kx)/w_v0)**2) )

        dax(ikx)=dax(ikx)*cmplx(cos(phsft),darg*sin(phsft))
     end if
  end do

end subroutine wemig_phs_correction
