module weiop_slx
  !! WEI slowness operator
  !! for FK-FX downward continuation
  !! w/ salt logic
  use wei_util
  use wei_process
  use wei_obj_s
  use wei_qksort
  use wei_hist
  use wei_dump
  !---------------
  implicit none

  real, dimension(:,:,:,:), pointer, private :: sxy
  real, dimension(:),       pointer, private :: s1d

  integer, private :: n1d    !! size of the 1D-packed slowness slice

  !! type of reference slowness
  character(len=3), private :: sloref
  !! for equally spaced intervals, the reference slowness can be : 
  !! - avg - average value
  !! - med -  median value
  !! - min -     min value
  !! - max -     max value 
  !! - mid -  middle of the interval
  !! - his - most popular value in the interval (based on a histogram)

  real,    private :: vsalt, ssalt
  real,    private :: minds

  real,    private :: slgamax
contains
  !----------------------------------------------------------------
  integer function weiop_slxN_init() result(st)
    integer :: er
    integer :: mpierr
#ifdef DEBUG
    character(len=128), parameter :: unit='weiop_slxN_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')

    n1d=pro%amx%n*pro%amy%n
    allocate(s1d(n1d),stat=er)
    if(er/=0) call seperr('cant allocate s1d')

    !! salt velocity
    if(pro%mpi_ithread==0) call from_param("vsalt",vsalt,0.)
    if(pro%mpi_nthread >0) call MPI_BCAST(vsalt, 1,MPI_REAL,0,MPI_COMM_WORLD,mpierr)
    ssalt=0.
    if(vsalt >0.) ssalt=1./vsalt
    if(pro%mpi_ithread==0 .and. pro%verb) call weimessage("vsalt",vsalt)
    if(pro%mpi_ithread==0 .and. pro%verb) call weimessage("ssalt",ssalt)

    !! min slowness variation
    !! =>   nfkeff = 1+ (smax-smin)/minds 
    if(pro%mpi_ithread==0) call from_param("minds",minds,0.0)
    if(pro%mpi_nthread >0) call MPI_BCAST(minds, 1,MPI_REAL,0,MPI_COMM_WORLD,mpierr)
    if(pro%mpi_ithread==0 .and. pro%debug) call weimessage("minds",minds)

    !! type of reference slowness
    if(pro%mpi_ithread==0) call from_param("sloref",sloref,'mid')
    if(pro%mpi_nthread >0) call MPI_BCAST(sloref, 3,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr)
    if(pro%mpi_ithread==0 .and. pro%verb) call weimessage("sloref",sloref)

    call weiqksort_init(n1d)
    if(sloref=='his') then
       call weihist_init()
    end if

    !! type 6 needs the max dip angle propagating correctly
    if(sloref=='geo') then
       if(pro%mpi_ithread==0) call from_param("slgamax",slgamax,60.) !! degrees
       if(pro%mpi_nthread >0) call MPI_BCAST(slgamax, 1,MPI_REAL,0,MPI_COMM_WORLD,mpierr)
       if(pro%mpi_ithread==0 .and. pro%debug) call weimessage("slgamax",slgamax)
    end if

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=0
  end function weiop_slxN_init
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weiop_slxN(X) result(st)
    type(weslow), intent(inout) :: X
    integer :: izs
    real    :: smin,smax, saltslow
    integer :: snfk !! nfk in sediments
#ifdef DEBUG
    character(len=128), parameter :: unit='weiop_slxN'
    if(pro%debug) call in(unit)
#endif

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

          izsloop: do izs=1,pro%nzs
       slo%si(:,:,:,:,izs)=1  !! by default, select wfld for nfk=1

       !! slowness(mx,my,hx,hy) @ izs
       sxy(:,:,:,:) &
       = .5* real(X%pba(:,:,:,:,izs)) &
       + .5*aimag(X%pba(:,:,:,:,izs))

       !! min and max slowness @ izs
       smin = slowness_min(X,izs)
       smax = slowness_max(X,izs)

       if(identify_salt(izs,smin)) then !! slowness w/  SALT
          !! find average salt slowness
          saltslow = find_saltslow(X,izs)
          !! setup ifk=1 w/ for the salt
          slo%ssz(izs,1) = saltslow
          slo%srz(izs,1) = saltslow
          where(sxy <= saltslow)
             slo%si(:,:,:,:,izs)=1
          end where
          !! end salt -------------------

          !! restrict smin to the sediments
          smin = find_nosalt_smin(X,izs)

          snfk = pro%nfk-1 !! no of f-k velocities in sediments
          slo%nrs(izs) = 1
       else                             !! slowness w/o SALT
          snfk = pro%nfk   !! no of f-k velocities in sediments
          slo%nrs(izs) = 0
       end if

       if(smax>smin) then

          if(sloref/='geo') then
             call wei_setfkslo(X,izs,smin,smax,snfk)
          else
             call wei_geomsslo(  izs,smin,smax,snfk)
          end if

       else

          slo%nrs(izs) = slo%nrs(izs) + 1
          slo%ssz(izs,slo%nrs(izs)) = smin
          slo%srz(izs,slo%nrs(izs)) = smin
          where(sxy >= smin)
             slo%si(:,:,:,:,izs)=slo%nrs(izs)
          end where
       end if

    end do izsloop

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  st=0
  end function weiop_slxN
  !----------------------------------------------------------------
  logical function identify_salt(izs,smin) result(findsalt)
    integer, intent(in) :: izs
    real,    intent(in) :: smin
#ifdef DEBUG
    character(len=128), parameter :: unit='identfy salt'
    if(pro%debug) call in(unit)
#endif

    !! do I want to apply the salt logic?
    !! is there salt in the section?
    !! make sure nfk allows salt logic
    if(vsalt>0..and. smin<=ssalt .and. pro%nfk>1) then      
       findsalt=.true.
    else
       findsalt=.false.
    end if

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end function identify_salt
  !----------------------------------------------------------------
  !! find the true slowness in the salt region
  real function find_saltslow(X,izs) result(saltslow)
    type(weslow), intent(inout) :: X
    integer, intent(in) :: izs
    real    :: ss, slocal
    integer :: ks
    integer :: imx, imy
#ifdef DEBUG
    character(len=128), parameter :: unit='find_saltslow'
    if(pro%debug) call in(unit)
#endif

    ss=0.
    ks=0.
    do imy=1,pro%amy%n
       do imx=1,pro%amx%n
          slocal = real(X%s_s%bin(imx,imy,izs))

          if(slocal <= ssalt) then
             ks = ks + 1
             ss = ss + slocal
          end if
       end do
    end do
    if(ks/=0) then
       ss = ss/ks 
    else
       ss = ssalt
    end if

    saltslow = ss

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end function find_saltslow
  !----------------------------------------------------------------  
  !! find the min slowness outside the salt region
  real function find_nosalt_smin(X,izs) result(smin)
    type(weslow), intent(inout) :: X
    integer, intent(in) :: izs
    real    :: ss, slocal
    integer :: ks
    integer :: imx, imy
#ifdef DEBUG
    character(len=128), parameter :: unit='find_nosalt_smin'
    if(pro%debug) call in(unit)
#endif

    ss=1.e+6 !! init w/ a huge value - can't be salt

    do imy=1,pro%amy%n
       do imx=1,pro%amx%n
          slocal = real(X%s_s%bin(imx,imy,izs))

          !! am I in the sediment region
          if(   slocal > ssalt) then
             if(slocal < ss   ) ss=slocal
          end if

       end do
    end do

    smin = ss
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end function find_nosalt_smin
  !----------------------------------------------------------------  
  subroutine wei_setfkslo(X,izs,smin,smax,snfk)
    type(weslow), intent(inout) :: X
    integer, intent(in) :: izs        !! depth level
    real,    intent(in) :: smin,smax  !! min/max slowness
    integer, intent(in) :: snfk       !! max nfk outside salt

    integer :: ioff, nfkeff, ifk
    real    :: ds,    dseff

    integer :: i1d

    real    :: ssl,ssh
    integer :: isl,ish
    real    :: ssum
    integer :: snum

#ifdef DEBUG
    character(len=128), parameter :: unit='wei_setfkslo'
    if(pro%debug) call in(unit)
#endif

    !! offset the index if salt was detected
    !! ioff=(pro%nfk-snfk) (0=no salt) (1=salt)
    ioff = pro%nfk-snfk

    !! adjust the effective nfk w.r.t. the min ds (minds) value
    ds = (smax-smin)/snfk !! true ds
    if(minds<ds) then
       dseff  = ds          !! effective ds
       nfkeff = snfk        !! effective nfk
    else
       dseff  = minds
       nfkeff = floor((smax-smin)/dseff); nfkeff=max(nfkeff,1)
    end if
    !    if(pro%mpi_ithread==0) write(0,*) izs,minds,ds,dseff,nfkeff

    !! set the # of fk ops @ izs
    slo%nrs(izs) = slo%nrs(izs) + nfkeff

    !! pack the slowness slice in a 1-D array
    s1d = pack(real(X%s_s%bin(:,:,izs)),.true.)
    !! quick sort the 1-D array
    call weiqksort(s1d,1,n1d)

    !! set the values of the ref slowness
    do ifk=ioff+1,ioff+nfkeff
       ssl = smin + (ifk-ioff-1.)*dseff !! lowest  slowness
       ssh = smin + (ifk-ioff   )*dseff !! highest slowness

       isl=1                            !! index for ssl
       do i1d=1,n1d
          if(s1d(i1d)<ssl) isl=isl+1
       end do
       ish=n1d                          !! index for ssh
       do i1d=n1d,1,-1
          if(s1d(i1d)>ssh) ish=ish-1
       end do

       select case(sloref)
       case('min')
          slo%ssz(izs,ifk) =    ssl
       case('mid')
          slo%ssz(izs,ifk) = .5*ssl+.5*ssh
       case('max')
          slo%ssz(izs,ifk) =           ssh
       case('avg')
          slo%ssz(izs,ifk) = sum(s1d(isl:ish))/(ish-isl+1)
       case('med')
          slo%ssz(izs,ifk) = s1d((isl+ish)/2) 
       case('his')
          slo%ssz(izs,ifk) = weihist(izs,s1d,isl,ish)
       end select

       slo%srz(izs,ifk) = slo%ssz(izs,ifk)
    end do

    !! interpolation map
    !    if(pro%mpi_ithread==0) write(0,*) izs,"nfkeff=",nfkeff,"ioff=",ioff,"dseff=",dseff
    where(sxy >= smin)
       slo%si(:,:,:,:,izs) = ioff+min(1+nint((sxy-smin)/dseff),nfkeff)
    end where

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine wei_setfkslo
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine wei_geomsslo(izs,smin,smax,snfk)
    integer, intent(in) :: izs
    real,    intent(in) :: smin,smax
    integer, intent(in) :: snfk
    real                :: sinamax, sinbmax

    integer             :: nfktot, nfkeff
    integer             :: ifk, ioff
#ifdef DEBUG
    character(len=128), parameter :: unit='wei_geomsslo'
    if(pro%debug) call in(unit)
#endif

    ioff = pro%nfk-snfk


    !! sin of the max angle propagated correctly
    sinamax = sin( pro%pi * slgamax / 180.)

    !! # of f-k operators required by the smin-smax spread
    nfktot = 1+ floor(log(smin/smax) / log(sinamax));

    !! # nfkeff can't exceed the sediment nfk 
    nfkeff = min(nfktot,snfk)

    !! sin of the effective max angle propagated correctly
    sinbmax = exp( log(smin/smax)/nfkeff )

    ioff = pro%nfk-snfk

    slo%ssz(izs,ioff+1) = smin
    slo%srz(izs,ioff+1) = slo%ssz(izs,1)
    do ifk=ioff+2,ioff+nfkeff
       slo%ssz(izs,ifk) = slo%ssz(izs,ifk-1) /sinbmax  
       slo%srz(izs,ifk) = slo%ssz(izs,ifk)
    end do

    slo%nrs(izs) = slo%nrs(izs) + nfkeff

    do ifk=ioff+1,ioff+nfkeff
       where( sxy >= slo%ssz(izs,ifk)) slo%si(:,:,:,:,izs)=ifk      
    end do

    slo%ssz(izs,:) = slo%ssz(izs,:) /sinbmax
    slo%srz(izs,:) = slo%srz(izs,:) /sinbmax
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine wei_geomsslo
  !----------------------------------------------------------------
end module weiop_slx










