subroutine aptspod (pow, xa, xb, np, xran, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTSPOD
c
c     call aptspod (pow, xa, xb, np, xran, nerr)
c
c     Version:  aptspod  Updated    1990 May 3 16:50.
c               aptspod  Originated 1990 May 3 16:50.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To randomly sample np values of xran in the interval from xa to
c               xb from the distribution function p(x) = k * x**pow.
c               The values of xa and xb must be non-negative, with xb no smaller
c               than xa.  Flag nerr indicates any input error.
c
c     Input:    pow, xa, xb, np.
c
c     Output:   xran, nerr.
c
c     Glossary:
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if xa .lt. 0.0 or xb .lt. xa.
c
c     np        Input    Size of array xran.
c
c     pow       Input    Power in probability distribution function:
c                          p(x) = k * x**pow.
c                          Negative values are acceptable, including -1.0.
c                          If pow = -1.0, and xa = 0.0, all xran will be 0.0.
c                          Magnitude must not be too large for computer.
c
c     xa, xb    Input    Limits of range for randomly sampling values of x.
c                          Restrictions:  xa .ge. 0.0, and xb .ge. xa.
c
c     xran      Output   Randomly sampled values of x in range from xa to xb.
c                          Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Randomly sampled value of x.
      dimension xran    (1)

c.... Local variables.

c---- Function xbpow - xapow.
      common /laptspod/ dxpow

c---- A very small number.
      common /laptspod/ fuz

c---- Index in local array.
      common /laptspod/ n
c---- First index of subset of data.
      common /laptspod/ n1
c---- Last index of subset of data.
      common /laptspod/ n2
c---- Index in external array.
      common /laptspod/ nn
c---- Size of current subset of data.
      common /laptspod/ ns
c---- Function 1.0 / powp.
      common /laptspod/ powin
c---- Function pow + 1.0.
      common /laptspod/ powp
c---- First random number.
      common /laptspod/ ranfp   (64)
c---- Function xa**powp.
      common /laptspod/ xapow
c---- Function xb**powp.
      common /laptspod/ xbpow
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean x.
cbug      common /laptspod/ xdev
cbugc---- Mean value of x.
cbug      common /laptspod/ xmean
cbug 9901 format (/ 'aptspod sampling from a power distribution.  np=',i3 /
cbug     &  '  pow,xa,xb=',1p3e22.14)
cbug      write (3, 9901) np, pow, xa, xb
cbugc***DEBUG ends.

c.... Initialize.

c---- A very small number.
      fuz = 1.e-99

c.... Test for input errors.

      nerr  = 0
      if (np .le. 0) then
        nerr = 1
        go to 210
      endif

      if ((xa .lt. 0.0) .or. (xb .lt. xa)) then
        nerr = 3
        go to 210
      endif

c.... initialize.

      powp  = pow + 1.0
      powin = 1.0 / (powp + fuz)
      xapow = xa**powp
      xbpow = xb**powp
      dxpow = xbpow - xapow

c.... Set up the indices of the first subset of data.

      n1 = 1
      n2 = min (np, 64)
  110 ns = n2 - n1 + 1

c.... Generate ns random numbers.

c---- Loop over subset of data.
      do 120 n = 1, ns
        ranfp(n) = ranf( )
c---- End of loop over subset of data.
  120 continue

c.... Sample values of x from the distribution p(x).

c---- Distribution not k / x.
      if (abs (powp) .gt. 0.01) then

c---- Loop over subset of data.
        do 130 n = 1, ns

          nn       = n + n1 - 1
          xran(nn) = (xapow + ranfp(n) * dxpow)**powin

c---- End of loop over subset of data.
  130   continue

c---- Use distribution p(x) = k / x.
      else

c---- Loop over subset of data.
        do 140 n = 1, ns

          nn       = n + n1 - 1
          xran(nn) = xa**ranfp(n) * xb**(1.0 - ranfp(n))

c---- End of loop over subset of data.
  140   continue

c---- Tested powp.
      endif

c.... See if all data subsets are done.

c---- Do another subset of data.
      if (n2 .lt. np) then
        n1 = n2 + 1
        n2 = min (np, n1 + 63)
        go to 110
      endif
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptspod results:  xran')
cbug      write ( 3, 9903)
cbug      call aptmean (xran, np, 1.e-11, xmean, xdev, nerr)
cbugc***DEBUG ends.

  210 return

c.... End of subroutine aptspod.      (+1 line.)
      end

UCRL-WEB-209832