subroutine aptripl (px, py, pz, ax, ay, az, bx, by, bz,
& cx, cy, cz, np, tol,
& dpmin, qx, qy, qz, wa, wb, wc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRIPL
c
c call aptripl (px, py, pz, ax, ay, az, bx, by, bz, cx, cy, cz,
c & np, tol, dpmin, qx, qy, qz, wa, wb, wc, nerr)
c
c Version: aptripl Updated 1990 November 29 15:30.
c aptripl Originated 1990 May 15 16:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np points p = (px, py, pz), the distance
c dpmin to the plane defined by the three points a = (ax, ay, az),
c b = (bx, by, bz) and c = (cx, cy, cz); the point
c q = (qx, qy, qz) in the plane nearest point "p"; and the local
c coordinates wa, wb, and wc of point "q" in triangle "abc".
c
c q = wa * a + wb * b + wc * c
c
c Flag nerr indicates any input error.
c
c Input: px, py, pz, ax, ay, az, bx, by, bz, cx, cy, cz, np, tol.
c
c Output: dpmin, qx, qy, qz, wa, wb, wc, nerr.
c
c Calls: aptvadd, aptvdis, aptvdot, aptvpln, aptvuna
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a" of the triangle.
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of the triangle.
c Size np.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of the triangle.
c Size np.
c
c dpmin Output Distance from point "p" to the nearest point "q" in the
c plane of triangle "abc". Truncated to zero if less
c than the estimated error in its calculation, based on
c tol. Positive when point "p" is on the side of the
c plane for which points "a", "b" and "c" are in
c counterclockwise order. Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c px,py,pz Input The x, y, z coordinates of point "p". Size np.
c
c qx,qy,qz Output The x, y, z coordinates of point "q", the point in the
c plane "abc" nearest point "p". Size np.
c
c tol Input Numerical tolerance limit for dpmin, wa, wb, wc.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c wa Output Fractional distance of point "q" from edge "bc" to
c vertex "a". Size np. Meaningless if triangle "abc"
c has zero area.
c
c wb Output Fractional distance of point "q" from edge "ca" to
c vertex "b". Size np. Meaningless if triangle "abc"
c has zero area.
c
c wc Output Fractional distance of point "q" from edge "ab" to
c vertex "c". Size np. Meaningless if triangle "abc"
c has zero area.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of point "a".
dimension ax (1)
c---- Coordinate y of point "a".
dimension ay (1)
c---- Coordinate z of point "a".
dimension az (1)
c---- Coordinate x of point "b".
dimension bx (1)
c---- Coordinate y of point "b".
dimension by (1)
c---- Coordinate z of point "b".
dimension bz (1)
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---- Distance from point "p" to point "q".
dimension dpmin (1)
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---- Coordinate x of point "q".
dimension qx (1)
c---- Coordinate y of point "q".
dimension qy (1)
c---- Coordinate z of point "q".
dimension qz (1)
c---- Fract. dist. of "q" from "bc" to "a".
dimension wa (1)
c---- Fract. dist. of "q" from "ca" to "b".
dimension wb (1)
c---- Fract. dist. of "q" from "ab" to "c".
dimension wc (1)
c.... Local variables.
c---- Component x of unit normal to "abc".
common /laptripl/ abcx (64)
c---- Component y of unit normal to "abc".
common /laptripl/ abcy (64)
c---- Component z of unit normal to "abc".
common /laptripl/ abcz (64)
c---- Component x of vector from "a" to "p".
common /laptripl/ apx (64)
c---- Component y of vector from "a" to "p".
common /laptripl/ apy (64)
c---- Component z of vector from "a" to "p".
common /laptripl/ apz (64)
c---- A very big number.
common /laptripl/ big
c---- Distance from point "a" to point "p".
common /laptripl/ dap (64)
c---- A very small number.
common /laptripl/ fuz
c---- Index in unit vector array.
common /laptripl/ n
c---- First index of subset of data.
common /laptripl/ n1
c---- Last index of subset of data.
common /laptripl/ n2
c---- Index in external array.
common /laptripl/ nn
c---- Size of current subset of data.
common /laptripl/ ns
c---- Component x of normal to "qab".
common /laptripl/ qabx (64)
c---- Component y of normal to "qab".
common /laptripl/ qaby (64)
c---- Component z of normal to "qab".
common /laptripl/ qabz (64)
c---- Component x of normal to "qbc".
common /laptripl/ qbcx (64)
c---- Component y of normal to "qbc".
common /laptripl/ qbcy (64)
c---- Component z of normal to "qbc".
common /laptripl/ qbcz (64)
c---- Component x of normal to "qca".
common /laptripl/ qcax (64)
c---- Component y of normal to "qca".
common /laptripl/ qcay (64)
c---- Component z of normal to "qca".
common /laptripl/ qcaz (64)
c---- Twice area of triangle "abc".
common /laptripl/ sabc (64)
c---- Value with sign of wc.
common /laptripl/ signab (64)
c---- Value with sign of wa.
common /laptripl/ signbc (64)
c---- Value with sign of wb.
common /laptripl/ signca (64)
c---- Twice area of triangle "qab".
common /laptripl/ sqab (64)
c---- Twice area of triangle "qbc".
common /laptripl/ sqbc (64)
c---- Twice area of triangle "qca".
common /laptripl/ sqca (64)
c---- Length of a vector.
common /laptripl/ vlen (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptripl finding local coordinates in triangle.' /
cbug & (i3,' px,py,pz=',1p3e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14))
cbug write ( 3, 9901) (n, px(n), py(n), pz(n), ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
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
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the unit vector normal to the plane of triangle "abc", and
c.... the area function of "abc".
call aptvpln (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& cx(n1), cy(n1), cz(n1), ns, tol,
& abcx, abcy, abcz, sabc, nerr)
call aptvuna (abcx, abcy, abcz, ns, 0.0, sabc, nerr)
c.... Find the distance vector from point "a" to point "p".
call aptvdis (ax(n1), ay(n1), az(n1), px(n1), py(n1), pz(n1),
& ns, tol, apx, apy, apz, dap, nerr)
c.... Find the distance from plane "abc" to point "p".
call aptvdot (abcx, abcy, abcz, apx, apy, apz, ns, tol,
& dpmin(n1), nerr)
c.... Find point "q" in the plane of triangle "abc" nearest point "p".
call aptvadd (px(n1), py(n1), pz(n1), -1., dpmin(n1),
& abcx, abcy, abcz, ns, tol,
& qx(n1), qy(n1), qz(n1), vlen, nerr)
c.... Find the normal vectors of the triangles formed by point "q" with the
c.... sides of the triangle, and the triangle area functions.
call aptvpln (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& qx(n1), qy(n1), qz(n1), ns, tol,
& qabx, qaby, qabz, sqab, nerr)
call aptvpln (qx(n1), qy(n1), qz(n1), bx(n1), by(n1), bz(n1),
& cx(n1), cy(n1), cz(n1), ns, tol,
& qbcx, qbcy, qbcz, sqbc, nerr)
call aptvpln (ax(n1), ay(n1), az(n1), qx(n1), qy(n1), qz(n1),
& cx(n1), cy(n1), cz(n1), ns, tol,
& qcax, qcay, qcaz, sqca, nerr)
c.... Find the signs of the "q" triangles relative to triangle "abc".
call aptvdot (abcx, abcy, abcz, qabx, qaby, qabz, ns, tol,
& signab, nerr)
call aptvdot (abcx, abcy, abcz, qbcx, qbcy, qbcz, ns, tol,
& signbc, nerr)
call aptvdot (abcx, abcy, abcz, qcax, qcay, qcaz, ns, tol,
& signca, nerr)
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
wa(nn) = sqbc(n) / (sabc(n) + fuz)
wb(nn) = sqca(n) / (sabc(n) + fuz)
wc(nn) = sqab(n) / (sabc(n) + fuz)
wa(nn) = sign (wa(nn), signbc(n))
wb(nn) = sign (wb(nn), signca(n))
wc(nn) = sign (wc(nn), signab(n))
if (sabc(n) .le. tol) then
wa(nn) = big
wb(nn) = big
wc(nn) = big
endif
c---- End of loop over subset of data.
120 continue
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 9902 format (/ 'aptripl results:' /
cbug & (i3,' dpmin= ',1pe22.14 /
cbug & ' qx,qy,qz=',1p3e22.14 /
cbug & ' wa,wb,wc=',1p3e22.14))
cbug write ( 3, 9902) (n, dpmin(n), qx(n), qy(n), qz(n),
cbug & wa(n), wb(n), wc(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptripl. (+1 line.)
end
UCRL-WEB-209832