subroutine aptnorm (nopt, xmean, xdev, np, xran, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTNORM
c
c call aptnorm (nopt, xmean, xdev, np, xran, nerr)
c
c Version: aptnorm Updated 1990 August 15 16:00.
c aptnorm Originated 1990 August 15 16:00.
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) from a Normal
c probability distribution function (PDF) with a mean value
c xmean(k) and a standard deviation xdev(k), where k = 1
c (nopt = 0), or k = n (nopt = 1):
c
c PDF = b * exp (-0.5 * y**2)
c
c where x ranges from minus infinity to plus infinity,
c b = 1.0 / (sqrt (2.0 * pi) * xdev), y = (x - xmean) / xdev,
c xdev = sqrt ( - **2), and <> indicates expected value.
c
c Flag nerr indicates any input error.
c
c Input: nopt, xmean, xdev, 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 the sampling option:
c 0: All np xran values will be randomly sampled from
c one Normal probability distribution function
c based on xmean and xdev (scalars).
c 1: Each value xran(n) will be randomly sampled from
c a Normal probability distribution function
c based on xmean(n) and xdev(n) (arrays).
c
c np Input Size of arrays xran, and if nopt = 1, xmean and xdev.
c
c xdev Input Standard deviation of x from xmean:
c xdev = sqrt ( - **2).
c Size 1 (nopt = 0) or np (nopt = 1).
c
c xmean Input Mean or expected value of x:
c xmean = .
c Size 1 (nopt = 0) or np (nopt = 1).
c
c xran Output Randomly sampled values of x. For xran(n), from
c the Normal probability distribution function with
c mean value xmean(k) and standard deviation xdev(k),
c where k = 1 (nopt = 0) or k = n (nopt = 1).
c Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Standard deviation of x from xmean.
dimension xdev (1)
c---- Mean value of x.
dimension xmean (1)
c---- Randomly sampled value of x.
dimension xran (1)
c.... Local variables.
c---- Index in local array.
common /laptnorm/ n
c---- First index of subset of data.
common /laptnorm/ n1
c---- Last index of subset of data.
common /laptnorm/ n2
c---- Function np - n1 + 2.
common /laptnorm/ nf1
c---- Function n1 - 1.
common /laptnorm/ nf2
c---- Index in xran.
common /laptnorm/ nn1
c---- Index in xran.
common /laptnorm/ nn2
c---- Function (np + 1) / 2.
common /laptnorm/ nph
c---- Size of current subset of data.
common /laptnorm/ ns
c---- Numerical constant, pi.
common /laptnorm/ pi
c---- A random variable.
common /laptnorm/ rad
c---- Random numbers.
common /laptnorm/ ranfp1 (64)
c---- Random numbers.
common /laptnorm/ ranfp2 (64)
c---- A random angle.
common /laptnorm/ theta
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean.
cbug common /laptnorm/ xdevr
cbugc---- Mean value.
cbug common /laptnorm/ xmeanr
cbug 9901 format (/ 'aptnorm sampling a normal distribution.',
cbug & ' nopt=',i2,' np=',i5)
cbug 9902 format (' xmean=',1pe22.14,' xdev=',1pe22.14)
cbug 9903 format (i5,' xmean,xdev=',1p2e22.14)
cbug write ( 3, 9901) nopt, np
cbug if (nopt .eq. 0) then
cbug write ( 3, 9902) xmean(1), xdev(1)
cbug elseif (nopt .eq. 1) then
cbug write (3, 9903) (n, xmean(n), xdev(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
c---- Numerical constant, pi.
pi = 3.14159265358979323
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.... Initialize.
nph = (np + 1) / 2
c.... Do data sets in groups of 64.
n1 = 1
n2 = min (nph, 64)
110 ns = n2 - n1 + 1
nf1 = np + 2 - n1
nf2 = n1 - 1
c.... Get needed random numbers.
c---- Loop over subset of data.
do 120 n = 1, ns
ranfp1(n) = ranf( )
c---- End of loop over subset of data.
120 continue
c---- Loop over subset of data.
do 130 n = 1, ns
ranfp2(n) = ranf( )
c---- End of loop over subset of data.
130 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 140 n = 1, ns
nn1 = nf2 + n
nn2 = nf1 - n
rad = sqrt (-2.0 * log (ranfp1(n)))
theta = 2.0 * pi * ranfp2(n)
xran(nn1) = xmean(1) + xdev(1) * rad * cos (theta)
xran(nn2) = xmean(1) + xdev(1) * rad * sin (theta)
c---- End of loop over subset of data.
140 continue
c---- One PDF for each xran.
else
c---- Loop over subset of data.
do 150 n = 1, ns
nn1 = nf2 + n
nn2 = nf1 - n
rad = sqrt (-2.0 * log (ranfp1(n)))
theta = 2.0 * pi * ranfp2(n)
xran(nn1) = xmean(nn1) + xdev(nn1) * rad * cos (theta)
xran(nn2) = xmean(nn2) + xdev(nn2) * rad * sin (theta)
c---- End of loop over subset of data.
150 continue
c---- Tested nopt.
endif
c.... See if all sets are done.
c---- Do another data set.
if (n2 .lt. nph) then
n1 = n2 + 1
n2 = min (nph, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptnorm results: xran')
cbug write ( 3, 9904)
cbug call aptmean (xran, np, 1.e-11, xmeanr, xdevr, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptnorm. (+1 line.)
end
UCRL-WEB-209832