module wei_cam_mod
  use wei_kmap_mod
  use wei_table_mod
  use wei_step_mod
  implicit none

  real, pointer,private:: khy(:,:),ckzko(:)
  real, private, save :: cdkxko
  integer,private, save :: nkxko
   


  contains 

  subroutine cam_init(wave,fk)
    type(wspace)   :: wave  ! wavefield description
    type(fktype)   :: fk    ! fk structure
    integer        :: n,m,i
    real           :: kxko,kmax

    allocate(fk%amp(wave%npad(1)))    !amplitude term
    allocate(fk%temp1(wave%npad(1)))  !before exp
    allocate(fk%temp2(wave%npad(1)))  !after  exp
    allocate(khy(wave%npad(1),wave%npad(3)))  !after  exp

    nkxko=2001
    allocate(ckzko(nkxko))
    kmax=2.

    cdkxko = kmax/(nkxko-1)
    m=nkxko/kmax

    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,nkxko
       ckzko(i) = 0
    end do

    ckzko = .5 * ckzko

  end subroutine


  subroutine wei_cam(wfld,step,dsr)
    !Passed in parameters
    complex     :: wfld(:,:,:,:,:)   !5-D wfld object
    real        :: dstep             !step size 
    real        :: slow              !slowness
    real        :: w                 !w   frequency
    type(ctable):: dsr               !dsr lookup table
    type(westep):: step              !current step
    !internal variables
    integer    :: i1,i2,i3,i4       !loop counting
    integer    :: itable            !table lookup
    real       :: k                 !wave number
    complex    :: ikz               !multiplication term
    integer     :: i5                !current i5 slice
    integer     :: ith,i
    real        :: ko,kmy,k_y
    real        :: sign,inv_norm_ko
    complex,  pointer:: amp(:)
    real, pointer:: temp1(:)
    real, pointer:: temp2(:)
    

#ifdef DEBUG
     write(0,*) "IN WEI_WEM"
#endif
   amp=>step%fk%amp
   temp1=>step%fk%temp1
   temp2=>step%fk%temp2
   i5=step%iw
   ith=step%ithread
   ko= 2* step%sref * step%w

   do i=1,size(step%ctable)
     k=ko*step%dstep*dsr%phase(i)
     step%ctable(i)=dsr%amp(i)*cmplx(cos(k),sin(k))
   end do


    if(step%source) then
      sign=-1
    else
      sign=1
    end if

    inv_norm_ko=1./(ko*cdkxko)

    call calc_khy(ko,inv_norm_ko,step%wave)



    do i4=1,size(wfld,4)
      do i3=1,size(wfld,3)
        do i2=1,size(wfld,2)
          kmy=step%wave%kmy(i2)
          do i1=1,size(wfld,1)
            k_y =kmy+ sign*khy(i1,i3)*kmy
            k = sqrt(step%kx(i1,i3)**2 + k_y**2)
            itable =max(1, min(int(1 + k/ko / dsr%d) , dsr%n))
            wfld(i1,i2,i3,i4,i5)=wfld(i1,i2,i3,i4,i5)*step%ctable(itable)
          end do
        end do
      end do
    end do
        
#ifdef DEBUG
   write(0,*) "OUT OF WEI_WEM"
#endif
  end subroutine

  subroutine calc_khy(ko,inv_norm_ko,wave)
    integer :: is,ir,ihx,imx
    real :: kzs,kzr,ko,inv_norm_ko
    type(wspace) :: wave

!write(0,*) "1",ko,inv_norm_ko,size(khy,1),size(khy,2)
!write(0,*) "2",size(wave%ksx,1),size(wave%ksx,2),size(wave%krx,1),size(wave%krx,2)

    do ihx=1,size(khy,2)
      do imx=1,size(khy,1)
        is = min(int(1 + abs(wave%ksx(imx,ihx)*inv_norm_ko)) , nkxko)
        ir = min(int(1 + abs(wave%krx(imx,ihx)*inv_norm_ko)) , nkxko)

!write(0,*) ihx,imx,is,ir

        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
