module nbi_random

  use iso_c_binding, only: fp=>c_double
  use rng
  
  implicit none
  private :: fp
  save
  
  integer :: ilunrng = 6  ! for messages
  
  integer, dimension(rng_s) :: rng_seed    ! rng seed
  type(rng_state) :: rng_curstate             ! rng current state

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

  real(fp), dimension(rng_batch_size) :: nbi_srng    ! random no.s
  real(fp), dimension(rng_batch_size) :: nbi_sgauss  ! random no.s from gaussian

  integer :: nbi_srng_ptr = 0          ! points to current random no.
  integer :: nbi_sgauss_ptr = 0        ! points to current r.gauss no.

contains

  real(fp) function nbrandm(idum)
    
    integer idum              ! dummy argument (artifact)
    
    if(nbi_srng_ptr.eq.0) then
       write(ilunrng,*) ' ?nbrandm:  nbi_random module not initialized!'
       nbrandm=0.0_fp
    else
       nbrandm=nbi_srng(nbi_srng_ptr)
       nbi_srng_ptr=nbi_srng_ptr+1
       if(nbi_srng_ptr.gt.rng_batch_size) then
          call rng_number(rng_curstate,nbi_srng)
          nbi_srng_ptr=1
       endif
    endif
    
  end function nbrandm
  
  real(fp) function nbgauss(zs,za)
    
    real(fp),intent(in) :: zs,za
    
    !c       this function samples from a gaussian of the
    !c       form exp(-(x-za)**2/(2.*zs**2))
    
    if(nbi_sgauss_ptr.eq.0) then
       write(ilunrng,*) ' ?nbgauss:  nbi_random module not initialized!'
       nbgauss=0.0_fp
    else
       nbgauss=zs * nbi_sgauss(nbi_sgauss_ptr) + za
       nbi_sgauss_ptr=nbi_sgauss_ptr+1
       if(nbi_sgauss_ptr.gt.rng_batch_size) then
          call rng_gauss(rng_curstate,nbi_sgauss)
          nbi_sgauss_ptr=1
       endif
    endif
    
  end function nbgauss
  
end module nbi_random
