subroutine aptxnup (nopt, tplanck, np, xnu, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTXNUP
c
c call aptxnup (nopt, tplanck, np, xnu, nerr)
c
c Version: aptxnup Updated 1990 December 3 16:20.
c aptxnup Originated 1990 January 31 17:00.
c
c Authors: Eugene H. Canfield, LLNL, L-298, Telephone (925) 422-4125.
c Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np frequencies xnu from a Planck or a Wien
c spectrum, for np temperatures tplanck.
c Flag nerr indicates any input error (np not positive).
c
c The Planck distribution function is:
c p(x) = (15.0 / pi**4) * x**3 / (exp (x) - 1.0),
c where x = xnu / tplanck, and = 3.83223...
c The Wien distribution function is:
c p(x) = x**3 * exp (-x) / 6.0.
c where = 4.0.
c
c History: See Barnett and Canfield, UCIR-473, June 1970.
c 1990 February 6 11:20. Truncated pifn, to limit number of terms
c required to converge.
c
c Input: nopt, tplanck, np.
c
c Output: xnu, nerr.
c
c Glossary:
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c nopt Input Indicates type of spectrum to sample from:
c 0 for a Planck spectrum.
c 1 for a Wien spectrum.
c
c np Input Size of arrays tplanck, xnu.
c
c tplanck Input Black body temperature of frequency distribution.
c Same units as xnu. Size np.
c
c xnu Output Frequency sampled randomly from spectrum. Size np.
c Same units as tplanck.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Blackbody temperature of spectrum.
dimension tplanck (1)
c---- Randomly sampled frequency.
dimension xnu (1)
c.... Local variables.
c---- Index in local array.
common /laptxnup/ n
c---- First index of subset of data.
common /laptxnup/ n1
c---- Last index of subset of data.
common /laptxnup/ n2
c---- Size of current subset of data.
common /laptxnup/ ns
c---- Constant pi**4 / 90.0, truncated.
common /laptxnup/ pifn
c---- An integer: 1, 2, 3, ...
common /laptxnup/ pln
c---- Reciprocal of an integer.
common /laptxnup/ plnr
c---- A random function.
common /laptxnup/ plrn
c---- 1 + (1/2)**4 + (1/3)**4 + ...
common /laptxnup/ plsum
c---- Random number (0.0 to 1.0).
common /laptxnup/ rn1
c---- Random number (0.0 to 1.0).
common /laptxnup/ rn2
c---- Random number (0.0 to 1.0).
common /laptxnup/ rn3
c---- Random number (0.0 to 1.0).
common /laptxnup/ rn4
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean value.
cbug common /laptxnup/ xnudev
cbugc---- Mean value of xnu.
cbug common /laptxnup/ xnumean
cbug 9901 format (/ 'aptxnup finding frequencies. nopt=',i2,' np=',i3 /
cbug & (i3,' tplanck=',1pe22.14))
cbug write ( 3, 9901) nopt, np, (n, tplanck(n), n = 1, np)
cbugc***DEBUG ends.
c.... Test for input errors.
nerr = 0
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Initialize. Truncate pifn, to limit number of terms in distribution.
c---- Constant pi**4 / 90.0, truncated.
pifn = 1.082323
c.... Do data sets in groups of 64.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Sample frequencies from a Planck or Wien spectrum.
c---- Loop over subset of data.
do 140 n = 1, ns
pln = 1.0
plnr = 1.0
plsum = 1.0
plrn = pifn * ranf( )
if (nopt .ne. 0) then
plrn = plrn - 1.0
endif
120 if (plsum .gt. plrn) go to 130
pln = pln + 1.0
plnr = 1.0 / pln
plsum = plsum + plnr**4
go to 120
130 rn1 = ranf( )
rn2 = ranf( )
rn3 = ranf( )
rn4 = ranf( )
xnu(n+n1-1) = -plnr * tplanck(n) *
& alog (rn1 * rn2 * rn3 * rn4)
c---- End of loop over subset of data.
140 continue
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
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptxnup results: xnu')
cbug write ( 3, 9902)
cbug call aptmean (xnu, np, 1.e-11, xnumean, xnudev, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptxnup. (+1 line.)
end
UCRL-WEB-209832