!!.. shot prifile migration for nws frequencies and nzs depth steps
!!.. the output is the image model(nx, ny, nzs)
!!.. VLop = function dealing with velocity
!!.. EPop = extrapolator
!!.. IGin = imaging condition 
module wei_shotmig_kernel_mod
  use  wei_global_defi_mod, only: vref, vmod, model, wfld_d, wfld_u,   &
                                  kx2, ky2, pgrid, ogrid, lloyd, nref, &
                                  decon, eps, maxang, istep, astep, rstep, nstep
  use  wei_type_mod       , only: axis
  use  wei_taper_mod

  implicit none

  real   , dimension(:)      , pointer, private :: vref_slice
  real   , dimension(:,:)    , pointer, private :: vmod_slice
  real   , dimension(:,:,:,:), pointer, private :: model_slice
  complex, dimension(:,:)    , pointer, private :: wfld_slice_d
  complex, dimension(:,:)    , pointer, private :: wfld_slice_u
  integer                             , private :: nws, nzs
  real                                , private :: dws, dzs
  real                                , private :: ows, ozs


contains
!========================================================================
  integer function wei_shotmig_kernel_init(aws_in, azs_in, VLin, FXin, WIin, EPin, IGin) result(stat)
    type(axis)                      :: aws_in
    type(axis)                      :: azs_in
    optional                        :: WIin
    optional                        :: EPin
    interface
      
      integer function VLin(lloyd_in, nref_in) result(stat)
        integer                     :: nref_in
        logical                     :: lloyd_in
      end function

      integer function FXin(kx2_in, ky2_in, dz_in, maxang_in) result(stat)
        real, dimension(:), pointer :: kx2_in
        real, dimension(:), pointer :: ky2_in
        real                        :: dz_in
        real                        :: maxang_in
      end function

      integer function WIin(nref_in) result(stat)
        integer                     :: nref_in
      end function

      integer function EPin(kx2_in, ky2_in, dz_in, maxang_in, nref_in, FXin, WIin) result(stat)
        real, dimension(:), pointer :: kx2_in
        real, dimension(:), pointer :: ky2_in
        real                        :: dz_in
        real                        :: maxang_in
        integer                     :: nref_in

        interface
          integer function FXin(kx2_in, ky2_in, dz_in, maxang_in) result(stat)
            real, dimension(:), pointer :: kx2_in
            real, dimension(:), pointer :: ky2_in
            real                        :: dz_in
            real                        :: maxang_in
          end function
          integer function WIin(nref_in) result(stat)
            integer                     :: nref_in
          end function
        end interface
      end function

      integer function IGin(decon_in, eps_in, ogrid_in, pgrid_in) result(stat)
        use wei_type_mod
        logical                     :: decon_in
        real                        :: eps_in
        type(out_grid)              :: ogrid_in
        type(pro_grid)              :: pgrid_in
      end function

    end interface

    nws = aws_in%n; nzs = azs_in%n
    dws = aws_in%d; dzs = azs_in%d
    ows = aws_in%o; ozs = azs_in%o

    stat = VLin(lloyd, nref)
    stat = EPin(kx2, ky2, dzs, maxang, nref, FXin, WIin)
    stat = IGin(decon, eps, ogrid, pgrid)

  end function
!===========================================================================
  integer function wei_shotmig_kernel_oper(VLop, FXop, WIop, EPop, IGop) result(stat)
    optional                               :: WIop
    optional                               :: EPop 
    interface

      integer function VLop(vmod_slice, vref_slice) result(stat)
        real, dimension(:,:), pointer      :: vmod_slice
        real, dimension(:)  , pointer      :: vref_slice
      end function

      integer function FXop(adj, opsign, ws, vref, vmod_slice, wfld_slice) result(stat)
        logical                            :: adj
        integer                            :: opsign
        real                               :: ws
        real                               :: vref
        real   , dimension(:,:), pointer   :: vmod_slice
        complex, dimension(:,:), pointer   :: wfld_slice
      end function

      integer function WIop(adj, add, vref_slice, vmod_slice, wfld_slice, wfld_cube) result(stat)
        logical                            :: adj, add
        real   , dimension(:)    , pointer :: vref_slice
        real   , dimension(:,:)  , pointer :: vmod_slice
        complex, dimension(:,:)  , pointer :: wfld_slice
        complex, dimension(:,:,:), pointer :: wfld_cube
      end function
  
      integer function EPop(adj, opsign, ws, vref_slice, vmod_slice, wfld_slice, FXop, WIop) result(stat)
        logical                           :: adj
        integer                           :: opsign
        real                              :: ws
        real,    dimension(:)  , pointer  :: vref_slice
        real,    dimension(:,:), pointer  :: vmod_slice
        complex, dimension(:,:), pointer  :: wfld_slice
        interface
          integer function FXop(adj, opsign, ws, vref, vmod_slice, wfld_slice) result(stat)
            logical                            :: adj
            integer                            :: opsign
            real                               :: ws
            real                               :: vref
            real   , dimension(:,:), pointer   :: vmod_slice
            complex, dimension(:,:), pointer   :: wfld_slice
          end function
          integer function WIop(adj, add, vref_slice, vmod_slice, wfld_slice, wfld_cube) result(stat)
            logical                            :: adj, add
            real   , dimension(:)    , pointer :: vref_slice
            real   , dimension(:,:)  , pointer :: vmod_slice
            complex, dimension(:,:)  , pointer :: wfld_slice
            complex, dimension(:,:,:), pointer :: wfld_cube
          end function
        end interface
      end function

      integer function IGop(wfld_slice_d, wfld_slice_u, model_slice) result(stat)
        complex, dimension(:,:)    , pointer  :: wfld_slice_d
        complex, dimension(:,:)    , pointer  :: wfld_slice_u
        real   , dimension(:,:,:,:), pointer  :: model_slice
      end function

    end interface

    integer         :: izs, iws
    real            :: ws, rvref
    integer         :: percent


!    write(0,*)"nzs, nws, dws, ows, dz", nzs, nws, dws, ows, dzs

    do izs=1, nzs

       model_slice => model(:,:,:,:,izs)
       vmod_slice  => vmod(:,:,izs)
       vref_slice  => vref(:,izs)

       stat = VLop(vmod_slice, vref_slice)
       !rvref = sum(vmod_slice)/size(vmod_slice)

       do iws=1, nws
          ws = ows + (iws-1)*dws

          wfld_slice_d => wfld_d(:,:,iws)
          wfld_slice_u => wfld_u(:,:,iws)

          stat = EPop(.true., +1, ws, vref_slice, vmod_slice, wfld_slice_d, FXop, WIop)
          stat = EPop(.true., -1, ws, vref_slice, vmod_slice, wfld_slice_u, FXop, WIop)

          !stat = FXop(.true., +1, ws, rvref, vmod_slice, wfld_slice_d)
          !stat = FXop(.true., -1, ws, rvref, vmod_slice, wfld_slice_u)

          stat = IGop(wfld_slice_d, wfld_slice_u, model_slice)

          if(associated(wfld_slice_d)) nullify(wfld_slice_d)
          if(associated(wfld_slice_u)) nullify(wfld_slice_u)

          istep = istep + 1
 
          if (istep == astep) then
             percent = nint(real(astep)/real(nstep)*100.)
             write(0,*)"Finished:",percent,"%"
             astep = astep + rstep
          end if

       end do

       if(associated(model_slice)) nullify(model_slice)
       if(associated(vmod_slice))  nullify(vmod_slice)
       if(associated(vref_slice))  nullify(vref_slice)

    end do

    stat = 0

  end function
!====================================================================
  integer function wei_shotmig_kernel_clos(VLnd, FXnd, WInd, EPnd, IGnd) result(stat)
    optional        :: WInd
    optional        :: EPnd
    interface 
      integer function VLnd() result(stat)
      end function
      integer function FXnd() result(stat)
      end function
      integer function WInd() result(stat)
      end function
      integer function EPnd(FXnd, WInd) result(stat)
        interface
          integer function FXnd() result(stat)
          end function
          integer function WInd() result(stat)
          end function
        end interface
      end function  
      integer function IGnd() result(stat)
      end function
    end interface

    stat = VLnd()
    stat = EPnd(FXnd, WInd)
    stat = IGnd()

  end function

end module








