module weiop_cam
  !! Wave-Equation Migration
  !! Common-azimuth

  use wei_util
  use wei_pro
  use wei_par
  use wei_S
  use wei_kmap
  use wei_comp
  !---------------
  implicit none

  real,                           private :: dkxko, cdkxko
  integer,                        private :: nkxko
  real,                           private ::  kxko
  real,    dimension(:,:), pointer, private :: ckzko 
  complex, dimension(:,:), pointer, private ::  kzko, tkzko
contains
  !----------------------------------------------------------------
  subroutine camdsr_init(max,n)
    real,    intent(in) :: max
    integer, intent(in) :: n
    integer             :: er
    integer             :: i,m
    real                :: kxko
    integer             :: iomp
    character(len=128), parameter :: unit='camdsr_init'

    call weimem_add(8,n*wpar%nomp,unit)
    allocate( kzko(   n,wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate  kzko')
    call weimem_add(8,n*wpar%nomp,unit)
    allocate(tkzko(   n,wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate tkzko')

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

    do iomp=1,wpar%nomp
       !!                    ____________
       !!   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,iomp) = - (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,iomp) = - (+1,0) * sqrt(kxko**2-1)
       end do

       kzko(:,iomp) = .5 * kzko(:,iomp)
    end do

  end subroutine camdsr_init
  !----------------------------------------------------------------
  subroutine comaz_init(n)
    integer, intent(in) :: n
    integer             :: er
    integer             :: i,m
    real                :: kxko
    integer             :: iomp
    character(len=128), parameter :: unit='comaz_init'

    call weimem_add(8,n*wpar%nomp,unit)
    allocate(ckzko(   n,wpar%nomp),stat=er)
    if(er/=0) call seperr('cant allocate ckzko')

    cdkxko = 1./(n-1)
    ckzko = 0.

    do iomp=1,wpar%nomp
       do i=1,n
          kxko     = 0. + 1.0 * (i-1)/(n-1) !! kx/ko=0...1
          ckzko(i,iomp) = sqrt(1-kxko**2)
       end do
       ckzko(n,iomp) = ckzko(n-1,iomp) !! avoid kzs=0 or kzr=0
       ckzko(:,iomp) = .5 * ckzko(:,iomp)
    end do
  end subroutine comaz_init
  !----------------------------------------------------------------
  integer function weiop_cam3_init() result(st)
    integer :: er
    integer :: mpierr
    character(len=128), parameter :: unit='weiop_cam3_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    !! init DSR
    if(wpar%impi==0) call from_param("nkxko",nkxko,2001)
    if(wpar%nmpi >0) call MPI_BCAST(nkxko, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)
    if(wpar%impi==0) call from_param( "kxko", kxko,2.)
    if(wpar%nmpi >0) call MPI_BCAST( kxko, 1,MPI_REAL   ,0,MPI_COMM_WORLD,mpierr)
    call camdsr_init(kxko,nkxko)
    call  comaz_init(     nkxko)

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_cam3_init
  !----------------------------------------------------------------
  integer function weiop_cam3(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 :: ii
    integer :: is,ir
    complex :: ikz
    real    :: khy,kzs,kzr
    real    :: ksy,kry
    real    :: kss,krr
    real    :: kk, ky
#ifdef DEBUG
    character(len=128), parameter :: unit='weiop_cam3'
    if(wpro%dbg>9) call in(unit)
#endif

    w  = wpro%bw%o + (iws-1)*wpro%bw%d
    ko = 2*w* srp(ith)%rs(izs,ifk)

    kos = 2*w* slo%rs_s(izs,ifk)
    kor = 2*w* slo%rs_r(izs,ifk)

    if(wpro%opsign==weiB2T) then
       tkzko(:,ith)=conjg(kzko(:,ith))
    else
       tkzko(:,ith)=kzko(:,ith)
    end if

    do ihy=1,wpro%ahy%n !! nhy should be 1
       do imy=1,wpro%amy%n
          do ihx=1,wpro%ahx%n
             do imx=1,wpro%amx%n

                is = min(int(1 + abs(kmap%ksx(imx,ihx)/kos) / cdkxko) , nkxko)
                ir = min(int(1 + abs(kmap%krx(imx,ihx)/kor) / cdkxko) , nkxko)
                kzs = kos*ckzko(is)
                kzr = kor*ckzko(ir)
		khy = kmap%kmy(imy) * (kzr-kzs)/(kzs+kzr)

                ky = kmap%kmy(imy) +  srp(ith)%sig * khy

                kk = sqrt(srp(ith)%kx(imx,ihx)**2 + ky**2)
                ii = max(1, min(int(1 + kk/ko / dkxko) , nkxko))
                ikz= ko*tkzko(ii,ith)

                wfld(imx,imy,ihx,ihy,iws) = &
                wfld(imx,imy,ihx,ihy,iws) * &
                cexp( ikz * wpro%az%d)

             end do
          end do
       end do
    end do

#ifdef DEBUG
    if(wpro%dbg>9) call out(unit)
#endif
    st=weiOK
  end function weiop_cam3
  !----------------------------------------------------------------
end module weiop_cam




