module nbi_random_nbkill

  use iso_c_binding, only: fp => c_double
  use rng

  implicit none
  private :: fp
  save

  integer :: ilunrng_nbkill = 6  ! for messages

  integer, dimension(rng_s) :: rng_seed_nbkill    ! rng seed
  type(rng_state) :: rng_curstate_nbkill             ! rng current state

  integer, parameter :: rng_batch_size_nbkill = 1000    ! rng batch array sizes

  real(fp), dimension(rng_batch_size_nbkill) :: nbi_srng_nbkill    ! random no.s
  real(fp), dimension(rng_batch_size_nbkill) :: nbi_sgauss_nbkill  ! random no.s from gaussian

  integer :: nbi_srng_ptr_nbkill = 0          ! points to current random no.
  integer :: nbi_sgauss_ptr_nbkill = 0        ! points to current r.gauss no.

  contains

    real(fp) function nbrandm_nbkill(idum_nbkill)

      integer idum_nbkill              ! dummy argument (artifact)

      if(nbi_srng_ptr_nbkill.eq.0) then
         write(ilunrng_nbkill,*) ' ?nbrandm_nbkill:  nbi_random_nbkill module not initialized!'
         nbrandm_nbkill=0.0_fp
      else
         nbrandm_nbkill=nbi_srng_nbkill(nbi_srng_ptr_nbkill)
         nbi_srng_ptr_nbkill=nbi_srng_ptr_nbkill+1
         if(nbi_srng_ptr_nbkill.gt.rng_batch_size_nbkill) then
            call rng_number(rng_curstate_nbkill,nbi_srng_nbkill)
            nbi_srng_ptr_nbkill=1
         endif
      endif

    end function nbrandm_nbkill
 
    real(fp) function nbgauss_nbkill(zs_nbkill,za_nbkill)

      real(fp),intent(in) :: zs_nbkill,za_nbkill

      !c       this function samples from a gaussian of the
      !c       form exp(-(x-za)**2/(2.*zs**2))

      if(nbi_sgauss_ptr_nbkill.eq.0) then
         write(ilunrng_nbkill,*) ' ?nbgauss_nbkill:  nbi_random_nbkill module not initialized!'
         nbgauss_nbkill=0.0_fp
      else
         nbgauss_nbkill=zs_nbkill * nbi_sgauss_nbkill(nbi_sgauss_ptr_nbkill) + za_nbkill
         nbi_sgauss_ptr_nbkill=nbi_sgauss_ptr_nbkill+1
         if(nbi_sgauss_ptr_nbkill.gt.rng_batch_size_nbkill) then
            call rng_gauss(rng_curstate_nbkill,nbi_sgauss_nbkill)
            nbi_sgauss_ptr_nbkill=1
         endif
      endif

    end function nbgauss_nbkill
            
end module nbi_random_nbkill
