c===========================================================
c     nurand:  Uses SGI-specific uniform-random number 
c     generator rand() to generate non-uniformly generated 
c     random numbers on the interval [xmin..xmax]. User
c     must supply probability distribution function 
c     having a header
c
c        subroutine pdf(x,pofx,maxpofx)
c
c     where 'x' is the input value, 'pofx' is the value 
c     of the PDF evaluated at 'x' and 'maxpofx' is the 
c     maximum value of the PDF (also a return argument).
c
c     Uses straight-forward algorithm based on area-under-
c     curve (PDF) idea--i.e. generate random point in 
c     rectangle [xmin..xmax] x [0..maxpofx], accept point 
c     and return x coordinate of point as random number
c     only if random point lies below PDF curve.
c===========================================================
      double precision function nurand(pdf,xmin,xmax)
         implicit      none

         external      pdf
         real*8        rand

         real*8        xmin,        xmax

         real*8        x,           y, 
     &                 pofx,        maxpofx

c-----------------------------------------------------------
c        Loop until a good deviate has been generated.
c        Note that we exit the loop via the 'return' 
c        statement---potentially this could be an infinite
c        loop, so for a "production" routine, it might
c        be wise to limit the number of iterations. 
c-----------------------------------------------------------
         do while( .true. ) 
c-----------------------------------------------------------
c           Generate a uniform number in the interval xmin 
c           to xmax.
c-----------------------------------------------------------
            x  =  xmin + rand() * (xmax - xmin)
c-----------------------------------------------------------
c           Evaluate PDF at x.
c-----------------------------------------------------------
            call pdf(x,pofx,maxpofx)
c-----------------------------------------------------------
c           Generate another uniform number in the interval
c           0 to maxpofx ... 
c-----------------------------------------------------------
            y = rand() * maxpofx
c-----------------------------------------------------------
c           ... and accept the original random number, x,
c           if y < pofx.
c-----------------------------------------------------------
            if( y .lt. pofx ) then
               nurand = x
               return
            end if
         end do

      end