module nbi_random_pll

  use iso_c_binding, only: fp => c_double
  use rng

  implicit none
  private :: fp
  save

  integer :: ilunrng_pll = 6  ! for messages

  integer, dimension(rng_s) :: rng_seed_pll    ! rng seed
  type(rng_state) :: rng_curstate_pll             ! rng current state

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

  real(fp), dimension(rng_batch_size_pll) :: nbi_srng_pll    ! random no.s
  real(fp), dimension(rng_batch_size_pll) :: nbi_sgauss_pll  ! random no.s from gaussian

  integer :: nbi_srng_ptr_pll = 0          ! points to current random no.
  integer :: nbi_sgauss_ptr_pll = 0        ! points to current r.gauss no.

  contains

    real(fp) function nbrandm_pll(idum_pll)

      integer idum_pll              ! dummy argument (artifact)

      if(nbi_srng_ptr_pll.eq.0) then
         write(ilunrng_pll,*) ' ?nbrandm_pll:  nbi_random_pll module not initialized!'
         nbrandm_pll=0.0_fp
      else
         nbrandm_pll=nbi_srng_pll(nbi_srng_ptr_pll)
         nbi_srng_ptr_pll=nbi_srng_ptr_pll+1
         if(nbi_srng_ptr_pll.gt.rng_batch_size_pll) then
            call rng_number(rng_curstate_pll,nbi_srng_pll)
            nbi_srng_ptr_pll=1
         endif
      endif

    end function nbrandm_pll
 
    real(fp) function nbgauss_pll(zs_pll,za_pll)

      real(fp),intent(in) :: zs_pll,za_pll

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

      if(nbi_sgauss_ptr_pll.eq.0) then
         write(ilunrng_pll,*) ' ?nbgauss_pll:  nbi_random_pll module not initialized!'
         nbgauss_pll=0.0_fp
      else
         nbgauss_pll=zs_pll * nbi_sgauss_pll(nbi_sgauss_ptr_pll) + za_pll
         nbi_sgauss_ptr_pll=nbi_sgauss_ptr_pll+1
         if(nbi_sgauss_ptr_pll.gt.rng_batch_size_pll) then
            call rng_gauss(rng_curstate_pll,nbi_sgauss_pll)
            nbi_sgauss_ptr_pll=1
         endif
      endif

    end function nbgauss_pll
            
end module nbi_random_pll
