subroutine aptptlc (pu, pv, au, av, bu, bv, np, tol, noptfd,
& dpmin, fdmin, nlim, itrun, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPTLC
c
c call aptptlc (pu, pv, au, av, bu, bv, np, tol, noptfd,
c & dpmin, fdmin, nlim, itrun, nerr)
c
c Version: aptptlc Updated 1990 November 29 10:50.
c aptptlc Originated 1989 December 29 16:40.
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 = (pu, pv) to the straight line through
c the points a = (au, av) and b = (bu, bv). All points are in
c the uv plane. Directions u, v and w are orthogonal.
c
c Option noptfd allows the fractional distance fdmin of the
c proximal point along line "ab" to be calculated, and allows
c the range of fdmin to be limited. Flag nlim indicates
c when such limitation has beem imposed. If nlim = 2, dpmin is
c the distance from point "p" to the nearest end of line
c segment "ab".
c
c The value of dpmin will be truncated to zero if less than the
c estimated error in its calculation, based on tol, and if so,
c itrun = 1 will be returned.
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 Flag nerr indicates any input error.
c
c History: 1990 February 12 15:20. Added input argument noptfd,
c optional output arguments fdmin, nlim.
c
c Input: pu, pv, au, av, bu, bv, np, tol, noptfd.
c
c Output: dpmin, fdmin, nlim, itrun, nerr.
c
c Calls: aptfdav, aptvaxc, aptvdic, aptvdoc
c
c
c Glossary:
c
c au, av Input The u and v coordinates of point "a" on line "ab".
c Size np. All points are in the uv plane.
c
c bu, bv Input The u and v coordinates of point "b" on line "ab".
c 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.
c If points "a" and "b" are coincident, dpmin is the
c distance from point "p" to point "a", and itrun = -1
c is returned, unless dpmin = 0.0.
c If noptfd is 2, and fdmin is initially not in the
c range from 0.0 to 1.0, dpmin is the distance from
c point "p" to the nearest of points "a" and "b", and
c nlim = 2 is returned.
c The value of dpmin is positive, if point "p" is
c to the left of the vector "ab" in the uv plane.
c
c fdmin Output Fractional distance between point "a" and point "b"
c of the point nearest point "p". See noptfd.
c Size np, if noptfd is not -1.
c
c itrun Output 0 if dpmin is not truncated to zero. Size np.
c 1 if dpmin is truncated to zero, when less than its
c estimated error, based on tol.
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 -1, 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, if noptfd is not -1.
c If 2, dpmin is the distance from point "p" to the
c nearest end of line segment "ab".
c
c noptfd Input Option to limit range of fdmin: 0 for no limit;
c -1 to not calculate fdmin or nlim.
c 0 to find fdmin, but impose no limits.
c 1 to increase fdmin to tol, if in the range from
c -tol to tol, and decrease fdmin to 1.0 - tol,
c if in the range from 1.0 - tol to 1.0 + tol.
c 2 to impose the limits for noptfd = 1, and then
c limit fdmin to the range from 0.0 to 1.0, and
c adjust the magnitude, but not the sign of dpmin
c if the later limit is imposed.
c
c np Input Size of arrays pu, pv, au, av, bu, bv, itrun.
c If noptfd is not -1, the size of arrays fdmin, nlim.
c
c pu, pv Input The u and v 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 u of point "a".
dimension au (1)
c---- Coordinate v of point "a".
dimension av (1)
c---- Coordinate u of point "b".
dimension bu (1)
c---- Coordinate v of point "b".
dimension bv (1)
c---- Distance from line "ab" to point "p".
dimension dpmin (1)
c---- Fractional distance along "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 u of point "p".
dimension pu (1)
c---- Coordinate v of point "p".
dimension pv (1)
c.... Local variables.
c---- Component u of vector "ab".
common /laptptlc/ abu (64)
c---- Component v of vector "ab".
common /laptptlc/ abv (64)
c---- Component u of vector "ap".
common /laptptlc/ apu (64)
c---- Component v of vector "ap".
common /laptptlc/ apv (64)
c---- Component u of vector "bp".
common /laptptlc/ bpu (64)
c---- Component v of vector "bp".
common /laptptlc/ bpv (64)
c---- Magnitude of vector "ab".
common /laptptlc/ dab (64)
c---- Magnitude of vector "ap".
common /laptptlc/ dap (64)
c---- Magnitude of vector "bp".
common /laptptlc/ dbp (64)
c---- Signed value of dap.
common /laptptlc/ dpmina
c---- Signed value of dbp.
common /laptptlc/ dpminb
c---- Component w of vector "ab" x "ap".
common /laptptlc/ dw (64)
c---- A very small number.
common /laptptlc/ fuz
c---- Index in arrays.
common /laptptlc/ n
c---- First index of subset of data.
common /laptptlc/ n1
c---- Last index of subset of data.
common /laptptlc/ n2
c---- Index in external array.
common /laptptlc/ nn
c---- Size of current subset of data.
common /laptptlc/ ns
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptptlc finding distance from points',
cbug & ' (noptfd=',i2,' tol=',1pe13.5,'):')
cbug 9902 format (
cbug & i3,' pu,pv=',1p2e22.14 /
cbug & ' to the line segment with endpoints:' /
cbug & ' au,av=',1p2e22.14 /
cbug & ' bu,bv=',1p2e22.14)
cbug write ( 3, 9901) noptfd, tol
cbug write ( 3, 9902) (n, pu(n), pv(n), au(n), av(n),
cbug & bu(n), bv(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. -1) .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 aptvdic (au(n1), av(n1), bu(n1), bv(n1), ns, tol,
& abu, abv, dab, nerr)
c.... Find the vector "ap" from point "a" to point "p", and its length, dap.
call aptvdic (au(n1), av(n1), pu(n1), pv(n1), ns, tol,
& apu, apv, dap, nerr)
c.... Find the parallelogram area for vectors "ab", "ap".
call aptvaxc (abu, abv, apu, apv, ns, tol, dw, nerr)
c.... Find the minimum distance dpmin from point "p" to line "ab".
c.... Recalculate dpmin if points "a" and "b" coincide.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
dpmin(nn) = dw(n) / (dab(n) + fuz)
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.
120 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 aptvdic.)
c---- Loop over subset of data.
do 130 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
endif
c---- End of loop over subset of data.
130 continue
c.... Find fdmin if noptfd is not -1.
c---- Find fdmin.
if (noptfd .ne. -1) then
c.... Find the scalar product of vectors "ab" and "ap".
call aptvdoc (abu, abv, apu, apv, ns, tol, fdmin(n1), nerr)
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
fdmin(nn) = fdmin(nn) / (dab(n)**2 + fuz)
c---- End of loop over subset of data.
140 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 for noptfd = 1 or 2.
endif
c---- Adjust dpmin when nlim = 2.
if (noptfd .eq. 2) then
c.... Find the vector "bp" from point "b" to point "p",
c.... and its length, dbp.
call aptvdic (bu(n1), bv(n1), pu(n1), pv(n1), ns, tol,
& bpu, bpv, dbp, nerr)
c.... Adjust dpmin if nlim = 2.
c---- Loop over subset of data.
do 150 n = 1, ns
nn = n + n1 - 1
dpmina = sign (dap(n), dpmin(nn))
dpminb = sign (dbp(n), dpmin(nn))
if ((nlim(nn) .eq. 2) .and. (fdmin(nn) .eq. 0.0)) then
dpmin(nn) = dpmina
endif
if ((nlim(nn) .eq. 2) .and. (fdmin(nn) .eq. 1.0)) then
dpmin(nn) = dpminb
endif
if ((itrun(nn) .eq. 1) .and. (dpmin(nn) .ne. 0.0)) then
itrun(nn) = 0
endif
c---- End of loop over subset of data.
150 continue
c---- Tested for noptfd = 2.
endif
c---- Tested for noptfd .ne. -1.
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 (/ 'aptptlc results:')
cbug 9904 format (i3,' itrun=',i2,' dpmin=',1pe22.14)
cbug 9905 format (' nlim= ',i2,' fdmin=',1pe22.14)
cbug write ( 3, 9903)
cbug do 160 n = 1, np
cbug write ( 3, 9904) n, itrun(n), dpmin(n)
cbug if (noptfd .ne. -1) then
cbug write ( 3, 9905) nlim(n), fdmin(n)
cbug endif
cbug 160 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptptlc. (+1 line.)
end
UCRL-WEB-209832