subroutine aptript (px, py, pz, ax, ay, az, bx, by, bz, & cx, cy, cz, np, tol, dabcp, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTRIPT c c call aptript (px, py, pz, ax, ay, az, bx, by, bz, c & cx, cy, cz, np, tol, dabcp, nerr) c c Version: aptript Updated 1990 January 18 14:20. c aptript 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 dabcp 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). 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: dabcp, nerr. c c Calls: aptvdis, aptvdot, aptvpln, aptvuna c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a". Size np. c c bx,by,bz Input The x, y, z coordinates of point "b". Size np. c c cx,cy,cz Input The x, y, z coordinates of point "c". Size np. c c dabcp Output Distance from the plane "abc" to point "p". Size np. c Positive if point "p" is on the side of the plane c for which points "a", "b" and "c" are in c counterclockwise order. c Zero if point "p" is in the plane of triangle "abc", c or if points "a", "b" and "c" are congruent or lie c in a straight line. c A value less than the estimated error in its c calculation is truncated to zero. 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 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 "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 plane. dimension dabcp (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.... Local variables. c---- Component x of vector "ap". common /laptript/ apx (64) c---- Component y of vector "ap". common /laptript/ apy (64) c---- Component z of vector "ap". common /laptript/ apz (64) c---- Length of vector "ap". common /laptript/ dap (64) c---- Index in local array. common /laptript/ n c---- First index of subset of data. common /laptript/ n1 c---- Last index of subset of data. common /laptript/ n2 c---- Size of current subset of data. common /laptript/ ns c---- Magnitude of normal vector. common /laptript/ vlen (64) c---- Component x of normal vector. common /laptript/ vnx (64) c---- Component y of normal vector. common /laptript/ vny (64) c---- Component z of normal vector. common /laptript/ vnz (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptript finding distance from abc to p.' / cbug & (i3,' ax,ay,az=',1p3e22.14 / cbug & ' bx,by,bz=',1p3e22.14 / cbug & ' cx,cy,cz=',1p3e22.14 / cbug & ' px,py,pz=',1p3e22.14)) cbug write ( 3, 9901) (n, ax(n), ay(n), az(n), bx(n), by(n), bz(n), cbug & cx(n), cy(n), cz(n), px(n), py(n), pz(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. 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 normal vector, which is upward from the side of the c.... plane for which points "a", "b" and "c" are counterclockwise. call aptvpln (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), & cx(n1), cy(n1), cz(n1), ns, tol, & vnx, vny, vnz, vlen, nerr) call aptvuna (vnx, vny, vnz, ns, 0.0, vlen, 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 (vnx, vny, vnz, apx, apy, apz, ns, tol, & dabcp(n1), nerr) 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 (/ 'aptript results:' / cbug & (i3,' dabcp=',1pe22.14)) cbug write ( 3, 9902) (n, dabcp(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptript. (+1 line.) end UCRL-WEB-209832