subroutine aptlnln (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, np, tol, dpmin, fracab, fraccd,
& ex, ey, ez, fx, fy, fz, itrun, ipar, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTLNLN
c
c call aptlnln (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, dpmin, fracab, fraccd,
c & ex, ey, ez, fx, fy, fz, itrun, ipar, nerr)
c
c Version: aptlnln Updated 1994 January 5 17:30.
c aptlnln 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 minimum
c distance dpmin between the line through the points
c a = (ax, ay, az) and b = (bx, by, bz), and the line through
c the points c = (cx, cy, cz) and d = (dx, dy, dz), and the
c point e = (ex, ey, ez) on line "ab", and the point
c f = (fx, fy, fz) on line "cd", at which the minimum distance
c dpmin occurs. If dpmin is smaller than the estimated error in
c its calculation, it will be truncated to zero, and itrun = 1
c will be returned.
c The fractional distance fracab of point "e" along line "ab",
c and the fractional distance fraccd of point "f" along line
c "cd", are also returned.
c If the lines are parallel, ipar = 1 will be returned.
c If a line segment is too short, ipar = 2, 3 or 4 will be
c returned.
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c 1994 January 5. Added truncation of coordinates of points e
c and f.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, tol.
c
c Output: dpmin, fracab, fraccd, ex, ey, ez, fx, fy, fz,
c itrun, ipar, nerr.
c
c Calls: aptvdis, aptvdot, aptvuna
c
c
c Glossary:
c
c ax,ay,az Input The first point on line segment "ab". Size np.
c
c bx,by,bz Input The second point on line segment "ab". Size np.
c
c cx,cy,cz Input The first point on line segment "cd". Size np.
c
c dpmin Output Minimum distance from line "ab" to line "cd".
c Distance from e = (ex, ey, ez) to f = (fx, fy, fz).
c Truncated to zero if less than the estimated error
c in its calculation. See tol. Size np.
c
c dx,dy,dz Input The second point on line segment "cd". Size np.
c
c ex,ey,ez Output The x, y, z coordinates of the point on line "ab"
c nearest line "cd". Size np.
c
c fracab Output Fractional distance of "e" along line "ab". Size np.
c Meaningless if ipar = 2 or 4.
c
c fraccd Output Fractional distance of "f" along line "cd". Size np.
c Meaningless if ipar = 3 or 4.
c
c fx,fy,fz Output The x, y, z coordinates of the point on line "cd"
c nearest line "ab". Size np.
c
c ipar Output 0 if lines are not parallel. 1 if they are, and
c points "e" and "f" can be moved arbitrarily by
c equal distances along the lines. Size np.
c 2 if line segment "ab" is too short.
c 3 if line segment "cd" is too short.
c 4 if "ab" and "cd" are both too short.
c
c itrun Output 0 if dpmin not truncated to zero, based on tol.
c 1 if dpmin is truncated to zero, based on tol.
c Size np.
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 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 line "ab" to line "cd".
dimension dpmin (1)
c---- Coordinate x of point "d".
dimension dx (1)
c---- Coordinate y of point "d".
dimension dy (1)
c---- Coordinate z of point "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---- Fractional distance of "e" on "ab".
dimension fracab (1)
c---- Fractional distance of "f" on "cd".
dimension fraccd (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---- 1 if lines parallel, 2 if bad.
dimension ipar (1)
c---- 1 if dpmin truncated to zero.
dimension itrun (1)
c.... Local variables.
c---- Component of "ac" in direction "ab".
common /laptlnln/ acdotab (64)
c---- Component of "ac" in direction "cd".
common /laptlnln/ acdotcd (64)
c---- Component x of vector "ac".
common /laptlnln/ acx (64)
c---- Component y of vector "ac".
common /laptlnln/ acy (64)
c---- Component z of vector "ac".
common /laptlnln/ acz (64)
c---- Cosine of angle between "ab" and "cd".
common /laptlnln/ costh (64)
c---- Distance from "a" to "b".
common /laptlnln/ distab (64)
c---- Distance from "a" to "e".
common /laptlnln/ distae
c---- Distance from "c" to "d".
common /laptlnln/ distcd (64)
c---- Distance from "c" to "f".
common /laptlnln/ distcf
c---- Value of dpmin before truncating.
common /laptlnln/ dpmins (64)
c---- Component x of vector "ef".
common /laptlnln/ efx (64)
c---- Component y of vector "ef".
common /laptlnln/ efy (64)
c---- Component z of vector "ef".
common /laptlnln/ efz (64)
c---- Estimated error in sinth2.
common /laptlnln/ errs
c---- Estimated error in ex or fx.
common /laptlnln/ errx
c---- Estimated error in ey or fy.
common /laptlnln/ erry
c---- Estimated error in ez or fz.
common /laptlnln/ errz
c---- A very small number.
common /laptlnln/ fuz
c---- Index in arrays.
common /laptlnln/ n
c---- First index of subset of data.
common /laptlnln/ n1
c---- Last index of subset of data.
common /laptlnln/ n2
c---- Index in external array.
common /laptlnln/ nn
c---- Size of current subset of data.
common /laptlnln/ ns
c---- Sine**2 of angle between "ab" and "cd"
common /laptlnln/ sinth2
c---- Component x of unit vector along "ab".
common /laptlnln/ uabx (64)
c---- Component y of unit vector along "ab".
common /laptlnln/ uaby (64)
c---- Component z of unit vector along "ab".
common /laptlnln/ uabz (64)
c---- Component x of unit vector along "cd".
common /laptlnln/ ucdx (64)
c---- Component y of unit vector along "cd".
common /laptlnln/ ucdy (64)
c---- Component z of unit vector along "cd".
common /laptlnln/ ucdz (64)
c---- Magnitude of a vector.
common /laptlnln/ vlen (64)
c---- Maqnitude of a vector.
common /laptlnln/ vlene
c---- Maqnitude of a vector.
common /laptlnln/ vlenf
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptlnln finding distance from the line through:')
cbug 9902 format (i3,' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' to the line through:' /
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 "uab" along line "ab".
call aptvdis (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), ns,
& tol, uabx, uaby, uabz, distab, nerr)
call aptvuna (uabx, uaby, uabz, ns, 0., vlen, nerr)
c.... Find the unit vector "ucd" along line "cd".
call aptvdis (cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1), ns,
& tol, ucdx, ucdy, ucdz, distcd, nerr)
call aptvuna (ucdx, ucdy, ucdz, ns, 0., vlen, nerr)
c.... Find the vector ac from point "a" to point "c".
call aptvdis (ax(n1), ay(n1), az(n1), cx(n1), cy(n1), cz(n1), ns,
& tol, acx, acy, acz, vlen, nerr)
c.... Find the components of vector "ac" in the directions of the lines
c.... (the scalar products of "ac" with "uab" and "ucd").
call aptvdot (acx, acy, acz, uabx, uaby, uabz, ns, tol,
& acdotab, nerr)
call aptvdot (acx, acy, acz, ucdx, ucdy, ucdz, ns, tol,
& acdotcd, nerr)
c.... Find the scalar (dot) product of the two unit direction vectors.
call aptvdot (uabx, uaby, uabz, ucdx, ucdy, ucdz, ns, tol,
& costh, nerr)
c.... Find the points of nearest approach on each line, "e" and "f",
c.... and their fractional distances along lines "ab" and "cd".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
sinth2 = (1.0 - costh(n)) * (1.0 + costh(n))
if ((distab(n) * distcd(n)) .eq. 0.0) then
sinth2 = 0.0
endif
errs = tol * (1.0 + costh(n)**2)
if (abs (sinth2) .ge. errs) then
ipar(nn) = 0
else
ipar(nn) = 1
endif
if (distab(n) .eq. 0.0) then
ipar(nn) = 2
endif
if (distcd(n) .eq. 0.0) then
ipar(nn) = ipar(nn) + 2
endif
distae = (acdotab(n) - acdotcd(n) * costh(n)) /
& (sinth2 + fuz)
distcf = (acdotab(n) * costh(n) - acdotcd(n)) /
& (sinth2 + fuz)
if (abs (sinth2) .le. errs) then
distae = 0.5 * acdotab(n)
distcf = -0.5 * acdotcd(n)
endif
ex(nn) = ax(nn) + distae * uabx(n)
ey(nn) = ay(nn) + distae * uaby(n)
ez(nn) = az(nn) + distae * uabz(n)
errx = tol * (abs (ax(nn)) + abs (bx(nn)))
erry = tol * (abs (ay(nn)) + abs (by(nn)))
errz = tol * (abs (az(nn)) + abs (bz(nn)))
if (abs (ex(nn)) .le. errx) ex(nn) = 0.0
if (abs (ey(nn)) .le. erry) ey(nn) = 0.0
if (abs (ez(nn)) .le. errz) ez(nn) = 0.0
fx(nn) = cx(nn) + distcf * ucdx(n)
fy(nn) = cy(nn) + distcf * ucdy(n)
fz(nn) = cz(nn) + distcf * ucdz(n)
errx = tol * (abs (cx(nn)) + abs (dx(nn)))
erry = tol * (abs (cy(nn)) + abs (dy(nn)))
errz = tol * (abs (cz(nn)) + abs (dz(nn)))
if (abs (fx(nn)) .le. errx) fx(nn) = 0.0
if (abs (fy(nn)) .le. erry) fy(nn) = 0.0
if (abs (fz(nn)) .le. errz) fz(nn) = 0.0
fracab(nn) = distae / (distab(n) + fuz)
fraccd(nn) = distcf / (distcd(n) + fuz)
c---- End of loop over subset of data.
120 continue
c.... Find the minimum distance between the lines.
call aptvdis (ex(n1), ey(n1), ez(n1), fx(n1), fy(n1), fz(n1), ns,
& 0.0, efx, efy, efz, dpmin(n1), nerr)
c.... See if truncation to zero is allowed.
c---- Truncation is allowed.
if (tol .gt. 0.0) then
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
dpmins(n) = dpmin(nn)
c---- End of loop over subset of data.
130 continue
call aptvdis (ex(n1), ey(n1), ez(n1), fx(n1), fy(n1), fz(n1),
& ns, tol, efx, efy, efz, dpmin(n1), nerr)
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
if (dpmin(nn) .eq. dpmins(n)) then
itrun(nn) = 0
else
itrun(nn) = 1
endif
c---- End of loop over subset of data.
140 continue
c---- Tested tol.
endif
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 9903 format (/ 'aptlnln results:' /
cbug & (i3,' dpmin= ',1pe22.14 /
cbug & ' fracab= ',1pe22.14,' fraccd=',1pe22.14 /
cbug & ' itrun=',i2,' ipar=',i2 /
cbug & ' ex,ey,ez=',1p3e22.14 /
cbug & ' fx,fy,fz=',1p3e22.14))
cbug write ( 3, 9903) (n, dpmin(n), fracab(n), fraccd(n), itrun(n),
cbug & ipar(n), ex(n), ey(n), ez(n), fx(n), fy(n), fz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptlnln. (+1 line.)
end
UCRL-WEB-209832