subroutine aptlocs (rin, rout, ax, ay, az, np, tol, & bx, by, bz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTLOCS c c call aptlocs (rin, rout, ax, ay, az, np, tol, c & bx, by, bz, nerr) c c Version: aptlocs Updated 1990 November 29 11:40. c aptlocs Originated 1990 May 4 10:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np points b = (bx, by, bz) from a uniform c spatial distribution inside the spherical shell with inner c radius rin, outer radius rout, centered at point c a = (ax, ay, az). c Coordinates of "b" 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 History: 1990 May 7 16:00. Deleted arguments tolx, toly, tolz from call c to aptscat. c c Input: rin, rout, ax, ay, az, np, tol. c c Output: bx, by, bz, nerr. c c Calls: aptscat, aptspod c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a" at the center of c the spherical shell of inner radius rin, outer c radius rout. c c bx,by,bz Output The coordinates of spatial points "b" sampled randomly c from a uniform distribution over the volume of the c spherical shell. Size np. Values less than the c estimated error in their calculation, based on tol, c will be truncated to zero. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c 2 if rin is negative, or rout is less than rin. c c np Input Size of arrays bx, by, bz. c Number of points "b" to sample. c c rin Input The inner radius of the spherical shell. c c rout Input The outer radius of the spherical shell. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Coordinate x of point "b". dimension bx (1) c---- Coordinate y of point "b". dimension by (1) c---- Coordinate z of point "b". dimension bz (1) c.... Local variables. c---- Estimated error in bx. common /laptlocs/ bxerr c---- Estimated error in by. common /laptlocs/ byerr c---- Estimated error in bz. common /laptlocs/ bzerr c---- Component x of vector "c". common /laptlocs/ cx (64) c---- Component y of vector "c". common /laptlocs/ cy (64) c---- Component z of vector "c". common /laptlocs/ cz (64) c---- Index in local array. common /laptlocs/ n c---- First index of subset of data. common /laptlocs/ n1 c---- Last index of subset of data. common /laptlocs/ n2 c---- Index in external array. common /laptlocs/ nn c---- Size of current subset of data. common /laptlocs/ ns c---- Randomly sampled radius. common /laptlocs/ rsam (64) cbugc***DEBUG begins. cbug common /laptlocs/ avgx, avgy, avgz, devx, devy, devz cbug common /laptlocs/ avgx2, avgy2, avgz2 cbug common /laptlocs/ avgr, avgr2, devr, rad, rad2 cbug 9901 format (/ 'aptlocs sampling inside a spherical shell.', cbug & ' np=',i3 / cbug & ' rin,rout=',1p2e22.14 / cbug & ' ax,ay,az=',1p3e22.14) cbug write (3, 9901) np, rin, rout, ax, ay, az cbugc***DEBUG ends. c.... initialize. nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif if ((rout .lt. rin) .or. (rin .lt. 0.0)) then nerr = 2 go to 210 endif c.... Set up the indices of the first subset of data. n1 = 1 n2 = min (np, 64) c.... Loop over the data subsets. 110 ns = n2 - n1 + 1 c.... Randomly sample unit vectors "c" in the radial direction of the sphere. call aptscat (ns, cx, cy, cz, nerr) c.... Randomly sample radii from a quadratic distribution from rin to rout. call aptspod (2.0, rin, rout, ns, rsam, nerr) c.... Find points "b" sampled uniformly over the spherical shell. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 bx(nn) = ax + rsam(n) * cx(n) by(nn) = ay + rsam(n) * cy(n) bz(nn) = az + rsam(n) * cz(n) c---- End of loop over subset of data. 120 continue c.... See if truncation to zero is allowed. c---- Truncation to zero allowed. if (tol .gt. 0.0) then c---- Loop over subset of data. do 130 n = 1, ns nn = n + n1 - 1 bxerr = tol * (abs (ax) + rsam(n) * abs (cx(n))) byerr = tol * (abs (ay) + rsam(n) * abs (cy(n))) bzerr = tol * (abs (az) + rsam(n) * abs (cz(n))) if (abs(bx(nn)) .lt. bxerr) then bx(nn) = 0.0 endif if (abs(by(nn)) .lt. byerr) then by(nn) = 0.0 endif if (abs(bz(nn)) .lt. bzerr) then bz(nn) = 0.0 endif c---- End of loop over subset of data. 130 continue 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 cbug 9902 format (/ 'aptlocs results:') cbug 9903 format (i3,' bx,by,bz=',1p3e22.14 / cbug & ' rad= ',1pe22.14) cbug 9904 format (/ 'aptlocs mean and deviation:' / cbug & ' avg(x,y,z)=',1p3e22.14 / cbug & ' dev(x,y,z)=',1p3e22.14 / cbug & ' avgr,devr= ',1p2e22.14) cbug cbug write ( 3, 9902) cbug cbug avgr = 0.0 cbug avgx = 0.0 cbug avgy = 0.0 cbug avgz = 0.0 cbug avgr2 = 0.0 cbug avgx2 = 0.0 cbug avgy2 = 0.0 cbug avgz2 = 0.0 cbug cbug do 180 n = 1, np cbug rad2 = (bx(n) - ax)**2 + (by(n) - ay)**2 + cbug & (bz(n) - az)**2 cbug rad = sqrt (rad2) cbug avgr = avgr + rad cbug avgx = avgx + bx(n) cbug avgy = avgy + by(n) cbug avgz = avgz + bz(n) cbug avgr2 = avgr2 + rad2 cbug avgx2 = avgx2 + bx(n)**2 cbug avgy2 = avgy2 + by(n)**2 cbug avgz2 = avgz2 + bz(n)**2 cbug write ( 3, 9903) n, bx(n), by(n), bz(n), rad cbug 180 continue cbug cbug avgr = avgr / np cbug avgx = avgx / np cbug avgy = avgy / np cbug avgz = avgz / np cbug avgr2 = avgr2 / np cbug avgx2 = avgx2 / np cbug avgy2 = avgy2 / np cbug avgz2 = avgz2 / np cbug devr = sqrt (avgr2 - avgr**2) cbug devx = sqrt (avgx2 - avgx**2) cbug devy = sqrt (avgy2 - avgy**2) cbug devz = sqrt (avgz2 - avgz**2) cbug cbug write ( 3, 9904) avgx, avgy, avgz, devx, devy, devz, avgr, devr cbugc***DEBUG ends. 210 return c.... End of subroutine aptlocs. (+1 line.) end UCRL-WEB-209832