subroutine aptscnv (frad, ra, ax, ay, az, rb, bx, by, bz, & np, tol, cx, cy, cz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTSCNV c c call aptscnv (frad, ra, ax, ay, az, rb, bx, by, bz, np, tol, c & cx, cy, cz, nerr) c c Version: aptscnv Updated 1991 June 6 13:40. c aptscnv Originated 1991 June 6 13:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np points c = (cx, cy, cz) from a uniform c spatial distribution in the volume bounded by the surfaces of c two right circular cones with the same vertex and the same axis, c and by two parallel planes perpendicular to that axis. The two c planar boundaries are the annular disk centered at point c a = (ax, ay, az) with inner and outer radii frad * ra and ra, c and the annular disk centered at point b = (bx, by, bz) with c inner and outer radii frad * rb and rb. If frad = 1.0, the c points "c" will be sampled uniformly over the conical surface c between the two circles. c Coordinates of "c" 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 Notes: The sampled volume is equal to c (1.0 - frad**2) * pi * (ra**2 + rb**2 + ra * rb) * h / 3, c and the area of the outer conical surface between the disks is c pi * (ra + rb) * sqrt ((rb - ra)**2 + h**2), c where pi = 3.14159265358979..., and h is the distance between c points "a" and "b", equal to c sqrt ((bx - ax)**2 + (by - ay)**2 + (bz - az)**2). c The vertex of the cones is at point v = (vx, vy, vz), given by c v = (rb * a - ra * b) / (rb - ra), c where v, a and b are vectors. c The vertex half-angles of the inner and outer conical surfaces c are equal to c arctan (frad * (rb - ra) / h) and c arctan ((rb - ra) / h). c c Input: frad, ra, ax, ay, az, rb, bx, by, bz, np, tol. c c Output: cx, cy, cz, nerr. c c Calls: aptlocd, aptscad, aptslid, aptspod, aptvdis c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a" on the axis of c the two conical surfaces. c c bx,by,bz Input The x, y, z coordinates of point "b" on the axis of c the two conical surfaces. c c cx,cy,cz Output The coordinates of spatial points "c" sampled randomly c from a uniform distribution over the volume bounded c by two conical surfaces and two planes. Size np. c Values less than the estimated error in their c calculation, based on tol, will be truncated to zero. c c frad Input The ratio of the inner radius to the outer radius of c each of the two annular disks centered at points c "a" and "b". 0.0 .le. frad .le. 1.0. c Use 0.0 to randomly sample the entire volume, and c use 1.0 to randomly sample only on the outer surface. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c 2 if ra or rb is negative. c 3 if the cone is too short (point "a" is too c close to point "b"). c 4 if frad .lt. 0.0 or frad .gt. 1.0. c c np Input Size of arrays cx, cy, cz. c Number of points "c" to sample. c c ra Input The outer radius of the annular disk centered at point c "a", and perpendicular to the axis "ab". c The inner radius is frad * ra. c c rb Input The outer radius of the annular disk centered at point c "b", and perpendicular to the axis "ab". c The inner radius is frad * rb. 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 "c". dimension cx (1) c---- Coordinate y of point "c". dimension cy (1) c---- Coordinate z of point "c". dimension cz (1) c.... Local variables. c---- Component x of vector "ab". common /laptscnv/ abx c---- Component y of vector "ab". common /laptscnv/ aby c---- Component z of vector "ab". common /laptscnv/ abz c---- Estimated error in cx. common /laptscnv/ cxerr c---- Estimated error in cy. common /laptscnv/ cyerr c---- Estimated error in cz. common /laptscnv/ czerr c---- Distance between points "a" and "b". common /laptscnv/ dab c---- Difference 1.0 - frad. common /laptscnv/ df c---- Estimated error in df. common /laptscnv/ dferr c---- Difference rb - ra. common /laptscnv/ dr c---- Estimated error in dr. common /laptscnv/ drerr c---- Component x of vector "d". common /laptscnv/ dx (64) c---- Component y of vector "d". common /laptscnv/ dy (64) c---- Component z of vector "d". common /laptscnv/ dz (64) c---- Randomly sampled fractional distance. common /laptscnv/ fab (64) c---- A very small number. common /laptscnv/ fuz c---- Index in local array. common /laptscnv/ n c---- First index of subset of data. common /laptscnv/ n1 c---- Last index of subset of data. common /laptscnv/ n2 c---- Index in external array. common /laptscnv/ nn c---- Size of current subset of data. common /laptscnv/ ns c---- Relative probability of ra. common /laptscnv/ pa c---- Relative probability of rb. common /laptscnv/ pb c---- Radius of cone at point c. common /laptscnv/ rc (64) cbugc***DEBUG begins. cbug common /laptscnv/ avgr, avgx, avgy, avgz, devr, devx, devy, devz cbug common /laptscnv/ avgr2, avgx2, avgy2, avgz2 cbug 9901 format (/ 'aptscnv sampling in a conical annulus. np=',i3 / cbug & ' frad =',1pe22.14 / cbug & ' ra,rb =',1p2e22.14 / cbug & ' ax,ay,az=',1p3e22.14 / cbug & ' bx,by,bz=',1p3e22.14) cbug write (3, 9901) np, frad, ra, rb, ax, ay, az, bx, by, bz cbugc***DEBUG ends. c.... initialize. c---- A very small number. fuz = 1.e-99 nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif if ((ra .lt. 0.0) .or. (rb .lt. 0.0)) then nerr = 2 go to 210 endif c.... Find vector "ab" between points "a" and "b". call aptvdis (ax, ay, az, bx, by, bz, 1, tol, & abx, aby, abz, dab, nerr) if (dab .le. tol) then nerr = 3 go to 210 endif if ((frad .lt. 0.0) .or. (frad .gt. 1.0)) then nerr = 4 go to 210 endif c.... Find the nature of the sampling space. dr = rb - ra drerr = tol * (ra + rb) c---- Shape is a cylinder or a line. if (abs (dr) .le. drerr) then dr = 0.0 endif df = 1.0 - frad dferr = tol * frad c---- Shape is a surface or a line. if (abs (df) .le. dferr) then df = 0.0 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.... Find the fractional distances of points "c" along the axis "ab". c---- Shape is a cylinder or a line. if (dr .eq. 0.0) then c---- Loop over subset of data. do 120 n = 1, ns c---- Uniform distribution. fab(n) = ranf( ) c---- End of loop over subset of data. 120 continue c---- Loop over subset of data. do 130 n = 1, ns rc(n) = ra + fab(n) * dr c---- End of loop over subset of data. 130 continue c---- Sample on a surface. elseif (df .eq. 0.0) then pa = ra + fuz pb = rb + fuz c++++ Linear distribution. call aptslid (0.0, pa, 1.0, pb, ns, fab, nerr) c---- Loop over subset of data. do 140 n = 1, ns rc(n) = frad * (ra + fab(n) * dr) c---- End of loop over subset of data. 140 continue c---- Sample in volume, ra > rb. elseif (dr .lt. 0.0) then c++++ Quadratic distribution. call aptspod (2.0, rb, ra, ns, rc, nerr) c---- Loop over subset of data. do 150 n = 1, ns fab(n) = (rc(n) - ra) / dr c---- End of loop over subset of data. 150 continue c---- Sample in volume, rb > ra. else c++++ Quadratic distribution. call aptspod (2.0, ra, rb, ns, rc, nerr) c---- Loop over subset of data. do 160 n = 1, ns fab(n) = (rc(n) - ra) / dr c---- End of loop over subset of data. 160 continue c---- Tested dr for shape of cone. endif c.... Randomly sample vectors perpendicular to the axis of the cone. c---- On outer surface. if (df .eq. 0.0) then call aptscad (ns, abx, aby, abz, dx, dy, dz, nerr) c---- On annular disk. else call aptlocd (frad, 1.0, 0., 0., 0., abx, aby, abz, ns, tol, & dx, dy, dz, nerr) endif c.... _ _ __ _ c.... Find points "c". c = a + fab * ab + rc * d. c---- Loop over subset of data. do 170 n = 1, ns nn = n + n1 - 1 cx(nn) = ax + fab(n) * abx + rc(n) * dx(n) cy(nn) = ay + fab(n) * aby + rc(n) * dy(n) cz(nn) = az + fab(n) * abz + rc(n) * dz(n) c---- End of loop over subset of data. 170 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 180 n = 1, ns nn = n + n1 - 1 cxerr = tol * (abs (ax) + fab(n) * (abs (ax) + abs (bx)) + & 2.0 * (ra + fab(n) * (ra + rb)) * abs (dx(n))) cyerr = tol * (abs (ay) + fab(n) * (abs (ay) + abs (by)) + & 2.0 * (ra + fab(n) * (ra + rb)) * abs (dy(n))) czerr = tol * (abs (az) + fab(n) * (abs (az) + abs (bz)) + & 2.0 * (ra + fab(n) * (ra + rb)) * abs (dz(n))) if (abs(cx(nn)) .lt. cxerr) then cx(nn) = 0.0 endif if (abs(cy(nn)) .lt. cyerr) then cy(nn) = 0.0 endif if (abs(cz(nn)) .lt. czerr) then cz(nn) = 0.0 endif c---- End of loop over subset of data. 180 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 (/ 'aptscnv results:') cbug 9903 format (i4,' cx,cy,cz=',1p3e22.14) cbug 9904 format (/ 'aptscnv 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 190 n = 1, np cbug avgr = avgr + sqrt (cx(n)**2 + cy(n)**2 + cz(n)**2) cbug avgx = avgx + cx(n) cbug avgy = avgy + cy(n) cbug avgz = avgz + cz(n) cbug avgr2 = avgr2 + cx(n)**2 + cy(n)**2 + cz(n)**2 cbug avgx2 = avgx2 + cx(n)**2 cbug avgy2 = avgy2 + cy(n)**2 cbug avgz2 = avgz2 + cz(n)**2 cbug write ( 3, 9903) n, cx(n), cy(n), cz(n) cbug 190 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 (amax1 (0.0, avgr2 - avgr**2)) cbug devx = sqrt (amax1 (0.0, avgx2 - avgx**2)) cbug devy = sqrt (amax1 (0.0, avgy2 - avgy**2)) cbug devz = sqrt (amax1 (0.0, 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 aptscnv. (+1 line.) end UCRL-WEB-209832