module weiop_slo
  !! WEI slowness operator
  !! for FK-FX downward continuation

  use wei_util
  use wei_pro
  use wei_S
  use wei_attr
  use wei_dump
  !---------------
  implicit none

  real, dimension(:,:,:,:), pointer, private :: sxy
  real, private :: eds !! epsilon(ds)
contains
  !----------------------------------------------------------------
  !! one reference slowness
  integer function weiop_slo1_init() result(st)
    integer :: er
    character(len=128), parameter :: unit='weiop_slo_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_slo1_init
  !----------------------------------------------------------------  
  !! one reference slowness
  integer function weiop_slo1(X) result(st)
    type(weslow), intent(inout) :: X
    integer :: izs
    character(len=128), parameter :: unit='weiop_slo1'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    !! spread the slowness in its PBA
    call weiS_push(X)

    !! s(z): FKop slowness
    do izs=1,wpro%nzs
       slo%nr_s(izs) = 1 !! one reference slowness
       slo%nr_r(izs) = 1 !! one reference slowness

       slo%rs_s(izs,1) = weiS_mean(X,izs)
       slo%rs_r(izs,1) = weiS_mean(X,izs)
    end do

    slo%si_s=1.
    slo%si_r=1.

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_slo1
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weiop_sloN_init() result(st)
    integer :: er
    character(len=128), parameter :: unit='weiop_sloN_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    call weimem_add(4,wpro%amx%n*wpro%amy%n*wpro%ahx%n*wpro%ahy%n)
    allocate(sxy(     wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate sxy')

    eds = epsilon(1.)

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_sloN_init
  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function weiop_sloN(X) result(st)
    type(weslow), intent(inout) :: X
    integer :: izs,ifk
    real    :: ds, so, sz
    real    :: sl,su
    real    :: smin,smax
    character(len=128), parameter :: unit='weiop_sloN'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    !! spread the slowness in its PBAs    
    call weiS_push(X)

    !! min/max slowness in the interval
    izsloop: do izs=1,wpro%nzs
       slo%nr_s(izs) = wpro%nfk
       slo%nr_r(izs) = wpro%nfk

       smin = weiS_min(X,izs)
       smax = weiS_max(X,izs)
       ds = (smax-smin) / wpro%nfk

       !! s(z): FKop slowness
       do ifk=1,wpro%nfk
          slo%rs_s(izs,ifk) = smin + (ifk-.5)*ds
          slo%rs_r(izs,ifk) = smin + (ifk-.5)*ds
       end do

       !! interpolation map
       slo%si_s(:,:,:,:,izs)=wpro%nfk
       slo%si_r(:,:,:,:,izs)=wpro%nfk

       !! source
       sxy(:,:,:,:) = slo%ss_s(:,:,:,:,izs)

       do ifk=1,wpro%nfk                        !! [sl...sz...su]
          sz = slo%rs_s(izs,ifk)
          sl = sz-ds/2-eds
          su = sz+ds/2
          where((sxy >= sl) .and. (sxy < su))
             slo%si_s(:,:,:,:,izs)=ifk
          end where
       end do
       slo%sw_s(:,:,:,:,izs)=slo%si_s(:,:,:,:,izs)

       !! receiver
       sxy(:,:,:,:) = slo%ss_r(:,:,:,:,izs)

       do ifk=1,wpro%nfk                        !! [sl...sz...su]
          sz = slo%rs_r(izs,ifk)
          sl = sz-ds/2-eds
          su = sz+ds/2
          where((sxy >= sl) .and. (sxy < su))
             slo%si_r(:,:,:,:,izs)=ifk
          end where
       end do
       slo%sw_r(:,:,:,:,izs)=slo%si_r(:,:,:,:,izs)

    end do izsloop

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_sloN
  !----------------------------------------------------------------
end module weiop_slo



