!-----------------------------------------------------------------------
!  bbnpre -- prepare for monte carlo beam-beam neutron calculation
!
!  mod dmc 11 jan 1995 -- can now have situation where one dist fcn is
!  from mc and one from fpp.  save the mc parameters in temporary arrays
!  (passed) for later restoration.
!
!-----------------------------------------------------------------------
!  mod dmc 24 sep 1991 -- energy grid of fast ion dist. fcns now varies
!  with fast ion species, and it may not be divided up in constant steps
!  of d(energy) -- take this into account in the code
!-----------------------------------------------------------------------
!
!    dmc june 1990 -- this routine is used when the beam distribution
!    function is provided from fpp -- i.e. it is not sparse and has
!    wide dynamic range, s.t. the monte-carlo discretization must be
!    carefully done.  because of the large no. of small but non-zero
!    points expected (0 < f(i,j) << fmax), a cumulative distribution
!    rather than a rejection method will be used for the selection of
!    representative points.
!
!    this routine will compute an array of the number of points to
!    be used to represent the distributions in the beam-beam neutron
!    monte carlo integration
!
!    input namelist integer :: nbbcal = total no. of collisions desired
!    to model beam-beam neutrons over entire plasma
!
! dec.18, 2008
! bbnpre divided to bbnpre_bdens2(is1,is2)
! to   calculate bdens2 and nfbmpt for isb1, isb2
! if nlfbmflr=.true. then distribution function at partical location will be used
! if nlfbmflr=.false. then distribution function at gc will be used
! and  to bbnpre_ncoll(is1,is2) to calculate  total no. of collisions desired
!    to model beam-beam neutrons over entire plasma


subroutine bbnpre_nevents(is1,is2)
  !
  !  is1 -- index of first reagent species
  !  is2 -- index of 2nd reagent species
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: is2,is1,i,ichk,inum
  !============
  real(fp) :: znnv,zdnnv,znum2
  !============
  integer :: izntot
  !  output -- the # of pts to represent the local distribution functions,
  !    in nfbmpt(i,is1) for zone i
  ! -----------------------------------------------------
  !
  ! sum n**2*dv
  !
  !
  !
  ! compute the no. of collisions in each zone to shoot for...
  izntot=mimxba*mimxbe
  znnv=0.0_fp
  do i =1,nfbznsi
     znnv=znnv+(bdens2(i,is1)/bmvol(i))*bdens2(i,is2)
  end do
  ichk=0
  do i=1,nfbznsi
     nfbmpt(i,is1)=0
     zdnnv=(bdens2(i,is1)/bmvol(i))*bdens2(i,is2)
     if(zdnnv.eq.0.0_fp) cycle
     !
     znum2=real(nbbcal)*zdnnv/znnv
     !
     inum=1+int(sqrt(znum2))
     inum=min(izntot,max(10,inum))
     !
     nfbmpt(i,is1)=inum
     ichk=ichk+inum*inum
  end do
  !
  return
end subroutine bbnpre_nevents

subroutine bbnpre_bdens2(is1,is2)
  ! dec.18, 2008
  ! bbnpre divided to bbnpre_bdens2(is1,is2)
  ! to   calculate bdens2 and nfbmpt for isb1, isb2 species
  ! if nlfbmflr=.true. then distribution function at partical location will be used
  ! if nlfbmflr=.false. then distribution function at gc will be used
  ! bdens2 here is different from one in outptb2 on factor delt
  !
  !     is1 -- index of first reagent species
  !     is2 -- index of 2nd reagent species
  !
  use nbi_com
  !
  implicit none
  integer :: is2,is1,i,ia,ie
  real(fp) :: zdelebd,zdnw
  real(fp) :: fbm_flr            !function to choose distribution function
  ! depending on nlfbmflr, fbm if nlfbmflr=.false.
  ! fbm_ptcl if nlfbmflr=.true.
  !
  do i=1,nfbznsi
     !------------------
     bdens2(i,is1)=0.0_fp
     bdens2(i,is2)=0.0_fp
     nfbmpt(i,is1)=0.0_fp
     nfbmpt(i,is2)=0.0_fp

     do ia=1,nznbma
        do ie=1,nznbmea(is1)
           zdelebd=efbmb(ie+1,is1)-efbmb(ie,is1)
           zdnw=fbm_flr(ie,ia,i,is1,0)
           if(zdnw.gt.0.0_fp) then
              nfbmpt(i,is1)=nfbmpt(i,is1) + 1
              bdens2(i,is1)=bdens2(i,is1)+ &
                   zdnw*bmvol(i)*zdelebd* &
                   delt/real(nznbma)
           end if
        end do
        if(is1.ne.is2) then
           do ie=1,nznbmea(is2)
              zdelebd=efbmb(ie+1,is2)-efbmb(ie,is2)
              zdnw=fbm_flr(ie,ia,i,is2,0)
              if(zdnw.gt.0.0_fp) then
                 nfbmpt(i,is2)=nfbmpt(i,is2) + 1
                 bdens2(i,is2)=bdens2(i,is2)+ &
                      zdnw*bmvol(i)*zdelebd* &
                      delt/real(nznbma)
              end if
           end do
        else
           nfbmpt(i,is2)=nfbmpt(i,is1)
           bdens2(i,is2)=bdens2(i,is1)
        end if
     end do
  end do
  !     bdens2 has a *dv in it
  return
end subroutine bbnpre_bdens2







