module weiop_sll
  !! WEI slowness operator
  !! Lloyd algorithm

  use wei_util
  use wei_pro
  use wei_par
  use wei_S
  use lloyd_mod
  use wei_dump
  !---------------
  implicit none

  real, dimension(:,:,:,:), pointer, private :: sxy
  real, dimension(:),       pointer, private :: bound,ref_slow,s_slice
  integer,                           private :: nsavg
contains
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weiop_sllN_init() result(st)
    integer :: er
    integer :: mpierr
    character(len=128), parameter :: unit='weiop_sllN_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,unit)
    allocate(sxy(     wpro%amx%n,wpro%amy%n,wpro%ahx%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate sxy')

    call weimem_add(4,wpro%nfk,unit)
    allocate(bound(   wpro%nfk))
    call weimem_add(4,wpro%nfk,unit)
    allocate(ref_slow(wpro%nfk))
    call weimem_add(4,wpro%amx%n*wpro%amy%n,unit)
    allocate(s_slice( wpro%amx%n*wpro%amy%n))

    !! INIT Lloyd algorithm
    call init_lloyd()

    if(wpar%impi==0) call from_param("nsavg",nsavg,0)
    if(wpar%nmpi >0) call MPI_BCAST(nsavg, 1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_sllN_init
  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function weiop_sllN(X) result(st)
    type(weslow), intent(inout) :: X
    integer :: izs,i4,i3,i2,i1
    integer :: n1,n2,n3,n4
    character(len=128), parameter :: unit='weiop_sllN'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    n1=size(sxy,1)
    n2=size(sxy,2)
    n3=size(sxy,3)
    n4=size(sxy,4)

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

    izsloop: do izs=1,wpro%nzs
       s_slice=reshape(X%bin(:,:,izs),(/size(X%bin(:,:,izs))/))

       slo%nr_s(izs)=lloyd_go(s_slice,ref_slow,wpro%nfk,bound,wpar%impi==0)
       slo%nr_r(izs)=slo%nr_s(izs)

       slo%rs_s(izs,1:size(ref_slow)) = ref_slow
       slo%rs_r(izs,1:size(ref_slow)) = ref_slow

       !! source interpolation map
       sxy(:,:,:,:) = slo%ss_s(:,:,:,:,izs)
       do  i4=1,size(sxy,4)
          do  i3=1,size(sxy,3)
             do  i2=1,size(sxy,2)
                do  i1=1,size(sxy,1)
                   slo%si_s(i1,i2,i3,i4,izs)=1
                   do while(sxy(i1,i2,i3,i4) > bound(slo%si_s(i1,i2,i3,i4,izs)))
                      slo%si_s(i1,i2,i3,i4,izs)=&
                      slo%si_s(i1,i2,i3,i4,izs)+1
                   end do
                end do
             end do
          end do
       end do
       slo%sw_s(:,:,:,:,izs)=slo%si_s(:,:,:,:,izs)
       call weiavg(slo%sw_s,izs)

       !! receiver interpolation map
       sxy(:,:,:,:) = slo%ss_r(:,:,:,:,izs)
       do  i4=1,size(sxy,4)
          do  i3=1,size(sxy,3)
             do  i2=1,size(sxy,2)
                do  i1=1,size(sxy,1)
                   slo%si_r(i1,i2,i3,i4,izs)=1
                   do while(sxy(i1,i2,i3,i4) > bound(slo%si_r(i1,i2,i3,i4,izs)))
                      slo%si_r(i1,i2,i3,i4,izs)=&
                      slo%si_r(i1,i2,i3,i4,izs)+1
                   end do
                end do
             end do
          end do
       end do
       slo%sw_r(:,:,:,:,izs)=slo%si_r(:,:,:,:,izs)
       call weiavg(slo%sw_r,izs)

#ifdef DEBUG
       if(wpar%impi==0 .and. wpro%dbg>2) write(0,*) izs,slo%nr_s(izs),1./slo%rs_s(izs,1:slo%nr_s(izs))
       if(wpar%impi==0 .and. wpro%dbg>2) write(0,*) izs,slo%nr_r(izs),1./slo%rs_r(izs,1:slo%nr_r(izs))
#endif

    end do izsloop

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
    st=weiOK
  end function weiop_sllN
  !----------------------------------------------------------------
  subroutine weiavg(array,izs)
    real, dimension(:,:,:,:,:), pointer :: array
    integer,                 intent(in) :: izs
    integer :: n4,n3,n2,n1
    integer :: i4,i3,i2,i1
    integer :: k4,k3,k2,k1
    integer :: l4,l3,l2,l1
    integer :: ns,ne

    ns=nsavg
    ne=1

    n4=size(sxy,4)
    n3=size(sxy,3)
    n2=size(sxy,2)
    n1=size(sxy,1)

    sxy(:,:,:,:) = array(:,:,:,:,izs)

    do  i4=1,n4
       k4=i4-ns; k4=max(k4, 1)
       l4=i4+ns; l4=min(l4,n4)
       do  i3=1,n3
          k3=i3-ns; k3=max(k3, 1)
          l3=i3+ns; l3=min(l3,n3)
          do  i2=1,n2
             k2=i2-ns; k2=max(k2, 1)
             l2=i2+ns; l2=min(l2,n2)
             do  i1=1,n1
                k1=i1-ns; k1=max(k1, 1)
                l1=i1+ns; l1=min(l1,n1)
                ne=(l1-k1+1)*(l2-k2+1)*(l3-k3+1)*(l4-k4+1)

                array(i1,i2,i3,i4,izs) = sum( sxy(k1:l1,k2:l2,k3:l3,k4:l4)) / ne
             end do
          end do
       end do
    end do

  end subroutine weiavg
!----------------------------------------------------------------
end module weiop_sll












