module nbi_random_dep

  use iso_c_binding, only: fp => c_double
  use rng

  implicit none
  private :: fp
  save

  integer :: ilunrng_dep = 6  ! for messages

  integer, dimension(rng_s) :: rng_seed_dep    ! rng seed
  type(rng_state) :: rng_curstate_dep             ! rng current state

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

  real(fp), dimension(rng_batch_size_dep) :: nbi_srng_dep    ! random no.s
  real(fp), dimension(rng_batch_size_dep) :: nbi_sgauss_dep  ! random no.s from gaussian

  integer :: nbi_srng_ptr_dep = 0          ! points to current random no.
  integer :: nbi_sgauss_ptr_dep = 0        ! points to current r.gauss no.

  contains

    real(fp) function nbrandm_dep(idum_dep)

      integer idum_dep              ! dummy argument (artifact)

      if(nbi_srng_ptr_dep.eq.0) then
         write(ilunrng_dep,*) ' ?nbrandm_dep:  nbi_random_dep module not initialized!'
         nbrandm_dep=0.0_fp
      else
         nbrandm_dep=nbi_srng_dep(nbi_srng_ptr_dep)
         nbi_srng_ptr_dep=nbi_srng_ptr_dep+1
         if(nbi_srng_ptr_dep.gt.rng_batch_size_dep) then
            call rng_number(rng_curstate_dep,nbi_srng_dep)
            nbi_srng_ptr_dep=1
         endif
      endif

    end function nbrandm_dep
 
    real(fp) function nbgauss_dep(zs_dep,za_dep)

      real(fp),intent(in) :: zs_dep,za_dep

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

      if(nbi_sgauss_ptr_dep.eq.0) then
         write(ilunrng_dep,*) ' ?nbgauss_dep:  nbi_random_dep module not initialized!'
         nbgauss_dep=0.0_fp
      else
         nbgauss_dep=zs_dep * nbi_sgauss_dep(nbi_sgauss_ptr_dep) + za_dep
         nbi_sgauss_ptr_dep=nbi_sgauss_ptr_dep+1
         if(nbi_sgauss_ptr_dep.gt.rng_batch_size_dep) then
            call rng_gauss(rng_curstate_dep,nbi_sgauss_dep)
            nbi_sgauss_ptr_dep=1
         endif
      endif

    end function nbgauss_dep
            
end module nbi_random_dep
