module wei_wem_mod
  use wei_kmap_mod
  use wei_table_mod
  use wei_step_mod
  implicit none

!  type wem_local
!    complex,  pointer:: amp(:,:)
!    real, pointer:: temp1(:,:)
!    real, pointer:: temp2(:,:)
!  end type
      !  type(wem_local),private,save :: local(5)


  contains 

  subroutine wem_init(wave,fk)
    type(wspace)   :: wave  ! wavefield description
    type(fktype)   :: fk    ! fk structure

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

  end subroutine


  subroutine wei_wem(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
    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)
!     step%ctable(i)=exp(dsr%val(i)*ko*step%dstep)
     k=ko*step%dstep*dsr%phase(i)
     step%ctable(i)=dsr%amp(i)*cmplx(cos(k),sin(k))
   end do

    do i4=1,size(wfld,4)
      do i3=1,size(wfld,3)
        do i2=1,size(wfld,2)
          do i1=1,size(wfld,1)
            k = sqrt(step%kx(i1,i3)**2 + step%ky(i2,i4)**2)
            itable =max(1, min(int(1 + k/ko / dsr%d) , dsr%n))
!            ikz=    ko *dsr%val(itable) *step%dstep
!            temp1(i1)= real(ikz)
!            amp(i1)=cmplx(cos(aimag(ikz)),sin(aimag(ikz)))
!          end do
!#ifdef AMD
!          call vrsa_expf (size(temp2,1), temp1,temp2)
!#else
!          do i1=1,size(temp2,1)
!            temp2(i1)=exp(temp1(i1))
!          end do
!#endif
!          do i1=1,size(temp2,1)
            wfld(i1,i2,i3,i4,i5)=wfld(i1,i2,i3,i4,i5)*step%ctable(itable)
!if(i1==50) write(0,*) "CHECK THIS",itable,step%ctable(itable),temp2(i1),amp(i1)
!            wfld(i1,i2,i3,i4,i5)=wfld(i1,i2,i3,i4,i5)*temp2(i1)*amp(i1)
! write(0,*) "CHECK1THIS",i1,itable,step%ctable(itable),amp(i1),temp2(i1)
          end do
        end do
      end do
    end do
        
#ifdef DEBUG
   write(0,*) "OUT OF WEI_WEM"
#endif

  end subroutine

end module
