c===========================================================
c     usage: tnurand <xmin> <xmax> <n> [<nbin> <option>]
c===========================================================
c     tnurand: Driver program for nurand(). This driver
c     generates non-uniformly distributed random-numbers 
c     using a user-specified distribution function.  The 
c     program is currently set up with two distribution 
c     functions (see 'pdfs.f'):
c
c     option = 0            --> uniform
c     option = 1 (default)  --> unit-sigma Gaussian
c
c     The routine calls nurand() to generate n random
c     numbers, then writes binned counts (the interval
c     xmin ... xmax is divided into nbin equal width bins)
c
c       <i>   <count i>
c
c     i = 1 ... nbin, on standard output.
c
c     Note that nurand() uses rand(), so srand() can be 
c     called to "seed" nurand().
c===========================================================
      program         tnurand

      implicit        none

c-----------------------------------------------------------
c     External declarations for the user-defined PDFs and 
c     declaration of nurand.
c-----------------------------------------------------------
      external        puniform,     pgauss
      real*8          nurand

      integer         iargc,        i4arg
      real*8          r8arg
      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )
c-----------------------------------------------------------
c     Command-line arguments:
c
c     xmin:   Minimum, maximum values of deviates
c     xmax:  
c     n:      Number of deviates to generate
c     nbin:   Number of binning intervals
c     option: Selects probability distribution function
c-----------------------------------------------------------
      real*8          xmin,         xmax
      integer         n,            nbin,         option

      integer         max_nbin
      parameter     ( max_nbin = 10 000 )
      real*8          x(max_nbin),  count(max_nbin)

      real*8          dx,           rnum
      integer         i,            j

c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .lt. 1 ) go to 900

      xmin   = r8arg(1,r8_never)
      if( xmin .eq. r8_never ) go to 900
      xmax   = r8arg(2,r8_never)
      if( xmax .eq. r8_never ) go to 900
      n      = i4arg(3,-1)
      if( n .le. 0 ) go to 900
      nbin   = min(i4arg(4,1000),max_nbin)
      option = i4arg(5,1)

c-----------------------------------------------------------
c     Set up bins and bin-coordinates (mid-points of bin
c     intervals).
c-----------------------------------------------------------
      dx = (xmax - xmin) / nbin 
      do i = 1 , nbin
         count(i) = 0.0d0
         if( i .eq. 1 ) then
            x(1) = xmin + 0.5d0 * dx
         else
            x(i) = x(i-1) + dx
         end if
      end do

c-----------------------------------------------------------
c     Generate and bin random numbers.
c-----------------------------------------------------------
      do i = 1 , n
         if(      option .eq. 0 ) then
            rnum = nurand(puniform,xmin,xmax)
         else if( option .eq. 1 ) then
            rnum = nurand(pgauss,xmin,xmax)
         else 
            write(0,*) 'tnurand: Unimplemented option ', 
     &                 option
            stop
         end if
         j = min(int((rnum - xmin) / dx) + 1,nbin)
         count(j) = count(j) + 1.0d0
      end do

c-----------------------------------------------------------
c     Normalize bin counts.
c-----------------------------------------------------------
      do i = 1 , nbin
         count(i) = count(i) / (dx * n)
      end do

c-----------------------------------------------------------
c     Output bin counts.
c-----------------------------------------------------------
      call dvvto('-',x,count,nbin)

      stop

 900  continue
         write(0,*) 'tnurand: <xmin> <xmax> <n> '//
     &              '[<nbin> <option>]'
      stop

      end

c-----------------------------------------------------------
c     Writes two double precision vectors to file 'fname'.
c     If fname is the string '-' then the vectors are 
c     written to standard output.
c-----------------------------------------------------------

      subroutine dvvto(fname,v1,v2,n)
         implicit         none

         integer          getu,       indlnb

         character*(*)    fname
         integer          n
         real*8           v1(n),      v2(n)

         integer          ustdout
         parameter      ( ustdout = 6 )

         integer          i,      uto,       rc

         if( fname .eq. '-' ) then
            uto = ustdout
         else
            uto = getu()
            open(uto,file=fname(1:indlnb(fname)),
     &           form='formatted',iostat=rc)
            if( rc .ne. 0 ) then
               write(0,*) 'dvvto: Error opening ',
     &                    fname(1:indlnb(fname))
               return
            end if
         end if

         do i = 1 , n
            write(uto,*) v1(i),     v2(i)
         end do

         if( uto .ne. ustdout ) then
            close(uto)
         end if

         return

      end