subroutine aptlocy (rin, rout, ax, ay, az, bx, by, bz, np, tol,
& cx, cy, cz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTLOCY
c
c call aptlocy (rin, rout, ax, ay, az, bx, by, bz, np, tol,
c & cx, cy, cz, nerr)
c
c Version: aptlocy Updated 1990 November 29 11:40.
c aptlocy Originated 1990 May 3 11:10.
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 inside the right circular cylindrical
c shell with inner radius rin, outer radius rout, with ends
c centered at points a = (ax, ay, az) and b = (bx, by, bz).
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 Input: rin, rout, ax, ay, az, bx, by, bz, np, tol.
c
c Output: cx, cy, cz, nerr.
c
c Calls: aptscad, aptslid, aptvdis
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a" at the center of
c the circular annulus of inner radius rin, outer
c radius rout, bounding one end of the cylinder.
c
c bx,by,bz Input The x, y, z coordinates of point "b" at the center of
c the circular annulus of inner radius rin, outer
c radius rout, bounding the other end of the cylinder.
c
c cx,cy,cz Output The coordinates of spatial points "c" sampled randomly
c from a uniform distribution over the volume of the
c cylindrical 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 3 if the cylinder is too short (point "a" is too
c close to point "b").
c
c np Input Size of arrays cx, cy, cz.
c Number of points "c" to sample.
c
c rin Input The inner radius of the right circular cylindrical
c shell.
c
c rout Input The outer radius of the right circular cylindrical
c 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 "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 /laptlocy/ abx
c---- Component y of vector "ab".
common /laptlocy/ aby
c---- Component z of vector "ab".
common /laptlocy/ abz
c---- Estimated error in cx.
common /laptlocy/ cxerr
c---- Estimated error in cy.
common /laptlocy/ cyerr
c---- Estimated error in cz.
common /laptlocy/ czerr
c---- Distance between points "a" and "b".
common /laptlocy/ dab
c---- Component x of vector "d".
common /laptlocy/ dx (64)
c---- Component y of vector "d".
common /laptlocy/ dy (64)
c---- Component z of vector "d".
common /laptlocy/ dz (64)
c---- A very small number.
common /laptlocy/ fuz
c---- Index in local array.
common /laptlocy/ n
c---- First index of subset of data.
common /laptlocy/ n1
c---- Last index of subset of data.
common /laptlocy/ n2
c---- Index in external array.
common /laptlocy/ nn
c---- Size of current subset of data.
common /laptlocy/ ns
c---- Relative probability of rin.
common /laptlocy/ pin
c---- Relative probability of rout.
common /laptlocy/ pout
c---- Random number (0.0 to 1.0).
common /laptlocy/ ranfp (64)
c---- Randomly sampled radius.
common /laptlocy/ rsam (64)
cbugc***DEBUG begins.
cbug common /laptlocy/ avgx, avgy, avgz, devx, devy, devz
cbug common /laptlocy/ avgx2, avgy2, avgz2
cbug 9901 format (/ 'aptlocy sampling inside a cylinder. np=',i3 /
cbug & ' rin,rout=',1p2e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14)
cbug write (3, 9901) np, rin, rout, 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 ((rout .lt. rin) .or. (rin .lt. 0.0)) then
nerr = 2
go to 210
endif
c.... Find vector "ab" along the axis of the cylinder from "a" to "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
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 "d" in the radial direction of the cylinder.
call aptscad (ns, abx, aby, abz, dx, dy, dz, nerr)
c.... Randomly sample radii from a linear distribution from rin to rout.
pin = rin + fuz
pout = rout + fuz
call aptslid (rin, pin, rout, pout, ns, rsam, nerr)
c.... Generate the needed 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.... Find points "c" sampled uniformly over the cylindrical shell.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
cx(nn) = ax + ranfp(n) * abx + rsam(n) * dx(n)
cy(nn) = ay + ranfp(n) * aby + rsam(n) * dy(n)
cz(nn) = az + ranfp(n) * abz + rsam(n) * dz(n)
c---- End of loop over subset of data.
130 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 140 n = 1, ns
nn = n + n1 - 1
cxerr = tol * (abs (ax) + ranfp(n) * abs (abx) +
& rsam(n) * abs (dx(n)))
cyerr = tol * (abs (ay) + ranfp(n) * abs (aby) +
& rsam(n) * abs (dy(n)))
czerr = tol * (abs (az) + ranfp(n) * abs (abz) +
& rsam(n) * 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.
140 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 (/ 'aptlocy results:')
cbug 9903 format (i3,' cx,cy,cz=',1p3e22.14)
cbug 9904 format (/ 'aptlocy mean and deviation:' /
cbug & ' avg(x,y,z)=',1p3e22.14 /
cbug & ' dev(x,y,z)=',1p3e22.14)
cbug
cbug write ( 3, 9902)
cbug
cbug avgx = 0.0
cbug avgy = 0.0
cbug avgz = 0.0
cbug avgx2 = 0.0
cbug avgy2 = 0.0
cbug avgz2 = 0.0
cbug
cbug do 180 n = 1, np
cbug avgx = avgx + cx(n)
cbug avgy = avgy + cy(n)
cbug avgz = avgz + cz(n)
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 180 continue
cbug
cbug avgx = avgx / np
cbug avgy = avgy / np
cbug avgz = avgz / np
cbug avgx2 = avgx2 / np
cbug avgy2 = avgy2 / np
cbug avgz2 = avgz2 / np
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
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptlocy. (+1 line.)
end
UCRL-WEB-209832