subroutine aptspar (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, np, tol, px, py, pz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPAR
c
c call aptspar (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, px, py, pz, nerr)
c
c Version: aptspar Updated 1991 June 10 17:30.
c aptspar Originated 1991 June 10 17:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np points p = (px, py, pz) from a uniform
c spatial distribution in the the inside of the parallelopiped
c with vertices at points a = (ax, ay, az), b = (bx, by, bz),
c c = (cx, cy, cz) and d = (dx, dy, dz), with edges "ab", "ac",
c and "ad" meeting at point "a"; or on the planar surface of a
c parallelogram if one of the three points "b", "c" or "d"
c coincides with point "a"; or on a line if two of the three
c points "b", "c" and "d" coincide with point "a"; or the point
c "a" if all four points coincide.
c Coordinates of points "p" 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: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, tol.
c
c Output: px, py, pz, nerr.
c
c Calls: aptvdis
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a".
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" on edge "ab".
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" on edge "ac".
c
c dx,dy,dz Input The x, y, z coordinates of vertex "d" on edge "ad".
c
c px,py,pz Output The coordinates of spatial points "p" sampled randomly
c from a uniform distribution over the volume of the
c parallelopiped with edges "ab", "ac" and "ad" meeting
c at vertex "a".
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
c np Input Size of arrays px, py, pz.
c Number of points "p" to sample.
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 "p".
dimension px (1)
c---- Coordinate y of point "p".
dimension py (1)
c---- Coordinate z of point "p".
dimension pz (1)
c.... Local variables.
c---- Component x of vector "ab".
common /laptspar/ abx
c---- Component y of vector "ab".
common /laptspar/ aby
c---- Component z of vector "ab".
common /laptspar/ abz
c---- Component x of vector "ac".
common /laptspar/ acx
c---- Component y of vector "ac".
common /laptspar/ acy
c---- Component z of vector "ac".
common /laptspar/ acz
c---- Component x of vector "ad".
common /laptspar/ adx
c---- Component y of vector "ad".
common /laptspar/ ady
c---- Component z of vector "ad".
common /laptspar/ adz
c---- Index in local array.
common /laptspar/ n
c---- First index of subset of data.
common /laptspar/ n1
c---- Last index of subset of data.
common /laptspar/ n2
c---- Index in external array.
common /laptspar/ nn
c---- Size of current subset of data.
common /laptspar/ ns
c---- Random number from 0.0 to 1.0.
common /laptspar/ ranfp1 (64)
c---- Random number from 0.0 to 1.0.
common /laptspar/ ranfp2 (64)
c---- Random number from 0.0 to 1.0.
common /laptspar/ ranfp3 (64)
c---- Estimated error in px.
common /laptspar/ pxerr
c---- Estimated error in py.
common /laptspar/ pyerr
c---- Estimated error in pz.
common /laptspar/ pzerr
c---- Distance between points "a" and "b".
common /laptspar/ vlenab
c---- Distance between points "a" and "c".
common /laptspar/ vlenac
c---- Distance between points "a" and "d".
common /laptspar/ vlenad
cbugc***DEBUG begins.
cbug common /laptspar/ avgx, avgy, avgz, devx, devy, devz
cbug common /laptspar/ avgx2, avgy2, avgz2
cbug 9901 format (/ 'aptspar sampling in a parallelopiped. np=',i3 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' dx,dy,dz=',1p3e22.14)
cbug write (3, 9901) np, ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Find the vectors and distances along the edges "ab", "ac" and "ad".
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, vlenab, nerr)
call aptvdis (ax, ay, az, cx, cy, cz, 1, tol,
& acx, acy, acz, vlenac, nerr)
call aptvdis (ax, ay, az, dx, dy, dz, 1, tol,
& adx, ady, adz, vlenad, nerr)
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 needed random numbers.
c---- Loop over subset of data.
do 120 n = 1, ns
ranfp1(n) = ranf( )
c---- End of loop over subset of data.
120 continue
c---- Loop over subset of data.
do 130 n = 1, ns
ranfp2(n) = ranf( )
c---- End of loop over subset of data.
130 continue
c---- Loop over subset of data.
do 140 n = 1, ns
ranfp3(n) = ranf( )
c---- End of loop over subset of data.
140 continue
c.... Find points "p" in the parallelopiped.
c---- Loop over subset of data.
do 150 n = 1, ns
nn = n + n1 - 1
px(nn) = ax + ranfp1(n) * abx + ranfp2(n) * acx +
& ranfp3(n) * adx
py(nn) = ay + ranfp1(n) * aby + ranfp2(n) * acy +
& ranfp3(n) * ady
pz(nn) = az + ranfp1(n) * abz + ranfp2(n) * acz +
& ranfp3(n) * adz
c---- End of loop over subset of data.
150 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
pxerr = tol * (ranfp1(n) * abs (abx) +
& ranfp2(n) * abs (acx) +
& ranfp3(n) * abs (adx) + abs (ax))
pyerr = tol * (ranfp1(n) * abs (aby) +
& ranfp2(n) * abs (acy) +
& ranfp3(n) * abs (ady) + abs (ay))
pzerr = tol * (ranfp1(n) * abs (abz) +
& ranfp2(n) * abs (acz) +
& ranfp3(n) * abs (adz) + abs (az))
if (abs(px(nn)) .lt. pxerr) then
px(nn) = 0.0
endif
if (abs(py(nn)) .lt. pyerr) then
py(nn) = 0.0
endif
if (abs(pz(nn)) .lt. pzerr) then
pz(nn) = 0.0
endif
c---- End of loop over subset of data.
180 continue
c---- Tested tol.
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 (/ 'aptspar results:')
cbug 9903 format (i4,' px,py,pz=',1p3e22.14)
cbug 9904 format (/ 'aptspar 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 190 n = 1, np
cbug avgx = avgx + px(n)
cbug avgy = avgy + py(n)
cbug avgz = avgz + pz(n)
cbug avgx2 = avgx2 + px(n)**2
cbug avgy2 = avgy2 + py(n)**2
cbug avgz2 = avgz2 + pz(n)**2
cbug write ( 3, 9903) n, px(n), py(n), pz(n)
cbug 190 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 (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
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptspar. (+1 line.)
end
UCRL-WEB-209832