subroutine apttlod (ax, ay, az, bx, by, bz, cx, cy, cz,
& np, kloc, tol, px, py, pz, wa, wb, wc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTTLOD
c
c call apttlod (ax, ay, az, bx, by, bz, cx, cy, cz,
c & np, kloc, tol, px, py, pz, wa, wb, wc, nerr)
c
c Version: apttlod Updated 1990 December 3 16:20.
c apttlod Originated 1990 February 8 14:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np spatial points p = (px, py, pz) from a
c uniform distribution over the triangle in 3-D space with
c vertices a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz),
c in any order around the triangle, and optionally (kloc = 1), to
c return the local coordinates wa, wb and wc for each point "p":
c
c vn = (b - a) X (c - b) = (c - b) X (a - c) = (a - c) X (b - a)
c
c wa = (vn . (c - b) X (p - c)) / vn**2
c wb = (vn . (a - c) X (p - a)) / vn**2
c wb = (vn . (b - a) X (p - b)) / vn**2
c
c p = wa * a + wb * b + wc * c
c
c where "X" indicates the vector product, and "." indicates the
c dot product.
c
c The relationship of the triangle local coordinates wa, wb and wc
c to the quadrangle local coordinates fdk and fdl is as follows:
c
c wa = 1 - fdk fdk = wb + wc
c wb = fdk * (1 - fdl) fdl = wc / (wb + wc)
c wc = fdk * fdl
c
c where fdk is the fractional distance from edge "aa" to edge
c "bc", and fdl is the fractional distance from edge "ab" to
c edge "ca", and any point "p" is given by:
c
c p = a + fdk * (b - a) + fdk * fdl * (c - b)
c
c Flag nerr indicates any input error.
c
c History: 1990 May 10 16:10. Replaced arguments fdk, fdl with wa, wb, wc,
c and made calculation and return of wa, wb and wc optional,
c based on new input argument kloc.
c 1990 May 11 10:00. Replaced apttlod arguments fdk, fdl with
c kloc, twa, twb, twc.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, np, kloc, tol.
c
c Output: px, py, pz, wa, wb, wc, nerr.
c
c Calls: aptslid
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a" of triangle.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of triangle.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of triangle.
c
c kloc Input Indicates local coordinates wa, wb and wc are to be
c returned, if 1. Otherwise, 0.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kloc is not 0 or 1.
c
c np Input Size of arrays wa, wb, wc, px, py, pz.
c Number of points "p" to sample.
c
c px,py,pz Output Sampled point "p". Size np. Coordinates may be
c truncated to zero, if less than the estimated error
c in their calculation, based on tol.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c wa Output Fractional distance of point "p" to vertex "a" from
c edge "bc". Size np.
c
c wb Output Fractional distance of point "p" to vertex "b" from
c edge "ca". Size np.
c
c wc Output Fractional distance of point "p" to vertex "c" from
c edge "ab". Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cooodinate x of a sampled point "p".
dimension px (1)
c---- Coordinate y of a sampled point "p".
dimension py (1)
c---- Coordinate z of a sampled point "p".
dimension pz (1)
c---- Fractional distance to "a" from "bc".
dimension wa (1)
c---- Fractional distance to "b" from "ca".
dimension wb (1)
c---- Fractional distance to "c" from "ab".
dimension wc (1)
c.... Local variables.
c---- Fractional distance from "a" to "bc".
common /lapttlod/ fdk (64)
c---- Fractional distance from "ab" to "ca".
common /lapttlod/ fdl (64)
c---- Index in local array.
common /lapttlod/ n
c---- First index of subset of data.
common /lapttlod/ n1
c---- Last index of subset of data.
common /lapttlod/ n2
c---- Index in external array.
common /lapttlod/ nn
c---- Size of current subset of data.
common /lapttlod/ ns
c---- Estimated error in px.
common /lapttlod/ pxerr
c---- Estimated error in py.
common /lapttlod/ pyerr
c---- Estimated error in pz.
common /lapttlod/ pzerr
cbugc***DEBUG begins.
cbug common /lapttlod/ avgx, avgy, avgz, devx, devy, devz
cbug common /lapttlod/ avgx2, avgy2, avgz2
cbug common /lapttlod/ avgwa, avgwb, avgwc
cbug common /lapttlod/ devwa, devwb, devwc
cbug common /lapttlod/ avgwa2, avgwb2, avgwc2
cbug 9901 format (/ 'apttlod sampling from a triangle. np=',i3 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14)
cbug write (3, 9901) np, ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((kloc .lt. 0) .or. (kloc .gt. 1)) then
nerr = 2
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.... Sample ns values of fdk from a linear distribution, with
c.... probability 0.0 at vertex "a", probability 2.0 at edge "bc".
call aptslid (0.0, 0.0, 1.0, 2.0, ns, fdk, nerr)
c.... Sample ns values of fdl from a uniform distribution, with
c.... equal probabilities at sides "ab" and "ca".
c---- Loop over subset of data.
do 120 n = 1, ns
fdl(n) = ranf( )
c---- End of loop over subset of data.
120 continue
c.... Find the coordinates of the points.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
px(nn) = ax + fdk(n) * (bx - ax) + fdk(n) * fdl(n) * (cx - bx)
py(nn) = ay + fdk(n) * (by - ay) + fdk(n) * fdl(n) * (cy - by)
pz(nn) = az + fdk(n) * (bz - az) + fdk(n) * fdl(n) * (cz - bz)
c---- End of loop over subset of data.
130 continue
c.... See if truncation 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
pxerr = tol * (fdk(n) * abs (bx - ax) +
& fdk(n) * fdl(n) * abs (cx - bx))
pyerr = tol * (fdk(n) * abs (by - ay) +
& fdk(n) * fdl(n) * abs (cy - by))
pzerr = tol * (fdk(n) * abs (bz - az) +
& fdk(n) * fdl(n) * abs (cz - bz))
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.
140 continue
c---- Tested tol.
endif
c.... See if the local triangular coordinates are to be returned.
c---- Find local triangular coordinates.
if (kloc .eq. 1) then
c---- Loop over subset of data.
do 150 n = 1, ns
nn = n + n1 - 1
wa(nn) = 1.0 - fdk(n)
wb(nn) = fdk(n) * (1.0 - fdl(n))
wc(nn) = fdk(n) * fdl(n)
c---- End of loop over subset of data.
150 continue
c---- Tested kloc.
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 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 160 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 160 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 9902 format (/ 'apttlod results:')
cbug 9903 format (i3,' px,py,pz=',1p3e22.14)
cbug 9904 format (/ 'apttlod mean and deviation:' /
cbug & ' avg(x,y,z)=',1p3e22.14 /
cbug & ' dev(x,y,z)=',1p3e22.14)
cbug
cbug write ( 3, 9902)
cbug write ( 3, 9903) (n, px(n), py(n), pz(n), n = 1, np)
cbug write ( 3, 9904) avgx, avgy, avgz, devx, devy, devz
cbug
cbugc---- Local coordinates returned.
cbug if (kloc .eq. 1) then
cbug
cbug avgwa = 0.0
cbug avgwb = 0.0
cbug avgwc = 0.0
cbug avgwa2 = 0.0
cbug avgwb2 = 0.0
cbug avgwc2 = 0.0
cbug
cbug do 170 n = 1, np
cbug avgwa = avgwa + wa(n)
cbug avgwb = avgwb + wb(n)
cbug avgwc = avgwc + wc(n)
cbug avgwa2 = avgwa2 + wa(n)**2
cbug avgwb2 = avgwb2 + wb(n)**2
cbug avgwc2 = avgwc2 + wc(n)**2
cbug 170 continue
cbug
cbug avgwa = avgwa / np
cbug avgwb = avgwb / np
cbug avgwc = avgwc / np
cbug avgwa2 = avgwa2 / np
cbug avgwb2 = avgwb2 / np
cbug avgwc2 = avgwc2 / np
cbug devwa = sqrt (avgwa2 - avgwa**2)
cbug devwb = sqrt (avgwb2 - avgwb**2)
cbug devwc = sqrt (avgwc2 - avgwc**2)
cbug
cbug 9905 format (i3,' wa,wb,wc=',1p3e22.14)
cbug 9906 format (/ 'apttlod mean and deviation:' /
cbug & ' avgw(a,b,c)=',1p3e22.14 /
cbug & ' devw(a,b,c)=',1p3e22.14)
cbug
cbug write ( 3, 9905) (n, wa(n), wb(n), wc(n), n = 1, np)
cbug write ( 3, 9906) avgwa, avgwb, avgwc,
cbug & devwa, devwb, devwc
cbug
cbugc---- Tested kloc.
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine apttlod. (+1 line.)
end
UCRL-WEB-209832