subroutine aptplpl (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, np, tol, ex, ey, ez, fx, fy, fz,
& ux, uy, uz, ipar, dpmin, itrun, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPLPL
c
c call aptplpl (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, ex, ey, ez, fx, fy, fz,
c & ux, uy, uz, ipar, dpmin, itrun, nerr)
c
c Version: aptplpl Updated 1990 November 28 14:50.
c aptplpl Originated 1989 November 9 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of a set of input data, the line of
c intersection of the plane through point a = (ax, ay, az)
c with normal vector b = (bx, by, bz), and the plane through
c point c = (cx, cy, cz) with normal vector d = (dx, dy, dz),
c if any. Otherwise, if the planes are parallel, to find
c the distance between them. For nonparallel planes,
c the points e = (ex, ey, ez) and f = (fx, fy, fz) will be the
c points on the line of intersection nearest points "a" and "c",
c respectively, unit vector u = (ux, uy, uz) will be the
c direction of the line of intersection, and ipar will be zero.
c
c If the planes are parallel, based on tol, ipar will be 1,
c and dpmin will be the the distance between the planes. If the
c planes are coincident, dpmin will be zero, itrun will be 1,
c and unit vector "u" will be in the direction of the line "ac".
c
c Flag nerr indicates any input error, such as np not positive.
c Flag ipar will be 2, 3, or 4 if vector "b" or "d", or both,
c are too short, based on tol. If so, "e", "f" and "u" will be
c meaningless.
c
c Method: The line of intersection "ef" lies in both planes, therefore
c is perpendicular to both normal vectors, thus parallel to b x d.
c The line "ae" is perpendicular to line "ef" and "b", thus
c parallel to b x (b x d). The line "cf" is perpendicular to line
c "ef" and "d", thus parallel to d x (b x d). The vector path
c "ac" equals the vector path "aefc". Taking components parallel
c to "ae", "ef", and "cf" provides equations for the unknown
c distances "ae", "ef", and "cf". If the planes are parallel,
c b x d is zero, and dpmin is the component of "ac" in the
c direction of the normal vector "b".
c
c Note: Subroutine aptvpln may be used to find the vector normal to
c a plane for which at least 3 points are known.
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 points "e" and "f" to "a" and "b", resp.
c when either vector "b" or "d" is too small. No effect when
c input data is good.
c 1990 April 2. Fixed bug in indexing. Caused error in problems
c with np .gt. 64.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, tol.
c
c Output: ex, ey, ez, fx, fy, fz, ux, uy, uz, ipar, dpmin, itrun, nerr.
c
c Calls: aptvadd, aptvdis, aptvdot, aptvuna, aptvunb, aptvxun
c
c
c Glossary:
c
c ax,ay,az Input A point in plane "a". Size np.
c
c bx,by,bz Input A vector normal to plane "a". Size np.
c If too short, based on tol, ipar = 2 or 4.
c
c cx,cy,cz Input A point in plane "c". Size np.
c
c dpmin Output The distance between planes "a" and "c", if ipar = 1.
c Otherwise, zero. Size np.
c
c dx,dy,dz Input A vector normal to plane "c". Size np.
c If too short, based on tol, ipar = 3 or 4.
c
c ex,ey,ez Output The x, y, z coordinates of the point on the line of
c intersection of planes "a" and "c" nearest point "a".
c Meaningless, but point "a", if ipar is not 0.
c Size np.
c
c fx,fy,fz Output The x, y, z coordinates of the point on the line of
c intersection of planes "a" and "c" nearest point "c".
c Meaningless, but point "c", if ipar is not 0.
c Size np.
c
c ipar Output Indicates relative orientation of planes "a" and "c":
c 0 if the planes intersect.
c 1 if the planes are parallel, based on tol.
c 2 if vector "b" is too short, based on tol.
c 3 if vector "d" is too short, based on tol.
c 4 if vectors "b" and "d" are both too short.
c Orientation is indeterminate if ipar = 2, 3 or 4.
c Size np.
c
c itrun Output If 1, indicates planes are parallel and coincident.
c Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c ux,uy,uz Output The x, y, z components of the unit vector parallel to
c the line of intersection of planes "a" and "c",
c if ipar = 0. Meaningless, but parallel to the
c line "ac" if ipar is not zero. 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---- Component x of vector "b".
dimension bx (1)
c---- Component y of vector "b".
dimension by (1)
c---- Component z of vector "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 between parallel planes.
dimension dpmin (1)
c---- Component x of vector "d".
dimension dx (1)
c---- Component y of vector "d".
dimension dy (1)
c---- Component z of vector "d".
dimension dz (1)
c---- Coordinate x of point "e".
dimension ex (1)
c---- Coordinate y of point "e".
dimension ey (1)
c---- Coordinate z of point "e".
dimension ez (1)
c---- Coordinate x of point "f".
dimension fx (1)
c---- Coordinate y of point "f".
dimension fy (1)
c---- Coordinate z of point "f".
dimension fz (1)
c---- Indicates orientation of planes.
dimension ipar (1)
c---- 1 if planes coincident.
dimension itrun (1)
c---- Component x of unit vector "u".
dimension ux (1)
c---- Component y of unit vector "u".
dimension uy (1)
c---- Component z of unit vector "u".
dimension uz (1)
c.... Local variables.
c---- Cosine of angle between "b" and "d".
common /laptplpl/ costh (64)
c---- Distance from point "a" to point "c".
common /laptplpl/ dac (64)
c---- Dot product of vectors "vac", "uae".
common /laptplpl/ dacae (64)
c---- Dot product of vectors "vac", "ucf".
common /laptplpl/ daccf (64)
c---- Distance from point "a" to point "e".
common /laptplpl/ dae (64)
c---- Distance from point "c" to point "f".
common /laptplpl/ dcf (64)
c---- A very small number.
common /laptplpl/ fuz
c---- Index in arrays.
common /laptplpl/ n
c---- First index of subset of data.
common /laptplpl/ n1
c---- Last index of subset of data.
common /laptplpl/ n2
c---- Index in external array.
common /laptplpl/ nn
c---- Size of current subset of data.
common /laptplpl/ ns
c---- Sine of angle between "b" and "d".
common /laptplpl/ sinth (64)
c---- Component x of unit vector "uae".
common /laptplpl/ uaex (64)
c---- Component y of unit vector "uae".
common /laptplpl/ uaey (64)
c---- Component z of unit vector "uae".
common /laptplpl/ uaez (64)
c---- Component x of unit vector "ub".
common /laptplpl/ ubx (64)
c---- Component y of unit vector "ub".
common /laptplpl/ uby (64)
c---- Component z of unit vector "ub".
common /laptplpl/ ubz (64)
c---- Component x of unit vector "ucf".
common /laptplpl/ ucfx (64)
c---- Component y of unit vector "ucf".
common /laptplpl/ ucfy (64)
c---- Component z of unit vector "ucf".
common /laptplpl/ ucfz (64)
c---- Component x of unit vector "ud".
common /laptplpl/ udx (64)
c---- Component y of unit vector "ud".
common /laptplpl/ udy (64)
c---- Component z of unit vector "ud".
common /laptplpl/ udz (64)
c---- Component x of vector from "a" to "c".
common /laptplpl/ vacx (64)
c---- Component y of vector from "a" to "c".
common /laptplpl/ vacy (64)
c---- Component z of vector from "a" to "c".
common /laptplpl/ vacz (64)
c---- Magnitude of a vector.
common /laptplpl/ vlen (64)
c---- Magnitude of normal vector "b".
common /laptplpl/ vlenb (64)
c---- Magnitude of normal vector "d".
common /laptplpl/ vlend (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptplpl finding the line of intersection' /
cbug & ' of the plane through point a with normal vector b' /
cbug & ' and the plane through point c with normal vector d:')
cbug 9902 format (i3,' ax,ay,az= ',1p3e22.14 /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 /
cbug & ' dx,dy,dz= ',1p3e22.14)
cbug write ( 3, 9901)
cbug write ( 3, 9902) (n, ax(n), ay(n), az(n), bx(n), by(n), bz(n),
cbug & cx(n), cy(n), cz(n), dx(n), dy(n), dz(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)
c.... Loop over subsets of data.
110 ns = n2 - n1 + 1
c.... Find the unit vector "ub" normal to plane "a".
call aptvunb (bx(n1), by(n1), bz(n1), ns, 0.,
& ubx, uby, ubz, vlenb, nerr)
c.... Find the unit vector "ud" normal to plane "c".
call aptvunb (dx(n1), dy(n1), dz(n1), ns, 0.,
& udx, udy, udz, vlend, nerr)
c.... Find the unit vector "u" parallel to the line of intersection.
call aptvxun (ubx, uby, ubz, udx, udy, udz, ns, tol,
& ux, uy, uz, sinth, nerr)
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptplpl unit normal vectors to planes:' /
cbug & (i3 /
cbug & ' ubx,uby,ubz=',1p3e22.14 /
cbug & ' udx,udy,udz=',1p3e22.14 /
cbug & ' unit vector along intersection:' /
cbug & ' ux,uy,uz= ',1p3e22.14))
cbug write ( 3, 9903) (n, ubx(n), uby(n), ubz(n),
cbug & udx(n), udy(n), udz(n), ux(n), uy(n), uz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Find the unit vector "uae" in plane "a" perpendicular to the
c.... line of intersection "ef".
call aptvxun (ubx, uby, ubz, ux, uy, uz, ns, tol,
& uaex, uaey, uaez, vlen, nerr)
c.... Find the unit vector "ucf" in plane "a" perpendicular to the
c.... line of intersection "ef".
call aptvxun (udx, udy, udz, ux, uy, uz, ns, tol,
& ucfx, ucfy, ucfz, vlen, nerr)
c.... Find the scalar (dot) product costh of vectors "uae" and "ucf".
c.... Equal to the scalar (dot) product of vectors "ub" and "ud".
call aptvdot (ubx, uby, ubz, udx, udy, udz, ns, tol,
& costh, nerr)
c.... Find the vector vac between the planar points "a" and "c".
call aptvdis (ax(n1), ay(n1), az(n1), cx(n1), cy(n1), cz(n1),
& ns, tol, vacx, vacy, vacz, dac, nerr)
c.... Find the scalar (dot) product dacae of vectors "vac" and "uae".
call aptvdot (vacx, vacy, vacz, uaex, uaey, uaez, ns, tol,
& dacae, nerr)
c.... Find the scalar (dot) product daccf of vectors "vac" and "ucf".
call aptvdot (vacx, vacy, vacz, ucfx, ucfy, ucfz, ns, tol,
& daccf, nerr)
c.... Find the distances dae from "a" to "e", and dcf from "c" to "f".
c.... sinth**2 = 1.0 - costh**2.
c---- Loop over subset of data.
do 120 n = 1, ns
dae(n) = (dacae(n) - costh(n) * daccf(n)) / (sinth(n)**2 + fuz)
dcf(n) = (daccf(n) - costh(n) * dacae(n)) / (sinth(n)**2 + fuz)
if ((vlenb(n) .le. tol) .or.
& (vlend(n) .le. tol)) then
dae(n) = 0.0
dcf(n) = 0.0
endif
c---- End of loop over subset of data.
120 continue
c.... Find the coordinates of points "e" and "f".
call aptvadd (ax(n1), ay(n1), az(n1), 1., dae,
& uaex, uaey, uaez, ns, tol,
& ex(n1), ey(n1), ez(n1), vlen, nerr)
call aptvadd (cx(n1), cy(n1), cz(n1), -1., dcf,
& ucfx, ucfy, ucfz, ns, tol,
& fx(n1), fy(n1), fz(n1), vlen, nerr)
c.... Find the component of the distance "ac" in direction "b".
c.... (The distance between planes, if they are parallel.)
call aptvdot (vacx, vacy, vacz, ubx, uby, ubz, ns, tol,
& dpmin(n1), nerr)
c.... Find the unit vector in the direction of vector "ac".
c.... (The direction of unit vector "u", if the planes are coincident.)
call aptvuna (vacx, vacy, vacz, ns, 0., dac, nerr)
c.... Test for parallel planes, or bad input vectors "b" and/or "d".
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
if (sinth(n) .gt. 0.0) then
ipar(nn) = 0
else
ipar(nn) = 1
endif
if (vlenb(n) .le. tol) then
ipar(nn) = 2
endif
if (vlend(n) .le. tol) then
ipar(nn) = 3
endif
if ((vlenb(n) .le. tol) .or.
& (vlend(n) .le. tol)) then
ipar(nn) = 4
endif
if (ipar(nn) .ne. 1) then
dpmin(nn) = 0.0
endif
if ((dpmin(nn) .eq. 0.0) .and.
& (ipar(nn) .eq. 1)) then
itrun(nn) = 1
else
itrun(nn) = 0
endif
if (ipar(nn) .ne. 0) then
ux(nn) = vacx(n)
uy(nn) = vacy(n)
uz(nn) = vacz(n)
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 9904 format (/ 'aptplpl results:' /
cbug & (i3,' dpmin= ',1pe22.14,' ipar=',i2,' itrun=',i2 /
cbug & ' ex,ey,ez= ',1p3e22.14 /
cbug & ' fx,fy,fz= ',1p3e22.14 /
cbug & ' ux,uy,uz= ',1p3e22.14))
cbug write ( 3, 9904) (n, dpmin(n), ipar(n), itrun(n),
cbug & ex(n), ey(n), ez(n), fx(n), fy(n), fz(n), ux(n), uy(n), uz(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptplpl. (+1 line.)
end
UCRL-WEB-209832