subroutine aptexpl (nopt, sigma, xa, xb, np, xran, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTEXPL
c
c call aptexpl (nopt, sigma, xa, xb, np, xran, nerr)
c
c Version: aptexpl Updated 1994 June 27 16:20.
c aptexpl Originated 1994 June 27 16:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample each of np values xran(n), in the interval
c from xa(k) to xb(k), from the unnormalized exponential
c differential probability distribution function (PDF) with decay
c constant sigma(k), where k = 1 (nopt = 0), or k = n (nopt = 1):
c
c PDF = exp (-sigma(k) * x), xa(k) <= x <= xb(k)
c
c Flag nerr indicates any input error.
c
c Input: nopt, sigma, xa, xb, np.
c
c Output: xran, nerr.
c
c Glossary:
c
c nerr Output Indicates an input error, if not zero.
c 1 if nopt is not 0 or 1.
c 2 if np is not positive.
c
c nopt Input Indicates the sampling option:
c 0: All np xran values will be randomly sampled from
c one exponential probability distribution function
c with decay constant sigma (a scalar).
c 1: Each value xran(n) will be randomly sampled from
c an exponential probability distribution function
c with decay constant sigma(n) (an array).
c
c np Input Size of arrays xran, xa and xb, and if nopt = 1, sigma.
c
c sigma Input Decay constant. Array size 1 (nopt = 0) or np
c (nopt = 1).
c If x is a collision distance, sigma is a cross
c section (per unit distance) or 1 / mean-free-path.
c If x is a time, sigma is log (2) / half-life.
c If xa or xb is very large and positive, sigma must
c not be negative.
c
c xa Input Lower limit on range of xran. Array size np.
c Array size 1 (nopt = 0) or np (nopt = 1).
c
c xb Output Upper limit on range of xran. Use a very large number
c such as 40 / sigma, to approximate infinity.
c Array size 1 (nopt = 0) or np (nopt = 1).
c
c xran Output Randomly sampled values of x. For xran(n), from the
c exponential probability distribution function with
c decay constant sigma, between xa(k) and xb(k), where
c k = 1 (nopt = 0), or k = n (nopt = 1).
c Array size np.
c The expected value of xran is
c xa + (1 - f/(exp(f) - 1)) / sigma,
c with a standard deviation of
c sqrt (1 - f**2*exp(f)/(exp(f) - 1)**2) / sigma,
c where f = sigma * (xb - xa),
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Decay constant.
dimension sigma (1)
c---- Randomly sampled value of x.
dimension xa (1)
dimension xb (1)
dimension xran (1)
c.... Local variables.
c---- Index in local array.
common /laptexpl/ n
c---- First index of subset of data.
common /laptexpl/ n1
c---- Last index of subset of data.
common /laptexpl/ n2
c---- Index in external array.
common /laptexpl/ nn
c---- Size of current subset of data.
common /laptexpl/ ns
c---- Random numbers.
common /laptexpl/ ranfp (64)
c---- Argument of log function.
common /laptexpl/ arg (64)
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean.
cbug common /laptexpl/ xdev
cbugc---- Mean value.
cbug common /laptexpl/ xmean
cbug 9901 format (/ 'aptexpl sampling an exponential distribution.',
cbug & ' nopt=',i2,' np=',i5)
cbug 9902 format (' sigma= ',1pe22.14 /
cbug & ' xa, xb=',1p2e22.14 )
cbug 9903 format (i5,' sigma= ',1pe22.14 /
cbug & ' xa, xb=',1p2e22.14 )
cbug write ( 3, 9901) nopt, np
cbug if (nopt .eq. 0) then
cbug write ( 3, 9902) sigma(1), xa(1), xb(1)
cbug elseif (nopt .eq. 1) then
cbug write (3, 9903) (n, sigma(n), xa(n), xb(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
c.... Test for input errors.
nerr = 0
if ((nopt .ne. 0) .and. (nopt .ne. 1)) then
nerr = 1
go to 210
endif
if (np .le. 0) then
nerr = 2
go to 210
endif
c.... Do data sets in groups of 64.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Get needed 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.... Randomly sample values of x.
c---- One PDF for all xran.
if (nopt .eq. 0) then
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
if (sigma(1) .ne. 0.0) then
arg(n) = ranfp(n) + (1.0 - ranfp(n)) *
& exp (-sigma(1) * (xb(1) - xa(1)))
xran(nn) = xa(1) - alog (arg(n)) / sigma(1)
else
xran(nn) = xa(1) + ranfp(n) * (xb(1) - xa(1))
endif
c---- End of loop over subset of data.
130 continue
c---- One PDF for each xran.
else
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
if (sigma(nn) .ne. 0.0) then
arg(n) = ranfp(n) + (1.0 - ranfp(n)) *
& exp (-sigma(nn) * (xb(nn) - xa(nn)))
xran(nn) = xa(nn) - alog (arg(n)) / sigma(nn)
else
xran(nn) = xa(nn) + ranfp(n) * (xb(nn) - xa(nn))
endif
c---- End of loop over subset of data.
140 continue
c---- Tested nopt.
endif
c.... See if all sets are done.
c---- Do another data set.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
210 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptexpl results: nerr = ',i3,' xran:')
cbug write ( 3, 9904) nerr
cbug call aptmean (xran, np, 1.e-11, xmean, xdev, nerr)
cbugc***DEBUG ends.
return
c.... End of subroutine aptexpl. (+1 line.)
end
UCRL-WEB-209832