! get_wbav -- computes averaged weight of m.c particles (wbav) 
! for each species: if wbav (:,isb)=0.0 -- leave as is
!                   if wbav (1:k,isb)=0.0 -- wbav (1:k,isb)=wbav (k,isb)
!                   if wbav (mj-k:mj,isb)=0.0 -- wbav (mj-k:mj,isb)=wbav (mj-k-1,isb)
! fill any interior gaps, sp wbav(j,isb)>0.0 everywhere if any particles exist for species isb
! using linear interpolation
subroutine get_wbav
!
  use nbi_com
  implicit none
  integer :: j
  real(fp) :: swnmc_f  ! weithed sum sum over all zone [ wnmc(k,lsbeam) * fdep_k ]
  real(fp) :: fdep_k      ! ptcl weight profil, fdep_k = [ 1 - xi(k)] / wghta + x(k)
  real(fp) :: fdep_oa
!
!
! wbav_kill = wnmc/nmc
  call get_noholes_array(wnmc(1:mj,lsbeam),nmc(1:mj,lsbeam),xi(1:mj,2),mj,lsbeam,wbav_kill(1:mj,lsbeam))
!
!
!     estimate wbav if there is some power
!     if power is = 0.0 then wbav will be the same as one for the last time step with power >0.0
!     
!     wbav  will be zeroed there and new values will be set
            
  if(pinjs(lsbeam).gt.0.0_fp) then
     wbav(:,lsbeam) = 0._fp    
     
     if(wghtot.eq.0.0_fp) return 
     !
     
     swnmc_f = 0._fp
     do j=1,mj
        swnmc_f = swnmc_f + fdep_oa(sdep_oa(1:mj,lsbeam),mj,j)* fdep_k(xi(j,2),xwght0(lsbeam))
     enddo
     
     if(swnmc_f.gt.0.0_fp) then
        do j=1,mj
           wbav(j,lsbeam) =wghtot* fdep_k(xi(j,2),xwght0(lsbeam)) / swnmc_f
        enddo
     endif
  endif
end subroutine get_wbav


real(fp) function fdep_oa(sdepoa,mj,j_in)
  !calculates probability of deposition to orbit averaged zone <j_in>
  use iso_c_binding, only: fp => c_double
  implicit none
  real(fp), dimension(mj) :: sdepoa
  integer, intent(in) ::j_in,mj
  real(fp) :: sum_sdepoa
  integer :: isz

  sum_sdepoa= sum(sdepoa(1:mj))
  if(sum_sdepoa.eq.0.0_fp) then
     fdep_oa=0.0_fp
  else
     fdep_oa=sdepoa(j_in)/sum_sdepoa
  endif
  
end function fdep_oa


real(fp) function fdep_k(xi_j,wghta_1)
  use iso_c_binding, only: fp => c_double
  implicit none
  real(fp) :: xi_j,wghta_1 !wghta_1=1/wghta

  if(wghta_1.le.0.0_fp) then 
     ! not suppose to go here
     fdep_k = 0.0_fp
     return
  endif
  fdep_k = (1._fp - xi_j)*wghta_1 + xi_j
  
end function fdep_k



subroutine get_noholes_array(wnmc,nmc,xizc,mj,lsbeam,wbav)
!
!  patch "holes" (zero values) in profile,
!
!  done for mc radial profiles that are formed by summing of the
!  profile itself and its normalization
!
!    zprof(j) = [mc sum]f*dw
!    zwsum(j) = [mc sum]dw
!
!  and at normalization time
!
!    if (zwsum(j).gt.0.0) then
!         zprof(j)=zprof(j)/zwsum(j)
!    else
!         zprof(j)=0.0
!    endif
!
!  which can leave "holes" in the profile in case of lousy statistics.
!
!  *** so ***
!  fill these holes with flat extrapolation & linear interpolation
  use iso_c_binding, only: fp => c_double
  implicit none
  real(fp), intent(in), dimension(mj):: wnmc,nmc
  real(fp), intent(in), dimension(mj) ::xizc
  real(fp), intent(out), dimension(mj):: wbav
  integer, intent(in) :: mj,lsbeam
  integer i,k,j,i1,i2,iw,iw_right,iw_left
  integer, dimension(:), allocatable ::i_wbav0
  allocate(i_wbav0(mj)) ! storage for index where wbav=0.0
  i_wbav0=0
!  
  k=0 !# with value.eq.0
  do j=1,mj
     
     if(nmc(j).gt.0.0_fp) then
        wbav(j)=wnmc(j)/nmc(j)
     else
        k=k+1
        wbav(j)=0.0_fp
        i_wbav0(k)=j
     endif
  enddo
! if k=mj, no particles in any zone; wbav=0.0 -- leave as is
  if(k.lt.mj) then
     i1=1
     i2=k
     if(i_wbav0(1).eq.1) then
        j=1
        iw_right=i_wbav0(1+j)
        iw_left=i_wbav0(1)
        do while(iw_right-iw_left.eq.1)
           j=j+1
           iw_left= iw_right
           iw_right=i_wbav0(1+j)
        enddo
        wbav(1:j)=wbav(1+j)
        i1=j+1
     endif
     if(i_wbav0(k).eq.mj) then
        j=1
        iw_right=i_wbav0(k)
        iw_left=i_wbav0(k-j)
        do while(iw_right-iw_left.eq.1)
           j=j+1
           iw_right=iw_left
           iw_left=i_wbav0(k-j)
        enddo
        wbav(mj-j+1:mj)=wbav(mj-j)
        i2=k-j
     endif
     j=i1
     do while(i_wbav0(j).gt.0.and.i_wbav0(j).le.i_wbav0(i2))
        iw=i_wbav0(j)
        iw_left=iw-1
        iw_right=i_wbav0(j+1)
        do while(iw_right-iw.eq.1)
           j=j+1
           iw=iw_right
           iw_right=i_wbav0(j+1)
        enddo
        iw_right=iw+1
! fill in interior gap of wbav by linear interpolation
        do i=i1,j
           wbav(i_wbav0(i))=wbav(iw_left)+(xizc(i_wbav0(i))-xizc(iw_left))*&
                & (wbav(iw_right)-wbav(iw_left))/(xizc(iw_right)-xizc(iw_left))
        enddo
        j=j+1
        i1=j
     enddo
  endif
  deallocate(i_wbav0)
end subroutine get_noholes_array

