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

  use wei_util
  use wei_process
  use wei_obj_S
  use wei_kmap
  !---------------
  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

    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 camdsr_init
  !----------------------------------------------------------------
!!$=head1 NAME
!!$
!!$weicam_init -initialize  common azimuth migration
!!$
!!$=head1 SYNOPSIS
!!$
!!$ierr= weicam_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$Common Azimuth Migration Initialization operator
!!$
!!$Grabs kxko=2. and nkxko=2001 from par
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<weicam> L<weiwem_init>
!!$
!!$=cut
!!$
  subroutine comaz_init(n)
    integer, intent(in) :: n
    integer             :: er
    integer             :: i,m
    real                :: kxko

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

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

    !! init SSR
    call from_param("nkxko",nkxko,2001)
    call from_param( "kxko", kxko,2.)
    call camdsr_init(kxko,nkxko)
    call  comaz_init(     nkxko)

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weicam_init
!!$=head1 NAME
!!$
!!$weicam - one depth step of common azimuth migration
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weicam(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
!!$
!!$Common Azimuth Migration  operator
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$sepweif90
!!$
!!$=head1 SEE ALSO
!!$
!!$ L<wei_cam_init> L<weiwem>
!!$
!!$=cut
  !----------------------------------------------------------------
  integer function weicam_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    :: khy,kzs,kzr
    real    :: ksy,kry
    real    :: kss,krr
#ifdef DEBUG
    character(len=128), parameter :: unit='weicam'
    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 imy=1,pro%amy%n
 do ihx=pro%ihxb_beg(izs),pro%ihxb_end(izs)
             do imx=1,pro%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) * (kzs-kzr)/(kzs+kzr)

                ksy = kmap%kmy(imy)- khy
                kry = kmap%kmy(imy)+ khy

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

                is = min(int(1 + kss/kos / dkxko) , nkxko)
                ir = 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 weicam_fast
  !----------------------------------------------------------------
  integer function weicam(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    :: khy,kzs,kzr
    real    :: ksy,kry
    real    :: kss,krr
#ifdef DEBUG
    character(len=128), parameter :: unit='weicam'
    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 !! nhy should be 1
       do imy=1,pro%amy%n
          do ihx=1,pro%ahx%n
             do imx=1,pro%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) * (kzs-kzr)/(kzs+kzr)

                ksy = kmap%kmy(imy)- khy
                kry = kmap%kmy(imy)+ khy

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

                is = min(int(1 + kss/kos / dkxko) , nkxko)
                ir = 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 weicam
  !----------------------------------------------------------------
end module wei_cam
