subroutine aptscns (ra, ax, ay, az, rb, bx, by, bz, np, tol,
& cx, cy, cz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCNS
c
c call aptscns (ra, ax, ay, az, rb, bx, by, bz, np, tol,
c & cx, cy, cz, nerr)
c
c Version: aptscns Updated 1991 May 24 9:40.
c aptscns Originated 1991 May 24 9: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 on the curved surface of the frustrum of a
c right circular cone bounded by a disk of radius ra centered at
c point a = (ax, ay, az), and a disk of radius rb centered at
c point b = (bx, by, bz). The curved surface has an area equal
c to pi * (ra + rb) * sqrt ((rb - ra)**2 + h**2), where h is the
c distance from point "a" to point "b".
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: ra, ax, ay, az, rb, 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 disk of radius ra, bounding one end of
c the frustrum of a right circular cone.
c
c bx,by,bz Input The x, y, z coordinates of point "b" at the center of
c the circular disk of radius rb, bounding one end of
c the frustrum of a right circular cone.
c
c cx,cy,cz Output The coordinates of spatial points "c" sampled randomly
c from a uniform distribution over the curved surface
c of the frustrum of a right circular cone. Size np.
c Values less than the estimated error in their
c calculation, based on tol, 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 ra or rb is negative.
c 3 if the cone 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 ra Input The radius of the disk centered at point "a", forming
c one end of the frustrum of a right circular cone.
c
c rb Input The radius of the disk centered at point "b", forming
c one end of the frustrum of a right circular cone.
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 /laptscns/ abx
c---- Component y of vector "ab".
common /laptscns/ aby
c---- Component z of vector "ab".
common /laptscns/ abz
c---- Estimated error in cx.
common /laptscns/ cxerr
c---- Estimated error in cy.
common /laptscns/ cyerr
c---- Estimated error in cz.
common /laptscns/ czerr
c---- Distance between points "a" and "b".
common /laptscns/ dab
c---- Component x of vector "d".
common /laptscns/ dx (64)
c---- Component y of vector "d".
common /laptscns/ dy (64)
c---- Component z of vector "d".
common /laptscns/ dz (64)
c---- Randomly sampled fractional distance.
common /laptscns/ fab (64)
c---- A very small number.
common /laptscns/ fuz
c---- Index in local array.
common /laptscns/ n
c---- First index of subset of data.
common /laptscns/ n1
c---- Last index of subset of data.
common /laptscns/ n2
c---- Index in external array.
common /laptscns/ nn
c---- Size of current subset of data.
common /laptscns/ ns
c---- Relative probability of ra.
common /laptscns/ pa
c---- Relative probability of rb.
common /laptscns/ pb
c---- Distance of point "c" from axis.
common /laptscns/ rc (64)
cbugc***DEBUG begins.
cbug common /laptscns/ avgx, avgy, avgz, devx, devy, devz
cbug common /laptscns/ avgx2, avgy2, avgz2
cbug 9901 format (/ 'aptscns sampling on a conical frustrum. np=',i3 /
cbug & ' ra,rb =',1p2e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14)
cbug write (3, 9901) np, 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" along the axis of the cone 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 a fractional distance between points "a" and "b" from
c.... a linear distribution with weights ra and rb, respectively.
pa = ra + fuz
pb = rb + fuz
call aptslid (0.0, pa, 1.0, pb, ns, fab, nerr)
c.... Randomly sample unit vectors "d" in the radial direction of the cone.
call aptscad (ns, abx, aby, abz, dx, dy, dz, nerr)
c.... Find points "c" sampled uniformly over the conical surface.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
rc(n) = ra + fab(n) * (rb - ra)
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.
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) + fab(n) * (abs (ax) + abs (bx)) +
& (ra + fab(n) * (ra + rb)) * abs (dx(n)))
cyerr = tol * (abs (ay) + fab(n) * (abs (ay) + abs (by)) +
& (ra + fab(n) * (ra + rb)) * abs (dy(n)))
czerr = tol * (abs (az) + fab(n) * (abs (az) + abs (bz)) +
& (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.
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 (/ 'aptscns results:')
cbug 9903 format (i4,' cx,cy,cz=',1p3e22.14)
cbug 9904 format (/ 'aptscns 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 aptscns. (+1 line.)
end
UCRL-WEB-209832