      module wei_wem
  !! Wave-Equation Migration
  !! full 3-D prestack

  use wei_util
  use wei_process
  use wei_obj_S
  use wei_kmap
  !---------------
  implicit none

  real,                           private :: dkxko
  integer,                        private :: nkxko
  real,                           private ::  kxko

  complex, dimension(:), pointer, private ::  kzko, tkzko
contains
  !----------------------------------------------------------------
  subroutine wemdsr_init(max,n)
    real,    intent(in) :: max
    integer, intent(in) :: n
    integer             :: er
    integer             :: i,m
    real                :: kxko

    allocate( kzko(n),stat=er)
    if(er/=0) call seperr('cant allocate  kzko')
    allocate(tkzko(n),stat=er)
    if(er/=0) call seperr('cant allocate tkzko')

    dkxko = max/(n-1)
    nkxko = n
    m=n/max

    !!                    ____________
    !!   kz         1    /     | kx |2          kx
    !! i -- dz = -i -   /  1 - | -- |  dz   0 < -- < 1 
    !!   ko         2 \/       | ko |           ko
    !!
    do i=1,m
       kxko    = 0. +     1.0 * (i  -1)/(  m-1) !! kx/ko=0...1
       kzko(i) = - (0,+1) * sqrt(1-kxko**2)
    end do

    !!                    ____________
    !!   kz         1    / | kx |2              kx
    !! i -- dz = -  -   /  | -- | - 1  dz   1 < -- < max 
    !!   ko         2 \/   | ko |               ko
    !!
    do i=m+1,n
       kxko    = 1. + (max-1) * (i-m-1)/(n-m-1) !! kx/ko=1...max
       kzko(i) = - (+1,0) * sqrt(kxko**2-1)
    end do

    kzko = .5 * kzko
  end subroutine wemdsr_init
  !----------------------------------------------------------------
!!$=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
#ifdef DEBUG
    character(len=128), parameter :: unit='weiwem_init'
    if(pro%debug) call in(unit)
#endif

    !! init DSR
    call from_param("nkxko",nkxko,2001)
    call from_param( "kxko", kxko,2.)
    call wemdsr_init(kxko,nkxko)







     
#ifdef DEBUG
    if(pro%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_fast(iws,izs,ifk,ith,wfld) result(st)
    integer, intent(in) ::iws,izs,ifk,ith
    complex, dimension(:,:,:,:,:), pointer :: wfld
    real    :: w, kos,kor
    integer :: imx,imy,ihx,ihy 
    integer :: is,ir
    complex :: ikz


    real    :: kss,krr
#ifdef DEBUG
    character(len=128), parameter :: unit='weiwem'
    if(pro%depar) call in(unit)
#endif

    w  = pro%bw%o + (iws-1)*pro%bw%d
    kos = 2* slo%ssz(izs,ifk) * w
    kor = 2* slo%srz(izs,ifk) * w

    if(pro%opsign==goB2T) then
       tkzko=conjg(kzko)
    else
       tkzko=kzko
    end if

    
  if(pro%ihyb_beg(izs) >1) wfld(:,:,:,1:pro%ihxb_beg(izs)-1,iws)=0.
  if(pro%ihyb_end(izs) < pro%ahy%n) wfld(:,:,:,pro%ihyb_end(izs)+1:pro%ahy%n,iws)=0.
  if(pro%ihxb_beg(izs) >1) wfld(:,:,1:pro%ihxb_beg(izs)-1,pro%ihyb_beg(izs):pro%ihyb_end(izs),iws)=0.
  if(pro%ihxb_end(izs) < pro%ahx%n) wfld(:,:,pro%ihxb_end(izs)+1:pro%ahx%n,pro%ihyb_beg(izs):pro%ihyb_end(izs),iws)=0.
    do ihy=pro%ihyb_beg(izs),pro%ihyb_end(izs)
      do ihx=pro%ihxb_beg(izs),pro%ihxb_end(izs)
          do imy=1,pro%amy%n
             do imx=1,pro%amx%n



                kss = sqrt(kmap%ksx(imx,ihx)**2 + kmap%ksy(imy,ihy)**2)
                krr = sqrt(kmap%krx(imx,ihx)**2 + kmap%kry(imy,ihy)**2)

                is =max(1, min(int(1 + kss/kos / dkxko) , nkxko))
                ir =max(1, min(int(1 + krr/kor / dkxko) , nkxko))

                ikz= kos*tkzko(is) + kor*tkzko(ir)

                !! kz(w) = SSR(kos,|ks|) + SSR(kor,|kr|)
                wfld(imx,imy,ihx,ihy,iws) = &
                wfld(imx,imy,ihx,ihy,iws) * &
                cexp( ikz * pro%dzb(izs))

             end do
          end do
       end do
    end do

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
    st=OK
  end function weiwem_fast
  !----------------------------------------------------------------
  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
    integer :: imx,imy,ihx,ihy 
    integer :: is,ir
    complex :: ikz


    real    :: kss,krr
#ifdef DEBUG
    character(len=128), parameter :: unit='weiwem'
    if(pro%depar) call in(unit)
#endif

    w  = pro%bw%o + (iws-1)*pro%bw%d
    kos = 2* slo%ssz(izs,ifk) * w
    kor = 2* slo%srz(izs,ifk) * w

    if(pro%opsign==goB2T) then
       tkzko=conjg(kzko)
    else
       tkzko=kzko
    end if

    do ihy=1,pro%ahy%n
       do ihx=1,pro%ahx%n
          do imy=1,pro%amy%n
             do imx=1,pro%amx%n







                kss = sqrt(kmap%ksx(imx,ihx)**2 + kmap%ksy(imy,ihy)**2)
                krr = sqrt(kmap%krx(imx,ihx)**2 + kmap%kry(imy,ihy)**2)

                is =max(1, min(int(1 + kss/kos / dkxko) , nkxko))
                ir =max(1, min(int(1 + krr/kor / dkxko) , nkxko))

                ikz= kos*tkzko(is) + kor*tkzko(ir)

                !! kz(w) = SSR(kos,|ks|) + SSR(kor,|kr|)
                wfld(imx,imy,ihx,ihy,iws) = &
                wfld(imx,imy,ihx,ihy,iws) * &
                cexp( ikz * pro%dzb(izs))

             end do
          end do
       end do
    end do

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
    st=OK
  end function weiwem
  !----------------------------------------------------------------
end module wei_wem
