subroutine aptptln (px, py, pz, ax, ay, az, bx, by, bz, np, tol,
& noptfd, dpmin, fdmin, cx, cy, cz,
& nlim, itrun, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPTLN
c
c call aptptln (px, py, pz, ax, ay, az, bx, by, bz, np, tol,
c & noptfd, dpmin, fdmin, cx, cy, cz, nlim, itrun, nerr)
c
c Version: aptptln Updated 1990 November 29 10:50.
c aptptln 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 distance
c dpmin from the point p = (px, py, pz) to the straight line
c through the points a = (ax, ay, az) and b = (bx, by, bz),
c to find the coordinates c = (cx, cy, cz) of the proximate point
c on line "ab", and to find the fractional distance fdmin of that
c point along the line segment "ab". The value of
c dpmin will be truncated to zero if less than the estimated error
c in its calculation, based on tol, and if so itrun = 1 will be
c returned. Flag nerr indicates any input error.
c
c Option noptfd allows the line "ab" to be treated as a finite
c segment, by limiting the range of fdmin. Flag nlim indicates
c when such limitation has beem imposed.
c
c If the points "a" and "b" coincide, based on tol, dpmin will be
c the distance from point "a" to point "p", and if dpmin is not
c zero, itrun will be -1.
c
c History: 1990 February 12 16:00. Fixed bug affecting fdmin when
c np is greater than 64.
c
c Input: px, py, pz, ax, ay, az, bx, by, bz, np, tol, noptfd.
c
c Output: dpmin, fdmin, cx, cy, cz, nlim, itrun, nerr.
c
c Calls: aptfdav, aptvadd, aptvdis, aptvdot
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a" on line "ab".
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of point "b" on line "ab".
c Size np.
c
c cx,cy,cz Output The x, y, z coordinates of the point on the line "ab"
c nearest point "p". Point "p" if dpmin = 0. Size np.
c
c dpmin Output Distance from point "p" to the line "ab". Size np.
c Truncated to zero if less than the estimated error
c in its calculation, based on tol, and if so,
c itrun = 1 will be returned. Perpendicular distance,
c unless itrun = -1, or nlim = 2.
c
c fdmin Output Fractional distance between point "a" and point "b"
c of the point "c". Size np.
c
c itrun Output 0 if dpmin is not truncated to zero.
c 1 if dpmin is truncated to zero, when less than its
c estimated error, based on tol. Size np.
c -1 if dpmin is not zero, and points "a" and "b"
c coincide, based on tol.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if noptfd is not 0, 1, or 2.
c
c nlim Output 0 if no limit imposed on fdmin, 1 if the limit of
c noptfd = 1 is imposed, 2 if the limit of noptfd = 2
c is imposed. Size np.
c If the latter limit is imposed, dpmin will be the
c distance from point "p" to the nearest end of the
c line segment.
c
c noptfd Input Option to limit range of fdmin: 0 for no limit;
c 1 to increase fdmin to tol, if in the range from
c -tol to tol, and decrease fdmin to 1.0 - tol, if
c in the range from 1.0 - tol to 1.0 + tol; and
c 2 to impose the limits for noptfd = 1, and then
c limit fdmin to the range from 0.0 to 1.0.
c
c np Input Size of arrays px, py, pz, ax, ay, az, bx, by, bz,
c cx, cy, cz, fdmin, itrun.
c
c px,py,pz Input The x, y, z coordinates of point "p". 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 line "ab" to point "p".
dimension dpmin (1)
c---- Fractional distance of "c" on "ab"
dimension fdmin (1)
c---- Indicates dpmin zero, if 1.
dimension itrun (1)
c---- Indicates fdmin limited, if not 0.
dimension nlim (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.... Local variables.
c---- Component x of vector "ab".
common /laptptln/ abx (64)
c---- Component y of vector "ab".
common /laptptln/ aby (64)
c---- Component z of vector "ab".
common /laptptln/ abz (64)
c---- Component x of vector "ap".
common /laptptln/ apx (64)
c---- Component y of vector "ap".
common /laptptln/ apy (64)
c---- Component z of vector "ap".
common /laptptln/ apz (64)
c---- Magnitude of vector "ab".
common /laptptln/ dab (64)
c---- Magnitude of vector "ap".
common /laptptln/ dap (64)
c---- Component x of vector "ab" x "ap".
common /laptptln/ dx (64)
c---- Component y of vector "ab" x "ap".
common /laptptln/ dy (64)
c---- Component z of vector "ab" x "ap".
common /laptptln/ dz (64)
c---- A very small number.
common /laptptln/ fuz
c---- Index in arrays.
common /laptptln/ n
c---- First index of subset of data.
common /laptptln/ n1
c---- Last index of subset of data.
common /laptptln/ n2
c---- Index in external array.
common /laptptln/ nn
c---- Size of current subset of data.
common /laptptln/ ns
c---- Magnitude of a vector.
common /laptptln/ vlen (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptptln finding distance from points',
cbug & ' (noptfd=',i2,' tol=',1pe13.5,'):')
cbug 9902 format (
cbug & i3,' px,py,pz=',1p3e22.14 /
cbug & ' to the line segment with endpoints:' /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14)
cbug write ( 3, 9901) noptfd, tol
cbug write ( 3, 9902) (n, px(n), py(n), pz(n), ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(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
if ((noptfd .lt. 0) .or. (noptfd .gt. 2)) then
nerr = 2
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 vector "ab" along the line segment, and its length, dab.
call aptvdis (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), ns,
& tol, abx, aby, abz, dab, 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, dap, nerr)
c.... Find the scalar product of vectors "ab" and "ap".
call aptvdot (abx, aby, abz, apx, apy, apz, ns, tol,
& fdmin(n1), nerr)
c.... Find the fractional distance fdmin, along the extended line "ab",
c.... of the point "c" nearest point "p".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
fdmin(nn) = fdmin(nn) / (dab(n)**2 + fuz)
c---- End of loop over subset of data.
120 continue
c.... Test fdmin for limits if noptfd = 1 or 2.
c---- Impose limits on fdmin.
if (noptfd .ne. 0) then
call aptfdav (fdmin(n1), ns, noptfd, tol, nlim(n1), nerr)
c---- Tested noptfd.
endif
c.... Find the point "c" on the line "ab", nearest to point "p".
call aptvadd (ax(n1), ay(n1), az(n1), 1., fdmin(n1),
& abx, aby, abz, ns, tol,
& cx(n1), cy(n1), cz(n1), vlen, nerr)
c.... Find the distance dpmin from point "p" to point "c".
call aptvdis (px(n1), py(n1), pz(n1), cx(n1), cy(n1), cz(n1),
& ns, tol, dx, dy, dz, dpmin(n1), nerr)
c.... Recalculate dpmin if points "a" and "b" coincide.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
if (dab(n) .le. 0.0) then
dpmin(nn) = dap(n)
itrun(nn) = -1
else
itrun(nn) = 0
endif
c---- End of loop over subset of data.
130 continue
c.... Reduce dpmin to zero if it is small, relative to the line length,
c.... or less than the estimated error in its calculation.
c.... (May have already been done in aptvdis.)
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
if (abs (dpmin(nn)) .le. tol * dab(n)) then
itrun(nn) = itrun(nn) + 1
dpmin(nn) = 0.0
cx(nn) = px(nn)
cy(nn) = py(nn)
cz(nn) = pz(nn)
endif
c---- End of loop over subset of data.
140 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 9903 format (/ 'aptptln results:' /
cbug & (i3,' itrun=',i2,' nlim=',i2 /
cbug & ' dpmin= ',1pe22.14,' fdmin=',1pe22.14 /
cbug & ' cx,cy,cz=',1p3e22.14))
cbug write ( 3, 9903) (n, itrun(n), nlim(n), dpmin(n), fdmin(n),
cbug & cx(n), cy(n), cz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptptln. (+1 line.)
end
UCRL-WEB-209832