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

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

  real, dimension(:,:,:,:), pointer, private :: sxy
  real, private :: eds !! epsilon(ds)
contains
  !----------------------------------------------------------------
!!$
!!$=head1 NAME
!!$
!!$ weislo1_init - Initialize single slowness operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weislo1_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize single slowness operator
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !! one reference slowness
  integer function weislo1_init() result(st)
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weislo1_init'
    if(pro%debug) call in(unit)
#endif
    !! do nothing
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weislo1_init
!!$
!!$=head1 NAME
!!$
!!$ weisloN_init - Initialize multiple slowness operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weisloN_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize multiple slowness operator
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weislo1_init> L<weisloN>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------  
  !! one reference slowness
  integer function weislo1() result(st)
    integer :: izs
#ifdef DEBUG
    character(len=128), parameter :: unit='weislo1'
    if(pro%debug) call in(unit)
#endif

    !! spread the slowness in its PBA
    call weislow_push(S)

    !! s(z): FKop slowness
    do izs=1,pro%nzs
       slo%nrs(izs) = 1 !! one reference slowness
       slo%ssz(izs,1) = slowness_mean(S,izs)
       slo%srz(izs,1) = slowness_mean(S,izs)
    end do

    slo%si=1.

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weislo1
!!$
!!$=head1 NAME
!!$
!!$ weisloN - multiple slowness operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weisloN()
!!$
!!$=head1 DESCRIPTION
!!$
!!$  multiple slowness operator
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$

!!$
!!$=head1 NAME
!!$
!!$ weislo1 - single slowness operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weislo1()
!!$
!!$=head1 DESCRIPTION
!!$
!!$  single slowness operator
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weislo1_init> L<weisloN>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weisloN_init() result(st)
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weisloN_init'
    if(pro%debug) call in(unit)
#endif

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

    eds = epsilon(1.)

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weisloN_init
  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function weisloN() result(st)
    integer :: izs,ifk
    real    :: ds, so, sz
    real    :: sl,su
    real    :: smin,smax
#ifdef DEBUG
    character(len=128), parameter :: unit='weisloN'
    if(pro%debug) call in(unit)
#endif

    !! spread the slowness in its PBA    
    call weislow_push(S)

    !! min/max slowness in the interval
    do izs=1,pro%nzs
       slo%nrs(izs) = pro%nfk

       smin = slowness_min(S,izs)
       smax = slowness_max(S,izs)
       ds = (smax-smin) /pro%nfk

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

       !! interpolation map
       slo%si(:,:,:,:,izs)=pro%nfk

       sxy(:,:,:,:) &
       = .5* real(S%pba(:,:,:,:,izs)) &
       + .5*aimag(S%pba(:,:,:,:,izs))

       do ifk=1,pro%nfk                         !! [sl...sz...su]
          sz &
          = .5*slo%ssz(izs,ifk) &
          + .5*slo%srz(izs,ifk)                 !! s(z)
          sl = sz-ds/2-eds
          su = sz+ds/2
          where((sxy >= sl) .and. (sxy < su))
             slo%si(:,:,:,:,izs)=ifk
          end where
       end do
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weisloN
  !----------------------------------------------------------------
end module wei_slo



