subroutine apttrip (px, py, pz, ax, ay, az, bx, by, bz,
& cx, cy, cz, noptfd, tol,
& dpmin, wa, wb, wc, qx, qy, qz,
& nlima, nlimb, nlimc, itrun, nside, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTTRIP
c
c call apttrip (px, py, pz, ax, ay, az, bx, by, bz,
c & cx, cy, cz, noptfd, tol,
c & dpmin, wa, wb, wc, qx, qy, qz,
c & nlima, nlimb, nlimc, itrun, nside, nerr)
c
c Version: apttrip Updated 1990 October 23 15:30.
c apttrip Originated 1989 November 2 14:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the distance dpmin from a point p = (px, py, pz)
c to a plane defined by the three points a = (ax, ay, az),
c b = (bx, by, bz), and c = (cx, cy, cz), and the point
c q = (qx, qy, qz) nearest to point p, and in the plane,
c subject to constraints that may be imposed by option noptfd
c and the value of tol.
c
c Optionally, to find the fractional distances wa, wb and wc
c of point "q" along the triangle's altitudes:
c
c q = wa * a + wb * b + wc * c.
c
c Flags nlima, nlimb, nlimc indicate when wa, wb and wc have been
c restrained. Flag itrun indicates when dpmin has been truncated
c to zero. Flag nside indicates when the minimum point is inside
c the triangle. Flag nerr indicates any input error.
c
c Input: px, py, pz, ax, ay, az, bx, by, bz, cx, cy, cz, noptfd, tol.
c
c Output: dpmin, wa, wb, wc, qx, qy, qz,
c nlima, nlimb, nlimc, itrun, nside, nerr.
c
c Calls: aptfdad, aptptln, aptvdis, aptvpln
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a".
c
c bx,by,bz Input The x, y, z coordinates of point "b".
c
c cx,cy,cz Input The x, y, z coordinates of point "c".
c
c dpmin Output Distance from point "p" to the nearest (constrained)
c point in the plane defined by points "a", "b", "c".
c a value less than the estimated error in its
c calculation is truncated to zero (itrun = 1).
c dpmin is positive when the external point is on the
c side of the plane for which the three points are in
c counterclockwise order. See noptfd.
c
c itrun Output 0 if no change is made in the calculated value of
c dpmin, 1 if dpmin is changed to zero, when less than
c the estimated error in its calculation.
c
c nerr Output Indicates an input error, if not 0.
c 1 is added if noptfd is not between 0 and 2.
c 2 is added if 3 points representing triangle are
c colinear or congruent.
c
c nlima Output 0 if no limit imposed on wa, 1 if the limit of
c noptfd = 1 is imposed, 2 if the limit of noptfd = 2
c is imposed.
c
c nlimb,c Output Like nlima, but for wb, wc, respectively.
c
c noptfd Input Option to limit range of wa, wb, wc:
c -1 for no limit, no calculation of wa, wb, wc,
c 0 for no limit,
c 1 to increase to tol, if in the range from -tol
c to tol, and decrease to 1.0 - tol, if in the
c range from 1.0 - tol to 1.0 + tol (move a point
c near an edge slightly inside the triangle), and
c 2 to limit to the range from 0.0 to 1.0 (move a
c point outside the triangle to an edge).
c
c nside Output 0 if minimum point outside the triangle, 1 if inside.
c 0 if moved to edge, when noptfd = 2.
c
c px,py,pz Input The x, y, z coordinates of point "p".
c
c qx,qy,qz Output The x, y, z coordinate of the point in the plane
c nearest point "p". May be constrained by option
c noptfd.
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"
c from edge "bc" to vertex "a".
c
c wb Output Fractional distance of point "q"
c from edge "ca" to vertex "b".
c
c wc Output Fractional distance of point "q"
c from edge "ab" to vertex "c".
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c.... (None.)
c.... Local variables.
c---- A very big number.
common /lapttrip/ big
c---- Is dpmin after moving wa, wb, wc.
common /lapttrip/ dpminx
c---- Is dpmin after moving wa.
common /lapttrip/ dpmina
c---- Is dpmin after moving wb.
common /lapttrip/ dpminb
c---- Is dpmin after moving wc.
common /lapttrip/ dpminc
c---- Sign of dpmin.
common /lapttrip/ dpsign
c---- Is qx - px.
common /lapttrip/ dx
c---- Is qy - py.
common /lapttrip/ dy
c---- Is qz - pz.
common /lapttrip/ dz
c---- Fractional distance of p along edge.
common /lapttrip/ fds
c---- Component of vector "ap" along vn.
common /lapttrip/ fract
c---- 2 if point moved to end of edge.
common /lapttrip/ nlim
c---- Is qx after moving wa.
common /lapttrip/ qxa
c---- Is qx after moving wb.
common /lapttrip/ qxb
c---- Is qx after moving wc.
common /lapttrip/ qxc
c---- Is qy after moving wa.
common /lapttrip/ qya
c---- Is qy after moving wb.
common /lapttrip/ qyb
c---- Is qy after moving wc.
common /lapttrip/ qyc
c---- Is qz after moving wa.
common /lapttrip/ qza
c---- Is qz after moving wb.
common /lapttrip/ qzb
c---- Is qz after moving wc.
common /lapttrip/ qzc
c---- Sum wa + wb + wc.
common /lapttrip/ sum
c---- Magnitude of normal vector.
common /lapttrip/ vlen
c---- Square of length of normal vector.
common /lapttrip/ vnsq
c---- Component x of normal vector.
common /lapttrip/ vnx
c---- Component y of normal vector.
common /lapttrip/ vny
c---- Component z of normal vector.
common /lapttrip/ vnz
cbugc***DEBUG begins.
cbug 9900 format (/ 'apttrip finding distance from triangle',
cbug & ' through points:' /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' to point:' /
cbug & ' px,py,pz=',1p3e22.14)
cbug write ( 3, 9900) ax, ay, az, bx, by, bz, cx, cy, cz,
cbug & px, py, pz
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
itrun = 0
nerr = 0
nlima = 0
nlimb = 0
nlimc = 0
nside = 0
c.... Test for input errors.
c++++ Bad input.
if ((noptfd .lt. -1) .or. (noptfd .gt. 2)) then
nerr = nerr + 1
endif
c.... Find the components of the normal vector, which is upward from the
c.... side of the plane for which points (a, b, c) are counterclockwise.
call aptvpln (ax, ay, az, bx, by, bz, cx, cy, cz, 1, tol,
& vnx, vny, vnz, vlen, nerr)
vnsq = vlen**2
cbugc***DEBUG begins.
cbug 9901 format (/ 'apttrip' /
cbug & 'vnx=',1pe13.5,' vny=',1pe13.5,' vnz=',1pe13.5,' vlen=',1pe13.5)
cbug write ( 3, 9901) vnx, vny, vnz, vlen
cbugc***DEBUG ends.
c---- Null triangle.
if (vlen .le. tol) then
nerr = nerr + 2
endif
c---- An input error was found.
if (nerr .ne. 0) go to 210
c.... Find the fractional distance along the normal vector from the
c.... plane to the external point.
fract = ((px - ax) * vnx + (py - ay) * vny +
& (pz - az) * vnz) / vnsq
dpsign = sign (1.0, fract)
c.... Find the distance from the plane to (px, py, pz).
dpmin = fract * vlen
if (abs (fract) .lt. tol) then
itrun = 1
dpmin = 0.0
fract = 0.0
endif
c.... Find the point in the plane nearest to (px, py, pz).
qx = px - fract * vnx
qy = py - fract * vny
qz = pz - fract * vnz
cbugc***DEBUG begins.
cbug 9902 format ('fract=',1pe13.5 /
cbug & ' initial qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & ' initial dpmin=',1pe22.14,' itrun=',i2)
cbug write ( 3, 9902) fract, qx, qy, qz, dpmin, itrun
cbugc***DEBUG ends.
c.... See which option is specified for adjusting point "q".
if (noptfd .eq. -1) go to 210
c.... Find the fractional distances of the point along the altitudes.
wa = (((cy - by)*(qz - cz) - (cz - bz)*(qy - cy)) * vnx +
& ((cz - bz)*(qx - cx) - (cx - bx)*(qz - cz)) * vny +
& ((cx - bx)*(qy - cy) - (cy - by)*(qx - cx)) * vnz) /
& vnsq
wb = (((ay - cy)*(qz - az) - (az - cz)*(qy - ay)) * vnx +
& ((az - cz)*(qx - ax) - (ax - cx)*(qz - az)) * vny +
& ((ax - cx)*(qy - ay) - (ay - cy)*(qx - ax)) * vnz) /
& vnsq
wc = (((by - ay)*(qz - bz) - (bz - az)*(qy - by)) * vnx +
& ((bz - az)*(qx - bx) - (bx - ax)*(qz - bz)) * vny +
& ((bx - ax)*(qy - by) - (by - ay)*(qx - bx)) * vnz) /
& vnsq
cbugc***DEBUG begins.
cbug 9903 format (/ 'apttrip' /
cbug & 'wa=',1pe13.5,' wb=',1pe13.5,' wc=',1pe13.5,
cbug & ' before adjusting' / 'sum=',1pe22.14)
cbug sum = wa + wb + wc
cbug write ( 3, 9903) wa, wb, wc, sum
cbugc***DEBUG ends.
c.... Adjust the fractional distances according to noptfd.
c---- No adjustment of wa, wb, wc.
if (noptfd .eq. 0) then
nside = 1
if ((wa .lt. 0.0) .or. (wa .gt. 1.0)) then
nside = 0
endif
if ((wb .lt. 0.0) .or. (wb .gt. 1.0)) then
nside = 0
endif
if ((wc .lt. 0.0) .or. (wc .gt. 1.0)) then
nside = 0
endif
c---- Adjust points near triangle edges.
elseif (noptfd .eq. 1) then
nside = 1
call aptfdad (wa, 1, tol, nlima, nerr)
if ((wa .lt. tol) .or. (wa .gt. (1.0 - tol))) then
nside = 0
endif
call aptfdad (wb, 1, tol, nlimb, nerr)
if ((wb .lt. tol) .or. (wb .gt. (1.0 - tol))) then
nside = 0
endif
call aptfdad (wc, 1, tol, nlimc, nerr)
if ((wc .lt. tol) .or. (wc .gt. (1.0 - tol))) then
nside = 0
endif
c++++ Adjusted wa, wb, wc.
if ((nlima + nlimb + nlimc) .gt. 0) then
cbugc***DEBUG begins.
cbug 9904 format (/ 'apttrip' /
cbug & 'wa=',1pe13.5,' wb=',1pe13.5,' wc=',1pe13.5,' after adjusting')
cbug write ( 3, 9904) wa, wb, wc
cbugc***DEBUG ends.
c.... Renormalize wa, wb, wc.
sum = wa + wb + wc
wa = wa / sum
wb = wb / sum
wc = wc / sum
c.... Recalculate point "q".
qx = wa * ax + wb * bx + wc * cx
qy = wa * ay + wb * by + wc * cy
qz = wa * az + wb * bz + wc * cz
c.... Recalculate dpmin, with correct sign.
call aptvdis (px, py, pz, qx, qy, qz, 1, tol,
& dx, dy, dz, dpmin, itrun)
dpmin = dpsign * dpmin
cbugc***DEBUG begins.
cbug 9905 format (/ 'apttrip' /
cbug & 'wa=',1pe13.5,' wb=',1pe13.5,' wc=',1pe13.5,
cbug & ' after sum = 1.0' /
cbug & 10x,'qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & 9x,'dpmin=',1pe13.5)
cbug write ( 3, 9905) wa, wb, wc, qx, qy, qz, dpmin
cbugc***DEBUG ends.
c---- Adjusted wa, wb, wc.
endif
c---- Move outside point to nearest edge.
elseif (noptfd .eq. 2) then
dpminx = big
c---- Impose limit 0.0 .le. wa .le. 1.0.
if (wa .lt. 0.0) then
nlima = 2
call aptptln (px, py, pz, bx, by, bz, cx, cy, cz, 1, tol,
& noptfd, dpmina, fds, qxa, qya, qza,
& nlim, itrun, nerr)
if (dpmina .lt. dpminx) then
dpminx = dpmina
qx = qxa
qy = qya
qz = qza
endif
cbugc***DEBUG begins.
cbug 9906 format (/ 'apttrip' /
cbug & ' new wa. qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & ' new dpmin=',1pe22.14)
cbug write ( 3, 9906) qxa, qya, qza, dpmina
cbugc***DEBUG ends.
c---- Imposed limit on wa.
endif
c---- Impose limit 0.0 .le. wb .le. 1.0.
if (wb .lt. 0.0) then
nlimb = 2
call aptptln (px, py, pz, ax, ay, az, cx, cy, cz, 1, tol,
& noptfd, dpminb, fds, qxb, qyb, qzb,
& nlim, itrun, nerr)
if (dpminb .lt. dpminx) then
dpminx = dpminb
qx = qxb
qy = qyb
qz = qzb
endif
cbugc***DEBUG begins.
cbug 9907 format (/ 'apttrip' /
cbug & ' new wb. qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & ' new dpmin=',1pe22.14)
cbug write ( 3, 9907) qxb, qyb, qzb, dpminb
cbugc***DEBUG ends.
c---- Imposed limit on wb.
endif
c---- Impose limit 0.0 .le. wc .le. 1.0.
if (wc .lt. 0.0) then
nlimc = 2
call aptptln (px, py, pz, ax, ay, az, bx, by, bz, 1, tol,
& noptfd, dpminc, fds, qxc, qyc, qzc,
& nlim, itrun, nerr)
if (dpminc .lt. dpminx) then
dpminx = dpminc
qx = qxc
qy = qyc
qz = qzc
endif
cbugc***DEBUG begins.
cbug 9908 format (/ 'apttrip' /
cbug & ' new wc. qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & ' new dpmin=',1pe22.14)
cbug write ( 3, 9908) qxc, qyc, qzc, dpminc
cbugc***DEBUG ends.
c---- Imposed limit on wc.
endif
c.... Recalculate wa, wb, wc if nearest point moved to a triangle edge.
c.... Adjust as necessary to reduce truncation error.
nside = 1
c++++ Recalculate wa, wb, wc.
if ((nlima + nlimb + nlimc) .gt. 0) then
nside = 0
wa = (((cy - by)*(qz - cz) - (cz - bz)*(qy - cy)) * vnx +
& ((cz - bz)*(qx - cx) - (cx - bx)*(qz - cz)) * vny +
& ((cx - bx)*(qy - cy) - (cy - by)*(qx - cx)) * vnz)/
& vnsq
if (wa .lt. tol) wa = 0.0
wb = (((ay - cy)*(qz - az) - (az - cz)*(qy - ay)) * vnx +
& ((az - cz)*(qx - ax) - (ax - cx)*(qz - az)) * vny +
& ((ax - cx)*(qy - ay) - (ay - cy)*(qx - ax)) * vnz)/
& vnsq
if (wb .lt. tol) wb = 0.0
wc = (((by - ay)*(qz - bz) - (bz - az)*(qy - by)) * vnx +
& ((bz - az)*(qx - bx) - (bx - ax)*(qz - bz)) * vny +
& ((bx - ax)*(qy - by) - (by - ay)*(qx - bx)) * vnz)/
& vnsq
if (wc .lt. tol) wc = 0.0
c.... Renormalize wa, wb, wc.
sum = wa + wb + wc
wa = wa / sum
wb = wb / sum
wc = wc / sum
c.... Recalculate point "q".
qx = wa * ax + wb * bx + wc * cx
qy = wa * ay + wb * by + wc * cy
qz = wa * az + wb * bz + wc * cz
c.... Recalculate dpmin, with correct sign.
call aptvdis (px, py, pz, qx, qy, qz, 1, tol,
& dx, dy, dz, dpmin, itrun)
dpmin = dpsign * dpmin
cbugc***DEBUG begins.
cbug 9909 format (/ 'apttrip' /
cbug & 'wa=',1pe13.5,' wb=',1pe13.5,' wc=',1pe13.5,' after moving')
cbug 9910 format (6x,'new qx=',1pe13.5,' qy=',1pe13.5,' qz=',1pe13.5 /
cbug & ' new dpmin=',1pe22.14)
cbug write ( 3, 9909) wa, wb, wc
cbug write ( 3, 9910) qx, qy, qz, dpmin
cbugc***DEBUG ends.
c---- Recalculated wa, wb, wc.
endif
c---- Noptfd options.
endif
210 return
c.... End of subroutine apttrip. (+1 line.)
end
UCRL-WEB-209832