module wei_wem
 use wei_fk_share_mod

  !! Wave-Equation Migration
  !! full 3-D prestack
  implicit none

  real,                           private,save :: dkxko
  integer,                        private,save :: nkxko
  real,                           private,save ::  kxko,rdkzko,idkzko

  real   , dimension(:), pointer, private,save ::  exp_kzko,ikzko,ckzko,bin(:),bout(:)
  complex, dimension(:), pointer, save ::  wemkzko,sc(:)
  complex, dimension(:), pointer, private,save :: tkzko

contains



  !----------------------------------------------------------------
!!$=head1 NAME
!!$
!!$weiwem_init -initialize  wave equation migration 
!!$
!!$=head1 SYNOPSIS
!!$
!!$ierr= weiwem_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$Wave equation migration operator
!!$
!!$Grabs kxko=2. and nkxko=2001 from par
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weiwem> L<weiwem_init>
!!$
!!$=cut
!!$
  integer function weiwem_init() result(st)
    integer :: er
    integer ::i,imax
    real :: hxmax,hymax
    type(sep3d) :: wsep
#ifdef DEBUG
    character(len=128), parameter :: unit='weiwem_init'
    if(my_debug()) call in(unit)
#endif

    !! init DSR
    call from_param("nkxko",nkxko,2001)
    call from_param( "kxko", kxko,2.)
    call dsr_init(kxko,nkxko,wemkzko,tkzko,exp_kzko,ikzko,dkxko,nkxko,idkzko,rdkzko)
    call grab_current_wave_space_pad(wsep)
    allocate(bin(wsep%n(1)))
    allocate(bout(wsep%n(1)))
    allocate(sc(wsep%n(1)))
   


#ifdef DEBUG
    if(my_debug()) call out(unit)
#endif
    st=OK
  end function weiwem_init


!!$=head1 NAME
!!$
!!$weiwem - one depth step of  wave equation migration
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weiwem(iws,izs,ifk,ith,wfld)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item iws - integer
!!$
!!$  current frequency slice
!!$
!!$=item izs - integer
!!$
!!$  current depth slice
!!$
!!$=item ifk - integer
!!$
!!$  current fk operation
!!$
!!$=item ith - integer
!!$
!!$  current thread
!!$
!!$=item wfld - complex(:,:,:,:,:)
!!$
!!$  wavefield
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$Wave Equation Migration operator
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weiwem_init> L<weicam>
!!$
!!$=cut
  !----------------------------------------------------------------
  integer function weiwem(iws,izs,ifk,ith,wfld) result(st)
    integer, intent(in) ::iws,izs,ifk,ith
    complex, dimension(:,:,:,:,:), pointer :: wfld
    real    :: w, kos,kor,ko
    integer :: imx,imy,ihx,ihy 
    integer :: is,ir,i,iz,sign
    complex :: ikz
    real,pointer :: kx(:,:),ky(:,:)
    real    :: kss,krr,k,slow_0,fract,dstep,re,im
    type(sep3d) :: wsep
#ifdef DEBUG
    character(len=128), parameter :: unit='weiwem'
    if(my_debug()) call in(unit)
#endif

!  if(.not. timer_start("FK")) call seperr("trouble starting timer")
  call grab_current_omega(w)
  call grab_current_wave_space_pad(wsep)
  call grab_current_dstep(dstep)




   call calc_comp_depend_vars(w,slow_0,ko,kx,ky,sign,izs,ifk)
    call calc_sign_depend_vars(fract,tkzko,dstep,wemkzko,ko)

    do ihy=1,size(wfld,4)
      do ihx=1,size(wfld,3)
          do imy=1,wsep%n(2)
             do imx=1,wsep%n(1)
                k = sqrt(kx(imx,ihx)**2 + ky(imy,ihy)**2)
                i =max(1, min(int(1 + k/ko / dkxko) , nkxko))

                ikz= ko*tkzko(i) *dstep
                bin(imx)= real(ikz)
                sc(imx)=cmplx(cos(aimag(ikz)),sin(aimag(ikz)))


                !! kz(w) = SSR(kos,|ks|) 
!                wfld(imx,imy,ihx,ihy,iws) = &
!                  wfld(imx,imy,ihx,ihy,iws) * cexp( ikz * dstep)
!                  wfld(imx,imy,ihx,ihy,iws) * cexp( ikz * dstep)

             end do
#ifdef SEP_MKL
            call vsexp(size(wfld,1),bin,bout)
#else
             bout=exp(bin)
#endif
                wfld(:,imy,ihx,ihy,iws) = &
                  wfld(:,imy,ihx,ihy,iws) * bout*sc

          end do
       end do

    end do

!  if(.not. timer_stop("FK")) call seperr("trouble stopping timer")

#ifdef DEBUG
    if(my_debug()) call out(unit)
#endif
    st=OK
  end function weiwem

end module wei_wem
