subroutine aptlnpl (px, py, pz, sx, sy, sz, ax, ay, az, & vnx, vny, vnz, np, tol, dpmin, dint, fracps, & qx, qy, qz, ipar, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTLNPL c c call aptlnpl (px, py, pz, sx, sy, sz, ax, ay, az, c & vnx, vny, vnz, np, tol, dpmin, dint, fracps, c & qx, qy, qz, ipar, nerr) c c Version: aptlnpl Updated 1990 November 29 11:40. c aptlnpl 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, 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 point c a = (ax, ay, az) with normal vector vn = (vnx, vny, vnz). 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 vector "vn" is too short, based on tol, the result c will be the same 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 History: 1990 March 14. Changed tol to 0.0 in call to unit vector c subroutine. Allows small magnitudes. c 1990 March 15. Changed results when vector "vn" is too short. c Now gives same results as if line "ps" is in the plane. c c Input: px, py, pz, sx, sy, sz, ax, ay, az, vnx, vny, vnz, np, tol. c c Output: dpmin, dint, fracps, qx, qy, qz, ipar, nerr. c c Calls: aptvadd, aptvdis, aptvdot, aptvuna, aptvunb 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 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 vector "vn" is too short, based on tol. c 4 if "ps" and "vn" are both too short, 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 c vnx,y,z Input The x, y, z components of vector "vn" normal to the c plane. Magnitude must exceed tol. Size np. 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---- 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---- Component x of vector "vn". dimension vnx (1) c---- Component y of vector "vn". dimension vny (1) c---- Component z of vector "vn". dimension vnz (1) c.... Local variables. c---- Component x of vector "ap". common /laptlnpl/ apx (64) c---- Component y of vector "ap". common /laptlnpl/ apy (64) c---- Component z of vector "ap". common /laptlnpl/ apz (64) c---- Cosine of angle between "ps" and "vn". common /laptlnpl/ costh (64) c---- A very small number. common /laptlnpl/ fuz c---- Index in arrays. common /laptlnpl/ n c---- First index of subset of data. common /laptlnpl/ n1 c---- Last index of subset of data. common /laptlnpl/ n2 c---- Index in external array. common /laptlnpl/ nn c---- Size of current subset of data. common /laptlnpl/ ns c---- Component x of unit normal vector. common /laptlnpl/ unx (64) c---- Component y of unit normal vector. common /laptlnpl/ uny (64) c---- Component z of unit normal vector. common /laptlnpl/ unz (64) c---- Component x of unit vector along "ps". common /laptlnpl/ upsx (64) c---- Component y of unit vector along "ps". common /laptlnpl/ upsy (64) c---- Component z of unit vector along "ps". common /laptlnpl/ upsz (64) c---- Magnitude of a vector. common /laptlnpl/ vlen (64) c---- Magnitude of vector "ap". common /laptlnpl/ vlenap (64) c---- Magnitude of vector "ps". common /laptlnpl/ vlenps (64) c---- Magnitude of normal vector "vn". common /laptlnpl/ vlenvn (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptlnpl 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 & ' with normal vector:' / cbug & ' vn(x,y,z)=',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), vnx(n), vny(n), vnz(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. 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 "ups" along line "ps". call aptvdis (px(n1), py(n1), pz(n1), sx(n1), sy(n1), sz(n1), ns, & tol, upsx, upsy, upsz, vlenps, nerr) call aptvuna (upsx, upsy, upsz, ns, 0., vlen, nerr) c.... Find the unit vector normal to the plane. call aptvunb (vnx(n1), vny(n1), vnz(n1), ns, 0., & unx, uny, unz, vlenvn, nerr) c.... Find the vector "ap" 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, vlenap, nerr) c.... Find the perpendicular distance from point "p" to the plane. call aptvdot (apx, apy, apz, unx, uny, unz, ns, tol, & dpmin(n1), nerr) c.... Find the component of unit vector "ups" perpendicular to the plane. call aptvdot (upsx, upsy, upsz, unx, uny, unz, ns, tol, & costh, nerr) c.... Test for the line being parallel to the plane, and for "ps" too short, c.... and for "vn" too short. In the later case, set dpmin = 0.0. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 if (abs (costh(n)) .gt. tol) then ipar(nn) = 0 else ipar(nn) = 1 endif if (vlenps(n) .le. tol) then ipar(nn) = 2 endif if (vlenvn(n) .le. tol) then ipar(nn) = 3 endif if ((vlenps(n) .le. tol) .and. & (vlenvn(n) .le. tol)) then ipar(nn) = 4 endif if (vlenvn(n) .le. tol) then dpmin(nn) = 0.0 endif c---- End of loop over subset of data. 120 continue c.... Find the distance and fractional distance of the point of intersection c.... along line "ps", from point "p". c---- Loop over subset of data. do 130 n = 1, ns nn = n + n1 - 1 dint(nn) = -dpmin(nn) / (costh(n) + fuz) fracps(nn) = dint(nn) / (vlenps(n) + fuz) c---- End of loop over subset of data. 130 continue c.... Find the coordinates of the point of intersection. call aptvadd (px(n1), py(n1), pz(n1), 1., dint, & upsx, upsy, upsz, ns, tol, & qx(n1), qy(n1), qz(n1), vlen, 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 (/ 'aptlnpl 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 aptlnpl. (+1 line.) end UCRL-WEB-209832