subroutine aptbrks (ksys, iunit, u, du, v, dv, w, dw, np, tol, & pu, pv, pw, fu, fv, fw, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTBRKS c c call aptbrks (ksys, iunit, u, du, v, dv, w, dw, np, tol, c & pu, pv, pw, fu, fv, fw, nerr) c c Version: aptbrks Updated 1994 March 3 12:00. c aptbrks Originated 1994 March 3 12:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np points p = (pu, pv, pw) from a uniform c spatial distribution in the the inside of the volume element c defined in the coordinate system ksys (0 = rectangular, c 1 = cylindrical, 2 = spherical) with angle units iunit c (0 = degrees, 1 = radians), bounded by the coordinate surfaces c at u, u + du, v, v + dv, w and w + dw. The fractional distances c across the volume element, fu, fv and fw will also be returned. c Coordinates of points "p" less than the estimated error in their c calculation, based on tol, will be truncated to zero. c Flag nerr indicates any input error. c c Input: ksys, iunit, u, du, v, dv, w, dw, np, tol. c c Output: pu, pv, pw, fu, fv, fw, nerr. c c Calls: None c c Glossary: c c fu,fv,fw Output The fractional distance (on a volume basis) of a c randomly sampled point "p" across the volume element. c Size np. c c iunit Input Indicates the unit used for angles: c 0 if angles theta and phi are in degrees. c 1 if angles theta and phi are in radians. c Not used if ksys = 0. c c ksys Input Indicates the coordinate system type: c 0 for Cartesian coordinates. u = x, v = y, w = z. c 1 for cylindrical coordinates. u = radius from z c axis, v = angle theta in xy plane, counterclockwise c from x axis, w = z. c 2 for spherical coordinates. u = radius from origin, c v = angle theta in xy plane, counterclockwise from c x axis, w = angle phi from z axis. c c nerr Output Indicates an input error, if not 0. c 1 if ksys is not 0, 1 or 2. c 2 if iunit is not 0 or 1. c 3 if np is not positive. c c np Input Size of arrays pu, pv, pw, fu, fv, fw. c Number of points "p" to sample. Must be positive. c c pu,pv,pw Output The coordinates of a spatial point "p" sampled randomly c from a uniform distribution over the volume element c defined in coordinate system ksys, with angle units c iunit, bounded by the coordinate surfaces at u = u, c u = u + du, v = v, v = bvr + dv, w = w and c w = dw. Values less than the estimated error in c their calculation, based on tol, will be truncated to c zero. Size np. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c c u, du Input The coordinate and increment of the bounding coordinate c surfaces of the volume element in the u direction. c c v, dv Input The coordinate and increment of the bounding coordinate c surfaces of the volume element in the v direction. c c w, dw Input The coordinate and increment of the bounding coordinate c surfaces of the volume element in the w direction. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c____ Mathematical constant pi. Dimensionless. parameter (pi = 3.14159265358979323) c.... Dimensioned arguments. c---- Coordinate u of point "p". dimension pu (1) c---- Coordinate v of point "p". dimension pv (1) c---- Coordinate w of point "p". dimension pw (1) c---- Fractional u distance of point "p". dimension fu (1) c---- Fractional v distance of point "p". dimension fv (1) c---- Fractional w distance of point "p". dimension fw (1) c.... Local variables. c---- Cosine of angle phi. common /laptbrks/ cosw c---- Cosine of minimum angle phi. common /laptbrks/ cosw1 c---- Cosine of maximum angle phi. common /laptbrks/ cosw2 c---- Angle conversion factor. common /laptbrks/ fact c---- Index in output array. common /laptbrks/ n c---- Estimated error in pu. common /laptbrks/ puerr c---- Estimated error in pv. common /laptbrks/ pverr c---- Estimated error in pw. common /laptbrks/ pwerr c---- Square of minimum radius. common /laptbrks/ r1sq c---- Square of maximum radius. common /laptbrks/ r2sq cbugc***DEBUG begins. cbug common /laptbrks/ avgu, avgv, avgw, devu, devv, devw cbug common /laptbrks/ avgu2, avgv2, avgw2 cbug 9901 format (/ 'aptbrks sampling in a brick. np=',i3,' ksys=',i3, cbug & ' iunit=',i2 / cbug & ' tol=',1pe22.14 / cbug & ' u, du= ',1p2e22.14 / cbug & ' v, dv= ',1p2e22.14 / cbug & ' w, dw= ',1p2e22.14 ) cbug write ( 3, 9901) np, ksys, iunit, tol, u, du, v, dv, w, dw cbugc***DEBUG ends. c.... initialize. nerr = 0 c.... Test for input errors. if ((ksys .lt. 0) .or. (ksys .gt. 2)) then nerr = 1 go to 210 endif if ((iunit .ne. 0) .and. (iunit .ne. 1)) then nerr = 2 go to 210 endif if (np .le. 0) then nerr = 3 go to 210 endif c.... Find the angle conversion factor (degrees to radians). if (iunit .eq. 0) then fact = pi / 180.0 else fact = 1.0 endif c.... Find the needed random numbers. c---- Loop over data. do 120 n = 1, np fu(n) = ranf( ) c---- End of loop over data. 120 continue c---- Loop over data. do 130 n = 1, np fv(n) = ranf( ) c---- End of loop over data. 130 continue c---- Loop over data. do 140 n = 1, np fw(n) = ranf( ) c---- End of loop over data. 140 continue c.... Find points "p" in the brick. if (ksys .eq. 0) then c---- Loop over data. do 150 n = 1, np pu(n) = u + fu(n) * du pv(n) = v + fv(n) * dv pw(n) = w + fw(n) * dw c---- End of loop over data. 150 continue elseif (ksys .eq. 1) then c---- Loop over data. do 160 n = 1, np r1sq = u**2 r2sq = (u + du)**2 pu(n) = sqrt (r1sq + fu(n) * (r2sq - r1sq)) pv(n) = v + fv(n) * dv pw(n) = w + fw(n) * dw c---- End of loop over data. 160 continue elseif (ksys .eq. 2) then c---- Loop over data. do 170 n = 1, np r1cb = u**3 r2cb = (u + du)**3 pu(n) = (r1cb + fu(n) * (r2cb - r1cb))**(1.0 / 3.0) pv(n) = v + fv(n) * dv cosw1 = cos (fact * w) cosw2 = cos (fact * (w + dw)) cosw = cosw1 + fw(n) * (cosw2 - cosw1) pw(n) = acos (cosw) / fact c---- End of loop over data. 170 continue c---- Tested ksys. endif c.... See if truncation to zero is allowed. c---- Truncation to zero allowed. if (tol .gt. 0.0) then if (ksys .eq. 0) then c---- Loop over data. do 180 n = 1, np puerr = tol * (abs (u) + fu(n) * abs (du)) pverr = tol * (abs (v) + fv(n) * abs (dv)) pwerr = tol * (abs (w) + fw(n) * abs (dw)) if (abs (pu(n)) .lt. puerr) then pu(n) = 0.0 endif if (abs (pv(n)) .lt. pverr) then pv(n) = 0.0 endif if (abs (pw(n)) .lt. pwerr) then pw(n) = 0.0 endif c---- End of loop over data. 180 continue elseif (ksys .eq. 1) then c---- Loop over data. do 190 n = 1, np pverr = tol * (abs (v) + fv(n) * abs (dv)) pwerr = tol * (abs (w) + fw(n) * abs (dw)) if (abs (pv(n)) .lt. pverr) then pv(n) = 0.0 endif if (abs (pw(n)) .lt. pwerr) then pw(n) = 0.0 endif c---- End of loop over data. 190 continue elseif (ksys .eq. 2) then c---- Loop over data. do 200 n = 1, np cosw1 = cos (fact * w) cosw2 = cos (fact * (w + dw)) cosw = cosw1 + fw(n) * (cosw2 - cosw1) pverr = tol * (abs (v) + fv(n) * abs (dv)) cwerr = tol * (abs (cosw1) + & fw(n) * (abs (cosw2) + abs (cosw1))) if (abs (pv(n)) .lt. pverr) then pv(n) = 0.0 endif if (abs (cosw - 1.0) .lt. cwerr) then pw(n) = 0.0 endif if (abs (cosw) .lt. cwerr) then pw(n) = 0.5 * pi / fact endif if (abs (cosw + 1.0) .lt. cwerr) then pw(n) = pi / fact endif c---- End of loop over data. 200 continue c---- Tested ksys. endif c---- Tested tol. endif cbugc***DEBUG begins. cbug cbug 9902 format (/ 'aptbrks results:') cbug 9903 format (i3,' pu,pv,pw=',1p3e22.14) cbug 9904 format (i3,' fu,fv,fw=',1p3e22.14 ) cbug 9905 format (/ 'aptbrks mean and deviation:' / cbug & ' avg(u,v,w)= ',1p3e22.14 / cbug & ' dev(u,v,w)= ',1p3e22.14) cbug cbug write ( 3, 9902) cbug cbug avgu = 0.0 cbug avgv = 0.0 cbug avgw = 0.0 cbug avgu2 = 0.0 cbug avgv2 = 0.0 cbug avgw2 = 0.0 cbug cbug do 205 n = 1, np cbug avgu = avgu + pu(n) cbug avgv = avgv + pv(n) cbug avgw = avgw + pw(n) cbug avgu2 = avgu2 + pu(n)**2 cbug avgv2 = avgv2 + pv(n)**2 cbug avgw2 = avgw2 + pw(n)**2 cbug write ( 3, 9903) n, pu(n), pv(n), pw(n) cbug write ( 3, 9904) n, fu(n), fv(n), fw(n) cbug 205 continue cbug cbug avgu = avgu / np cbug avgv = avgv / np cbug avgw = avgw / np cbug avgu2 = avgu2 / np cbug avgv2 = avgv2 / np cbug avgw2 = avgw2 / np cbug devu = sqrt (amax1 (0.0, avgu2 - avgu**2)) cbug devv = sqrt (amax1 (0.0, avgv2 - avgv**2)) cbug devw = sqrt (amax1 (0.0, avgw2 - avgw**2)) cbug cbug write ( 3, 9905) avgu, avgv, avgw, devu, devv, devw cbugc***DEBUG ends. 210 return c.... End of subroutine aptbrks. (+1 line.) end UCRL-WEB-209832