# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module wem_mod
  !! Wave-Equation propagation
  !! full 3-D prestack SSR

  use util_mod
  use down_mod
  use external_mod
  use S_mod
  use kmap_mod
  use sep
  !---------------
  implicit none

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

  complex, dimension(:,:), allocatable, private :: ikz
  real, dimension(:), allocatable, private, save :: amp,phase

contains
  !----------------------------------------------------------------
  !!make amp and phase lookup tables
  integer function wem_init() result(st)
    integer :: er
    integer :: i,m,nmax
    real :: kmax=2.,kz,a





    call from_param("nkxko",nkxko,2001)
    call from_param('maxang',maxang,50.)
    down%maxang=maxang

    allocate( amp (nkxko),stat=er)
    if(er/=0) call seperr('cant allocate  amp: wem.F90')
    allocate( phase(nkxko),stat=er)
    if(er/=0) call seperr('cant allocate  phase: wem.F90')
    allocate( ikz (nkxko,run%o_nth),stat=er)
    if(er/=0) call seperr('cant allocate  ikz: wem.F90')

    dkxko = kmax/(nkxko-1)
    m=nkxko/2.

    !! ____________
    !! 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
       phase(i)=-sqrt(1-kxko**2)!*down%az%d
       amp(i)=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,nkxko
       kxko = 1. + (kmax-1.)*real(i-m-1)/real(nkxko-m-1) !! kx/ko=1...max
       amp(i) = exp(-sqrt(kxko**2-1))
       phase(i)=0.
! kzko(i) = - (+1.,0.) * sqrt(kxko**2-1)
    end do
!
! !! taper the high propagation angles if so desired (set parameter maxang<90.)
! kz=sin(maxang*atan(1.)/45.)
! m=kz/dkxko
! nmax=min(m+nint(nkxko/9.),nkxko)
! do i=m,nmax
! amp(i)=cos(atan(1.)*2.*real(i-m)/real(nmax-m))**4
! end do
! do i=nmax+1,nkxko
! amp(i)=0.
! end do
!
    if(down%opsign==goB2T) then
       do i=1,nkxko
          phase(i)=-phase(i)
       end do
    end if
    a=4*(size(amp)+size(phase)+2*size(ikz))
    down%mem=down%mem+a
    st=OK






  end function wem_init

!----------------------------------------------------------------
  integer function wem(izs,iws,ith,ivel,wfld) result(st)
    integer, intent(in) ::iws,izs,ivel,ith
    complex, dimension(:,:,:) :: wfld
    real :: w, ko, k
    integer :: ix,iy,i





    w = down%bw%o + (iws-1)*down%aw%d! omega
    ko = slo%sz(izs,ivel) * w!slowness*omega

    do i=1,nkxko !convert lookup tables to include dz,w
       k=ko*phase(i)*down%az%d !! dz applied during wem_init above
       ikz(i,ith)=cmplx(amp(i)*cos(k),amp(i)*sin(k))
    end do
    ky_loop:do iy=1,down%amy%n
       kx_loop:do ix=1,down%amx%n
          k=sqrt(kmap%kx(ix)*kmap%kx(ix) + kmap%ky(iy)*kmap%ky(iy))
          i= max(1, min(int(1 + k/ko / dkxko) , nkxko))
          wfld(ix,iy,iws)=wfld(ix,iy,iws)*ikz(i,ith)
       end do kx_loop
    end do ky_loop




    st=OK
  end function wem

end module wem_mod

! ! ikz = ko*kzko(i)*var%dzb(izs)
!
!! bin(ix,ith)=real(ikz)
! ! sc(ix,ith)=cmplx(cos(aimag(ikz)),sin(aimag(ikz)))
! end do kx_loop
! ! do ix=1,size(bout,1)
! ! bout(ix,ith)=exp(bin(ix,ith))
! ! end do
! ! do ix=1,size(bout,1)
! wfld(ix,iy,iws)=wfld(ix,iy,iws)*tkzko(i)
! ! end do
!! end do ky_loop
