subroutine apteqxn (x, pl, pr, nbins, nep, xep, pmin, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTEQXN c c call apteqxn (x, pl, pr, nbins, nep, xep, pmin, nerr) c c Version: apteqxn Updated 1990 August 13 13:50. c apteqxn Originated 1990 August 13 13:50. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find the nep equal probability (EP) bins, each bounded by c xep(n) on the left and xep(n+1) on the right, corresponding to c the nbins probability distribution function (PDF) bins, each c with relative probability pl(n) at x(n) on the left, and c relative probability pr(n) at x(n+1) on the right, varying c linearly. For a histogram PDF, pl(n) = pr(n). Also, to find c the minimum normalized probability pmin of any PDF bin. c Flag nerr indicates any input error. c c Random sampling of the random variable x is done as follows: c c n = 1 + ranf( ) * nep c x = xep(n) + ranf( ) * (xep(n+1) - xep(n)) c c Input: x, pl, pr, nbins, nep. c c Output: xep, pmin, nerr. c c Glossary: c c nbins Input Number of input PDF bins or intervals. c c nep Output Number of output EP bins. c c nerr Output Indicates an input error, if not 0. c 1 if nbins is not positive. c 2 if nep is not positive. c 3 if any pl(n) or pr(n) is negative. c 4 if any x(n+1) is less than x(n). c c pl Input For the PDF bin n, the unnormalized relative c probability of x(n), per unit x interval. c Size nbins. May not be negative for any n. c If the PDF is a histogram, pl(n) = pr(n). c If the PDF is continuous, pl(n) = pr(n-1), for c n = 2, nbins. c c pmin Output The minimum normalized probability of any of the c nbins PDF bins. The reciprocal of pmin is the c minimum value of nep needed to guarantee that no c EP bin contains more than one PDF bin boundary. c c pr Input For the PDF bin n, the unnormalized relative c probability of x(n+1), per unit x interval. c Size nbins. May not be negative for any n. c If the PDF is a histogram, pr(n) = pl(n). c If the PDF is continuous, pr(n) = pl(n+1), for c n = 1, nbins - 1. c c x Input The input PDF bin with index n is bounded at the left c by x(n) and at the right by x(n+1). Size nbins + 1. c For any n, x(n) must not be greater than x(n+1). c c xep Output The output EP bin with index n is bounded at the left c by xep(n) and at the right by xep(n+1). c Size nep + 1. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Left-hand probability in PDF bin. dimension pl (1) c---- Right-hand probability in PDF bin. dimension pr (1) c---- Bounding value of x in PDF bin. dimension x (1) c---- Bounding value of x in EP bin. dimension xep (1) cbugc***DEBUG begins. cbugc---- Final PDF bin used by EP bin. cbug dimension nbep (1000) cbugc***DEBUG ends. c.... Local variables. c---- Probability increment in EP bin. common /lapteqxn/ dpcum c---- A very small number. common /lapteqxn/ fuz c---- Index of PDF or EP bin. common /lapteqxn/ n c---- Index of a PDF bin. common /lapteqxn/ nbin c---- Relative probability of a PDF bin. common /lapteqxn/ pbin c---- Cumulative relative probability. common /lapteqxn/ pcum c---- Cum rel prob at left of PDF bin. common /lapteqxn/ pcuml c---- Cum rel prob at left of EP bin. common /lapteqxn/ pcumle c---- Maximum of pcuml, pcumle. common /lapteqxn/ pcumll c---- Cum rel prob at right of PDF bin. common /lapteqxn/ pcumr c---- Cum rel prob at right of EP bin. common /lapteqxn/ pcumre c---- Relative probability at xll. common /lapteqxn/ pll c---- Slope (pr - pl) / (xR - xL). common /lapteqxn/ slope c---- Maximum of left-side x, xep. common /lapteqxn/ xll cbugc***DEBUG begins. cbugc---- Minimum nep needed. cbug common /lapteqxn/ nepneed cbug 9901 format (/ 'apteqxn finding equal-probability bins.', cbug & ' nep=',i5) cbug 9902 format (i5,' xL,xR=',1p2e22.14 / ' pl,pr=',1p2e22.14) cbug write ( 3, 9901) nep cbug write ( 3, 9902) (n, x(n), x(n+1), pl(n), pr(n), n = 1, nbins) cbugc***DEBUG ends. c.... initialize. c---- A very small number. fuz = 1.e-99 nerr = 0 c.... Test for errors in scalar input. if (nbins .le. 0) then nerr = 1 go to 210 endif if (nep .le. 0) then nerr = 2 go to 210 endif c.... Test for errors in input arrays, and find total unnormalized probability. pcum = 0.0 c---- Loop over input PDF bins. do 120 n = 1, nbins if ((pl(n) .lt. 0.0) .or. (pr(n) .lt. 0.0)) then nerr = 3 go to 210 endif if (x(n+1) .lt. x(n)) then nerr = 4 go to 210 endif pbin = 0.5 * (pl(n) + pr(n)) * (x(n+1) - x(n)) pcum = pcum + pbin if (n .eq. 1) then pmin = pbin else pmin = amin1 (pmin, pbin) endif c---- End of loop over input PDF bins. 120 continue pmin = pmin / pcum c.... Integrate over input PDF bins to find output EP bins. nbin = 0 pcumr = 0.0 pcumre = 0.0 xep(1) = x(1) c---- Loop over output EP bins. do 140 n = 1, nep pcumle = pcumre pcumre = n * pcum / nep c---- Add one or more input PDF bins. 130 if (pcumre .gt. pcumr) then c---- Reached end of input PDF table. if (nbin .eq. nbins) then xep(n+1) = x(nbins+1) cbugc***DEBUG begins. cbug nbep(n) = nbins cbugc***DEBUG ends. go to 140 endif nbin = nbin + 1 pcuml = pcumr pbin = 0.5 * (pl(nbin) + pr(nbin)) * (x(nbin+1) - x(nbin)) pcumr = pcuml + pbin go to 130 c---- Tested current input PDF bin. endif c.... Right side of output EP bin is in this input PDF bin. slope = (pr(nbin) - pl(nbin)) / (x(nbin+1) - x(nbin) + fuz) pcumll = amax1 (pcuml, pcumle) dpcum = pcumre - pcumll xll = amax1 (x(nbin), xep(n)) pll = pl(nbin) + slope * (xll - x(nbin)) xep(n+1) = xll + 2.0 * dpcum / & (pll + sqrt (pll**2 + 2.0 * slope * dpcum)) cbugc***DEBUG begins. cbug nbep(n) = nbin cbugc###fx01 format (' n= ',i5,' xepL= ',1pe22.14,' xepR= ',1pe22.14 / cbugc### & 13x,'pcumle=',1pe22.14,' pcumre=',1pe22.14 / cbugc### & 13x,'dpcum= ',1pe22.14,' xll= ',1pe22.14) cbugc###fx02 format (' nbin=',i5,' xL= ',1pe22.14,' xR= ',1pe22.14 / cbugc### & 13x,'pcuml= ',1pe22.14,' pcumr= ',1pe22.14) cbugc### write ( 3, fx01) n, xep(n), xep(n+1), pcumle, pcumre, dpcum, xll cbugc### write ( 3, fx02) nbin, x(nbin), x(nbin+1), pcuml, pcumr cbugc***DEBUG ends. c---- End of loop over output EP bins. 140 continue cbugc***DEBUG begins. cbug 9903 format (/ 'apteqxn results. pmin=',1pe22.14,' nepneed=',i5) cbug 9904 format (i5,' xepL=',1pe22.14) cbug 9905 format (i5,' xepR=',1pe22.14,' pcumre=',1pe22.14) cbug nepneed = 1.0 / (pmin + fuz) cbug write ( 3, 9903) pmin, nepneed cbug n = 1 cbug write ( 3, 9904) n, xep(n) cbug do 170 n = 1, nep cbug pcumre = n * pcum / nep cbug write ( 3, 9905) n, xep(n+1), pcumre cbug 170 continue cbugc***DEBUG ends. 210 return c.... End of subroutine apteqxn. (+1 line.) end UCRL-WEB-209832