module wei_kxmig
  !! fk-fx migration
  !! receive wfld at (izs,iws,ith)
  !! return  wfld downward/upward continued
  !! needs FKop, FXop

  use wei_util
  use wei_process
  use wei_obj_S
  use wei_fft
  use wei_kmap
  !---------------
  implicit none

  !! temporary interpolation arrays
  complex, dimension(:,:,:,:,:), pointer, private :: tmp1,tmp2

contains
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine weikxmig(wfld,iws,izs,ith,KXMIGop,FKop,FXop)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer,external,optional ::    FKop,FXop
    integer :: izs,iws,ith
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weikxmig'
#endif
    interface !!-------------------------
       integer function KXMIGop(wfld,iws,izs,ith,FKop,FXop) result(st)
         complex, dimension(:,:,:,:,:), pointer :: wfld
         integer, intent(in)       :: iws,izs,ith
         integer,external,optional :: FKop,FXop
       end function KXMIGop
    end interface !!---------------------
#ifdef DEBUG
    if(pro%depar) call in(unit)
#endif

    st = KXMIGop(wfld,iws,izs,ith,FKop,FXop)

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
  end subroutine weikxmig
  !----------------------------------------------------------------
!!$
!!$=head1 NAME
!!$
!!$ weikxmig1_init - Single velocity KX-migration initialization
!!$
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr= weikxmig1_init()
!!$
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize a single-reference kx-migration velocity operation
!!$
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weikxmig1> L<weikxmigN_init>
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
        integer function weikxmig1_init() result(st)
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weikxmig1_init'
    if(pro%debug) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weikxmig1_init
!!$
!!$=head1 NAME
!!$
!!$ weikxmig1 - Single velocity KX-migration 
!!$
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr= weikxmig1(wfld,iws,izs,ith,FKOP,FXOP)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item wfld complex(:,:,:,:,:)
!!$
!!$  Wavefield
!!$
!!$=item iws - integer
!!$
!!$ Current frequency slice
!!$
!!$=item izs
!!$
!!$ Current depth slice
!!$
!!$=item ith
!!$
!!$ Current thread
!!$
!!$=item FKOP FK interface L<weiwem>
!!$
!!$=item FXOP FX interface L<weissf>
!!$
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Perform  a single-reference kx-migration velocity operation
!!$
!!$          st   = FXop(iws,izs,ifk,ith,wfld) !! f-x op
!!$
!!$          call weiift(iws,        ith,wfld)
!!$
!!$          st   = FKop(iws,izs,ifk,ith,wfld) !! f-k op
!!$
!!$          call weifft(iws,        ith,wfld)
!!$
!!$
!!$=head1 SEE ALSO
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weikxmigN> L<weikxmig1_init>
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  !! 1 reference slowness
  integer function weikxmig1(wfld,iws,izs,ith,FKop,FXop) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in)       :: iws,izs,ith
    integer                   :: ifk
#ifdef DEBUG
    character(len=128), parameter :: unit='weikxmig1'
#endif
    interface !!----------------------------------------------
       integer function FKop(iws,izs,ifk,ith,wfld)      result(st)
         integer, intent(in) :: iws,izs,ifk,ith
         complex, dimension(:,:,:,:,:), pointer :: wfld
       end function FKop
       integer function FXop(iws,izs,ifk,ith,wfld) result(st)
         integer, intent(in) :: iws,izs,ifk,ith
         complex, dimension(:,:,:,:,:), pointer :: wfld
       end function FXop
    end interface !!------------------------------------------
#ifdef DEBUG
    if(pro%depar) call in(unit)
#endif

    call redo_kmaps(pro%ihxb_beg(izs),pro%ihxb_end(izs), &
       pro%ihyb_beg(izs),pro%ihyb_end(izs))

    if(pro%adj) then !! MIGRATION
       st   = FXop(iws,izs,  1,ith,wfld) !! f-x op
       call weiift(iws,izs,    ith,wfld)
       st   = FKop(iws,izs,  1,ith,wfld) !! f-k op
       call weifft(iws,izs,    ith,wfld)
    else             !! MODELLING
       call weiift(iws,izs,    ith,wfld)
       st   = FKop(iws,izs,  1,ith,wfld) !! f-k op
       call weifft(iws,izs,    ith,wfld)
       st   = FXop(iws,izs,  1,ith,wfld) !! f-x op
    end if

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
    st=OK
  end function weikxmig1
!!$
!!$=head1 NAME
!!$
!!$ weikxmigN_init - Multi-velocity KX-migration initialization
!!$
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr= weikxmigN_init()
!!$
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize a multi-reference kx-migration velocity operation
!!$
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weikxmigN> L<weikxmig1_init>
!!$
!!$
!!$=head1 sepweif90
!!$
!!$ sepweif90
!!$
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weikxmigN_init() result(st)
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weikxmigN_init'
    if(pro%debug) call in(unit)
#endif

    !! allocate tmp arrays for interpolation of multiple fk ops
    allocate(tmp1(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nth),stat=er)
    if(er/=0) call seperr('cant allocate tmp1')
    allocate(tmp2(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nth),stat=er)
    if(er/=0) call seperr('cant allocate tmp2')

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weikxmigN_init
!!$
!!$=head1 NAME
!!$
!!$ weikxmigN - Multi- velocity KX-migration 
!!$
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr= weikxmigN(wfld,iws,izs,ith,FKOP,FXOP)
!!$
!!$=head1 PARAMETERS
!!$
!!$=over 4
!!$
!!$=item wfld complex(:,:,:,:,:)
!!$
!!$  Wavefield
!!$
!!$=item iws - integer
!!$
!!$ Current frequency slice
!!$
!!$=item izs
!!$
!!$ Current depth slice
!!$
!!$=item ith
!!$
!!$ Current thread
!!$
!!$=item FKOP FK interface L<weiwem>
!!$
!!$=item FXOP FX interface L<weissf>
!!$
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Perform  a single-reference kx-migration velocity operation
!!$
!!$          st   = FXop(iws,izs,ifk,ith,wfld) !! f-x op
!!$
!!$          call weiift(iws,izs, izs,ith,wfld)
!!$
!!$          st   = FKop(iws,izs,ifk,ith,wfld) !! f-k op
!!$
!!$          call weifft(iws,izs,ith,wfld)
!!$
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weikxmig1> L<weikxmigN_init>
!!$
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weikxmigN(wfld,iws,izs,ith,FKop,FXop) result(st)
    complex, dimension(:,:,:,:,:), pointer :: wfld
    integer, intent(in)       :: iws,izs,ith
    integer                   :: ifk
#ifdef DEBUG
    character(len=128), parameter :: unit='weikxmigN'
#endif
    interface !!----------------------------------------------
       integer function FKop(iws,izs,ifk,ith,wfld) result(st)
         integer, intent(in) :: iws,izs,ifk,ith
         complex, dimension(:,:,:,:,:), pointer :: wfld
       end function FKop
       integer function FXop(iws,izs,ifk,ith,wfld) result(st)
         integer, intent(in) :: iws,izs,ifk,ith
         complex, dimension(:,:,:,:,:), pointer :: wfld
       end function FXop
    end interface !!------------------------------------------
#ifdef DEBUG
    if(pro%depar) call in(unit)
#endif
    if(pro%adj) then !! MIGRATION
       !----------------------------------------------------------------       



       tmp2(:,:,:,:,ith) = wfld(:,:,:,:,iws)
       tmp1(:,:,:,:,ith) = 0.

        call redo_kmaps(pro%ihxb_beg(izs),pro%ihxb_end(izs), &
       pro%ihyb_beg(izs),pro%ihyb_end(izs))


   do ifk=1,slo%nrs(izs)


          where(slo%si(:,:,:,:,izs) == ifk)
             wfld(:,:,:,:,iws) = tmp2(:,:,:,:,ith)
          elsewhere
             wfld(:,:,:,:,iws) = 0.
          end where

          st   = FXop(iws,izs,ifk,ith,wfld) !! f-x op
          call weiift(iws,izs,    ith,wfld)
          st   = FKop(iws,izs,ifk,ith,wfld) !! f-k op
          call weifft(iws,izs,    ith,wfld)

          tmp1(:,:,:,:,ith) = &
          tmp1(:,:,:,:,ith) + &
          wfld(:,:,:,:,iws)

       end do
       wfld(:,:,:,:,iws) = tmp1(:,:,:,:,ith)
       !----------------------------------------------------------------
    else             !! MODELLING
       !----------------------------------------------------------------
       tmp1(:,:,:,:,ith) = wfld(:,:,:,:,iws)
       tmp2(:,:,:,:,ith) = 0.
   do ifk=1,slo%nrs(izs)
          wfld(:,:,:,:,iws) = &
          tmp1(:,:,:,:,ith)

          call weiift(iws,izs,    ith,wfld)
          st   = FKop(iws,izs,ifk,ith,wfld) !! f-k op
          call weifft(iws,izs,    ith,wfld)
          st   = FXop(iws,izs,ifk,ith,wfld) !! f-x op

          where(slo%si(:,:,:,:,izs) == ifk)
             tmp2(:,:,:,:,ith) = &
             tmp2(:,:,:,:,ith) + &
             wfld(:,:,:,:,iws)
          end where

       end do
       wfld(:,:,:,:,iws) = tmp2(:,:,:,:,ith)
       !----------------------------------------------------------------
    end if

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
    st=OK
  end function weikxmigN
  !----------------------------------------------------------------  
end module wei_kxmig


!! multiple fk ops
!!
!!       D(z)
!!       |
!!       v
!!       T1
!!
!!    D1 D2 D3 ...
!!
!!    S1 S2 S3 ...
!!     \ | /
!!       T2
!!       |
!!       v
!!       D(z+dz)



  
