subroutine aptrkpl (px, py, pz, vx, vy, vz, ax, ay, az,
& bx, by, bz, cx, cy, cz, np, tol,
& dpmin, dint, qx, qy, qz, ipar, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRKPL
c
c call aptrkpl (px, py, pz, vx, vy, vz, ax, ay, az,
c & bx, by, bz, cx, cy, cz, np, tol,
c & dpmin, dint, qx, qy, qz, ipar, nerr)
c
c Version: aptrkpl Updated 1990 November 29 17:40.
c aptrkpl Originated 1989 November 30 15:50.
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 q = (qx, qy, qz) of the linear track through point
c p = (px, py, pz) with direction vector v = (vx, vy, vz),
c and the plane through points a = (ax, ay, az), b = (bx, by, bz),
c and c = (cx, cy, cz), and the distance dint between points
c "p" and "q". The perpendicular distance dpmin from the plane to
c point "p" is also returned.
c
c If the vector "v" is parallel to the plane, flag ipar will
c be 1. If so, and dpmin is not zero, dint and the coordinates
c of point "q" will be very large. Otherwise, if dpmin is zero,
c dint will be zero, and the coordinates of point "q" will be
c those of point "p".
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 the values of dpmin, dint, and point "q"
c when the plane "abc" is undefined (ipar = 3 or 4). No effect
c on problems with good input data.
c
c Input: px, py, pz, vx, vy, vz, ax, ay, az, bx, by, bz, cx, cy, cz,
c np, tol.
c
c Output: dpmin, dint, qx, qy, qz, ipar, nerr.
c
c Calls: aptvadd, aptvdis, aptvdot, aptvpln, aptvuna, aptvunb
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a" in the plane.
c Must differ from "b" and "c", based on tol. Size np.
c
c bx,by,bz Input The x, y, z coordinates of point "b" in the plane.
c Must differ from "a" and "c", based on tol. Size np.
c
c cx,cy,cz Input The x, y, z coordinates of point "c" in the plane.
c Must differ from "a" and "b", based on tol. Size np.
c
c dint Output The distance of the point of intersection "q" from
c point "p", if ipar = 0. Positive if in the same
c direction as vector "v". Size np.
c Meaningless, but large if ipar = 1.
c Meaningless, but zero if ipar = 2, 3, or 4.
c
c dpmin Output The perpendicular distance to point "p" from the
c plane. Positive if point "p" is on the side of the
c plane for which points "a", "b", and "c" are
c counterclockwise. Meaningless if ipar = 2, 3, or 4.
c If less than the estimated error in its calculation,
c dpmin will be truncated to zero. Size np.
c
c ipar Output 0 if vector "v" is not parallel to the plane. Size np.
c 1 if it is. See dpmin, dint, qx, qy, qz.
c 2 if vector "v" is too short, based on tol.
c 3 if any of the points "a", "b", and "c" coincide.
c 4 if "v" is too short, and "a", "b" or "c" coincide.
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". Size np.
c
c qx,qy,qz Output The x, y, z coordinates of the point of intersection
c of the line through point "p" with unit direction
c vector "v", and the plane through points "a", "b",
c and "c", if ipar = 0. Size np.
c Meaningless, but large if ipar = 1.
c Meaningless, but zero if ipar = 2, 3 or 4.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c vx,vy,vz Input The x, y, z components of direction vector "v".
c 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---- 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---- Line 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---- Component x of vector "v".
dimension vx (1)
c---- Component y of vector "v".
dimension vy (1)
c---- Component z of vector "v".
dimension vz (1)
c.... Local variables.
c---- Component x of vector "ap".
common /laptrkpl/ apx (64)
c---- Component y of vector "ap".
common /laptrkpl/ apy (64)
c---- Component z of vector "ap".
common /laptrkpl/ apz (64)
c---- Cosine of angle between "v" and "un".
common /laptrkpl/ costh (64)
c---- A very small number.
common /laptrkpl/ fuz
c---- Index in arrays.
common /laptrkpl/ n
c---- First index of subset of data.
common /laptrkpl/ n1
c---- Last index of subset of data.
common /laptrkpl/ n2
c---- Index in external array.
common /laptrkpl/ nn
c---- Size of current subset of data.
common /laptrkpl/ ns
c---- Component x of unit normal vector.
common /laptrkpl/ unx (64)
c---- Component y of unit normal vector.
common /laptrkpl/ uny (64)
c---- Component z of unit normal vector.
common /laptrkpl/ unz (64)
c---- Component x of unit vector "u".
common /laptrkpl/ ux (64)
c---- Component y of unit vector "u".
common /laptrkpl/ uy (64)
c---- Component z of unit vector "u".
common /laptrkpl/ uz (64)
c---- Magnitude of a vector.
common /laptrkpl/ vlen (64)
c---- Magnitude of vector "ap".
common /laptrkpl/ vlenap (64)
c---- Magnitude of unit normal vector "un".
common /laptrkpl/ vlenun (64)
c---- Magnitude of vector "v".
common /laptrkpl/ vlenv (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrkpl finding intersection of the line through:' /
cbug & (i3,' px,py,pz=',1p3e22.14 /
cbug & ' with direction vector:' /
cbug & ' vx,vy,vz=',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), vx(n), vy(n), vz(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.
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 parallel to direction vector "v".
call aptvunb (vx(n1), vy(n1), vz(n1), ns, 0.,
& ux, uy, uz, vlenv, nerr)
c.... Find the unit vector normal to the plane "abc".
call aptvpln (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& cx(n1), cy(n1), cz(n1), ns, tol,
& unx, uny, unz, vlenun, nerr)
call aptvuna (unx, uny, unz, ns, 0., vlenun, 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.
c.... (Positive if point "p" is on the side of the plane for which
c.... points "a", "b" and "c" are counterclockwise.)
call aptvdot (apx, apy, apz, unx, uny, unz, ns, tol,
& dpmin(n1), nerr)
c.... Find the component of unit vector "u" perpendicular to the plane.
call aptvdot (ux, uy, uz, unx, uny, unz, ns, tol,
& costh, nerr)
c.... Find the distance of the point of intersection "q" from point "p",
c.... in the direction of unit vector "v".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
if (vlenun(n) .le. tol) then
dpmin(nn) = 0.0
endif
dint(nn) = -dpmin(nn) / (costh(n) + fuz)
c---- End of loop over subset of data.
120 continue
c.... Find the coordinates of the point of intersection.
call aptvadd (px(n1), py(n1), pz(n1), 1., dint,
& ux, uy, uz, ns, tol,
& qx(n1), qy(n1), qz(n1), vlen, nerr)
c.... Test for the direction vector "v" being parallel to the plane,
c.... and for the direction vector "v" being too short,
c.... and for 2 or 3 of the points "a", "b" and "c" being coincident,
c.... or all 3 being colinear.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
if (abs (costh(n)) .gt. tol) then
ipar(nn) = 0
else
ipar(nn) = 1
endif
if (vlenv(n) .le. tol) then
ipar(nn) = 2
endif
if (vlenun(n) .le. tol) then
ipar(nn) = 3
endif
if ((vlenv(n) .le. tol) .and.
& (vlenun(n) .le. tol) ) then
ipar(nn) = 4
endif
c---- End of loop over subset of data.
130 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 (/ 'aptrkpl results:' /
cbug & (i3,' dpmin= ',1pe22.14,' ipar=',i2 /
cbug & ' dint= ',1pe22.14 /
cbug & ' qx,qy,qz=',1p3e22.14))
cbug write ( 3, 9902) (n, dpmin(n), ipar(n), dint(n),
cbug & qx(n), qy(n), qz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptrkpl. (+1 line.)
end
UCRL-WEB-209832