subroutine aptlntr (px, py, pz, sx, sy, sz, ax, ay, az, & bx, by, bz, cx, cy, cz, np, tol, & dpmin, dint, fracps, qx, qy, qz, ipar, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTLNTR c c call aptlntr (px, py, pz, sx, sy, sz, ax, ay, az, c & bx, by, bz, cx, cy, cz, np, tol, c & dpmin, dint, fracps, qx, qy, qz, ipar, nerr) c c Version: aptlntr Updated 1990 May 25 16:20. c aptlntr Originated 1990 May 25 16:20. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for each of np sets of input data, the point of c intersection of the line through points p = (px, py, pz) c and s = (sx, sy, sz), with the plane through the points c a = (ax, ay, az), b = (bx, by, bz) and c = (cx, cy, cz). c The point of intersection will be defined by its distance dint c from point "p", its fractional distance fracps along the line c from "p" to "s", and its coordinates q = (qx, qy, qz). c The perpendicular distance dpmin from the plane to point "p" c is also returned. c c If point "p" coincides with point "s", based on tol, the result c will be the same as if line "ps" is parallel to the plane. c If points "a", "b" and "c" are colinear or congruent, based on c tol, the result will be as if line "ps" lies in the plane. c If line "ps" is parallel to the plane, ipar will be 1. If, in c addition, dpmin is not zero, dint, fracps, and the coordinates c of "q" will be very large. If the line is parallel to the plane c and dpmin is zero, then the line is in the plane, and dint and c fracps will be zero, and the coordinates of "q" will be c (px, py, pz). c c Flag nerr indicates any input error. c c Input: px, py, pz, sx, sy, sz, ax, ay, az, bx, by, bz, cx, cy, cz, c np, tol. c c Output: dpmin, dint, fracps, qx, qy, qz, ipar, nerr. c c Calls: aptlnpl, aptvpln c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a" in the plane. c Size np. c c bx,by,bz Input The x, y, z coordinates of point "b" in the plane. c Size np. c c cx,cy,cz Input The x, y, z coordinates of point "c" in the plane. c Size np. c c dint Output The distance of the point of intersection "q" from c point "p". Positive if in the same direction as c that from "p" to "s". Size np. c Meaningless if ipar is not zero. c c dpmin Output The perpendicular distance to point "p" from the c plane. Positive if point "p" is in the same c direction from the plane as the normal vector "vn". c If less than the estimated error in its calculation, c dpmin will be truncated to zero. Size np. c Meaningless if ipar = 2, 3, or 4. c c fracps Output Fractional distance of point "q" along the line c segment from point "p" to point "s". Size np. c May be negative or greater than 1. c Meaningless if ipar is not zero. c c ipar Output 0 if the line is not parallel to the plane. Size np. c 1 if it is. See dpmin, dint, fracps, qx, qy, qz. c 2 if line "ps" is too short, based on tol. c 3 if triangle "abc" is too small, based on tol. c 4 if line "ps" is too short, and triangle "abc" is c too small, based on tol. 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. c c px,py,pz Input The x, y, z coordinates of point "p" on the line. c Must differ from "s", based on tol. Size np. c c qx,qy,qz Output The x, y, z coordinates of the point of intersection c of the line through "p" and "s" with the plane c through point "a" with normal vector "vn". c Meaningless if ipar is not zero. Size np. c c sx,sy,sz Input The x, y, z coordinates of point "s" on the line. c Must differ from "p", based on tol. 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 point "q". dimension dint (1) c---- Distance from plane to point "p". dimension dpmin (1) c---- Fractional distance of "q" along "ps". dimension fracps (1) c---- Line is parallel to plane, if 1. dimension ipar (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---- Coordinate x of point "s". dimension sx (1) c---- Coordinate y of point "s". dimension sy (1) c---- Coordinate z of point "s". dimension sz (1) c.... Local variables. c---- First index of subset of data. common /laptlntr/ n1 c---- Last index of subset of data. common /laptlntr/ n2 c---- Size of current subset of data. common /laptlntr/ ns c---- Magnitude of normal vector "vn".. common /laptlntr/ vlen (64) c---- Component x of vector "vn". common /laptlntr/ vnx (64) c---- Component y of vector "vn". common /laptlntr/ vny (64) c---- Component z of vector "vn". common /laptlntr/ vnz (64) cbugc***DEBUG begins. cbugc---- Index in arrays. cbug common /laptlntr/ n cbug 9901 format (/ 'aptlntr finding intersection of the line through:' / cbug & (i3,' px,py,pz=',1p3e22.14 / cbug & ' sx,sy,sz=',1p3e22.14 / cbug & ' with the plane through:' / 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), sx(n), sy(n), sz(n), cbug & ax(n), ay(n), az(n), bx(n), by(n), bz(n), cbug & cx(n), cy(n), cz(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 vector "vn" normal to the plane of triangle "abc". 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) c.... Find the distances, point of intersection, etc. call aptlnpl (px(n1), py(n1), pz(n1), sx(n1), sy(n1), sz(n1), & ax(n1), ay(n1), az(n1), vnx, vny, vnz, ns, tol, & dpmin(n1), dint(n1), fracps(n1), & qx(n1), qy(n1), qz(n1), ipar(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 (/ 'aptlntr results:' / cbug & (i3,' dpmin= ',1pe22.14,' ipar=',i2 / cbug & ' dint= ',1pe22.14,' fracps=',1pe22.14 / cbug & ' qx,qy,qz=',1p3e22.14)) cbug write ( 3, 9902) (n, dpmin(n), ipar(n), dint(n), fracps(n), cbug & qx(n), qy(n), qz(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptlntr. (+1 line.) end UCRL-WEB-209832