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