module nbi_random_aep

  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_aep    ! rng seed
  type(rng_state) :: rng_curstate_aep             ! rng current state

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

  real(fp), dimension(rng_batch_size_aep) :: nbi_srng_aep    ! random no.s
  real(fp), dimension(rng_batch_size_aep) :: nbi_sgauss_aep  ! random no.s from gaussian

  integer :: nbi_srng_ptr_aep = 0          ! points to current random no.
  integer :: nbi_sgauss_ptr_aep = 0        ! points to current r.gauss no.

contains

  real(fp) function nbrandm_aep(idum_aep)

    integer idum_aep              ! dummy argument (artifact)

    if(nbi_srng_ptr_aep.eq.0) then
       write(ilunrng,*) ' ?nbrandm_aep:  nbi_random_aep module not initialized!'
       nbrandm_aep=0.0_fp
    else
       nbrandm_aep=nbi_srng_aep(nbi_srng_ptr_aep)
       nbi_srng_ptr_aep=nbi_srng_ptr_aep+1
       if(nbi_srng_ptr_aep.gt.rng_batch_size_aep) then
          call rng_number(rng_curstate_aep,nbi_srng_aep)
          nbi_srng_ptr_aep=1
       endif
    endif
    
  end function nbrandm_aep
  
  real(fp) function nbgauss_aep(zs_aep,za_aep)
    
    real(fp),intent(in) :: zs_aep,za_aep
    
    !c       this function samples from a gaussian of the
    !c       form exp(-(x-za)**2/(2.*zs**2))
    
    if(nbi_sgauss_ptr_aep.eq.0) then
       write(ilunrng,*) ' ?nbgauss_aep:  nbi_random_aep module not initialized!'
       nbgauss_aep=0.0_fp
    else
       nbgauss_aep=zs_aep * nbi_sgauss_aep(nbi_sgauss_ptr_aep) + za_aep
       nbi_sgauss_ptr_aep=nbi_sgauss_ptr_aep+1
       if(nbi_sgauss_ptr_aep.gt.rng_batch_size_aep) then
          call rng_gauss(rng_curstate_aep,nbi_sgauss_aep)
          nbi_sgauss_ptr_aep=1
       endif
    endif
    
  end function nbgauss_aep
  
end module nbi_random_aep
