subroutine aptlocd (rin, rout, ax, ay, az, bx, by, bz, np, tol,
& cx, cy, cz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTLOCD
c
c call aptlocd (rin, rout, ax, ay, az, bx, by, bz, np, tol,
c & cx, cy, cz, nerr)
c
c Version: aptlocd Updated 1991 May 28 17:40
c aptlocd 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 on the circular annulus with inner radius
c rin, outer radius rout, centered at point a = (ax, ay, az), and
c perpendicular to the vector b = (bx, by, bz) in 3-D space.
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
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, perpendicular to the vector "b".
c
c bx,by,bz Input The x, y, z components of the vector "b" perpendicular
c the plane of the circular annulus centered at point
c "a".
c
c cx,cy,cz Output The coordinates of spatial points "c" sampled randomly
c from a uniform distribution over the surface of the
c annulus of inner radius rin, outer radius rout,
c centered at point "a", in the plane perpendicular to
c vector "b". Size np. Values less than the estimated
c error in their calculation, based on tol, will be
c 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 vector "b" is too short.
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 circular annulus.
c
c rout Input The outer radius of the circular annulus.
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---- Square of length of vector "b".
common /laptlocd/ blen2
c---- Estimated error in cx.
common /laptlocd/ cxerr
c---- Estimated error in cy.
common /laptlocd/ cyerr
c---- Estimated error in cz.
common /laptlocd/ czerr
c---- Component x of vector "d".
common /laptlocd/ dx (64)
c---- Component y of vector "d".
common /laptlocd/ dy (64)
c---- Component z of vector "d".
common /laptlocd/ dz (64)
c---- A very small number.
common /laptlocd/ fuz
c---- Index in local array.
common /laptlocd/ n
c---- First index of subset of data.
common /laptlocd/ n1
c---- Last index of subset of data.
common /laptlocd/ n2
c---- Index in external array.
common /laptlocd/ nn
c---- Size of current subset of data.
common /laptlocd/ ns
c---- Relative probability of rin.
common /laptlocd/ pin
c---- Relative probability of rout.
common /laptlocd/ pout
c---- Randomly sampled radius.
common /laptlocd/ rsam (64)
cbugc***DEBUG begins.
cbug common /laptlocd/ avgx, avgy, avgz, devx, devy, devz
cbug common /laptlocd/ avgx2, avgy2, avgz2
cbug 9901 format (/ 'aptlocd sampling on an annulus. 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
blen2 = bx**2 + by**2 + bz**2
if (blen2 .le. tol**2) 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 annulus.
call aptscad (ns, bx, by, bz, 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.... Find points "c" sampled uniformly over the surface of the annulus.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
cx(nn) = ax + rsam(n) * dx(n)
cy(nn) = ay + rsam(n) * dy(n)
cz(nn) = az + rsam(n) * dz(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
cxerr = tol * (abs (ax) + rsam(n) * abs (dx(n)))
cyerr = tol * (abs (ay) + rsam(n) * abs (dy(n)))
czerr = tol * (abs (az) + 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.
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 (/ 'aptlocd results:')
cbug 9903 format (i3,' cx,cy,cz=',1p3e22.14)
cbug 9904 format (/ 'aptlocd 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 aptlocd. (+1 line.)
end
UCRL-WEB-209832