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

  implicit none

  real,private                            :: dkxko, cdkxko,idkzko,rdkzko
  integer     ,private                    :: nkxko
  real                            ::  kxko
  real,    dimension(:), pointer,private :: ckzko,exp_kzko,ikzko
  complex, dimension(:), pointer,private ::  kzko, tkzko
  complex, pointer, private :: w_extrap_1(:,:,:,:,:),w_extrap_2(:,:,:,:,:)
  real, dimension(:) , pointer,private :: khy(:,:)
contains
!----------------------------------------------------------------
!!$=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<weicamo> L<weiwem_init>
!!$
!!$=cut
!!$
  subroutine comaz_init(max,n)
    integer, intent(in) :: n
    real,    intent(in) :: max
    integer             :: er
    integer             :: i,m
    real                :: kxko

    allocate(khy(wsep%n(1),wsep%n(3)))
    call add_mem(size(khy))
    allocate(ckzko(n),stat=er)
    call add_mem(size(ckzko))
    if(er/=0) call seperr('cant allocate ckzko')


    cdkxko = max/(n-1)
    m=n/max

    do i=1,m
       kxko     = 1. * (i-1)/(m-1) !! kx/ko=0...1
       ckzko(i) = sqrt(1-kxko**2)
    end do
    do i=m+1,n
       ckzko(i) = 0
    end do

    ckzko = .5 * ckzko
  end subroutine comaz_init

  integer function weicam_init() result(st)
    integer :: er


    !! init SSR
    call from_param("nkxko",nkxko,2001)
    call from_param( "kxko", kxko,2.)
    call dsr_init(kxko,nkxko,kzko,tkzko,exp_kzko,ikzko,dkxko,nkxko,idkzko,rdkzko)
    call  comaz_init(kxko,nkxko)


    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_caminit> L<weiwem>
!!$
!!$=cut
  integer function weicam(iws,izs,ifk,ith,wfld) result(st)
    integer, intent(in) ::iws,izs,ifk,ith
    complex, dimension(:,:,:,:,:), pointer :: wfld
    real    :: w
    real    :: inv_norm_ko,inv_norm_ko2
    real    :: r_inv_norm_ko,i_inv_norm_ko
    integer :: imx,imy,ihx,ihy 
    integer :: is,ir,i,sign
    complex :: ikz
    real    :: kzs,kzr,k_y
    real    :: ksy,kry,k,j,kmy
    real    :: kss,krr,ko,fract,iko2,slow_0,dstep
    real,pointer :: kx(:,:),ky(:,:)
    type(sep3d) :: wsep


    call grab_current_wave_space(wsep)
    call grab_current_omega(w)
   call calc_comp_depend_vars(w,slow_0,ko,kx,ky,sign,izs,ifk)
    call calc_sign_depend_vars(fract,tkzko,dstep,kzko,ko)
  call calc_khy(ko,inv_norm_ko,izs)

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

                ikz= ko*tkzko(i)

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

             end do
          end do
       end do

    end do




    st=OK
  end function weicam
  !----------------------------------------------------------------

subroutine calc_khy(ko,inv_norm_ko,izs)
integer :: is,ir,ihx,imx,izs
real :: kzs,kzr,ko,inv_norm_ko
 do ihx=1,wsep%n(3)
   do imx=1,wsep%n(1)
        is = min(int(1 + abs(kmap%ksx(imx,ihx)*inv_norm_ko)) , nkxko)
        ir = min(int(1 + abs(kmap%krx(imx,ihx)*inv_norm_ko)) , nkxko)
        kzs = ko*ckzko(is)
        kzr = ko*ckzko(ir)
        if((kzs > 0.) .and. (kzr >0.))then 
           khy(imx,ihx) =(kzr-kzs)/(kzs+kzr)
        else
           khy(imx,ihx) = 0.
       end if
   end do
 end do
end subroutine



end module wei_cam
