subroutine aptlnlc (au, av, bu, bv, cu, cv, du, dv, np, tol,
& dpmin, fracab, fraccd, eu, ev, ipar, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTLNLC
c
c call aptlnlc (au, av, bu, bv, cu, cv, du, dv, np, tol,
c & dpmin, fracab, fraccd, eu, ev, ipar, nerr)
c
c Version: aptlnlc Updated 1990 November 28 10:00.
c aptlnlc Originated 1990 January 11 11: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 intersection
c e = (eu, ev), if any, between the straight line between points
c a = (au, av) and b = (bu, bv), and the straight line between
c points c = (cu, cv), and d = (du, dv), all in the uv plane,
c where u, v and w are orthogonal directions. The fractional
c distance fracab of point "e" along line "ab", and the fractional
c distance fraccd of point "e" along line "cd" are also returned.
c If the lines are parallel, ipar = 1 will be returned, and the
c distance dpmin between the lines will be returned. If dpmin is
c smaller than the estimated error in its calculation, it will be
c truncated to zero. If a line segment is too short, ipar = 2, 3
c or 4 will be returned. Flag nerr will be 1 if np is not
c positive.
c
c Input: au, av, bu, bv, cu, cv, du, dv, np, tol.
c
c Output: dpmin, fracab, fraccd, eu, ev, ipar, nerr.
c
c Calls: aptvdic, aptaxc
c
c
c Glossary:
c
c au, av Input The u and v coordinates of the first point on line
c segment "ab" in the uv plane. Size np.
c
c bu, bv Input The u and v coordinates of the second point on line
c segment "ab" in the uv plane. Size np.
c
c cu, cv Input The u and v coordinates of the first point on line
c segment "cd" in the uv plane. Size np.
c
c dpmin Output Distance from line "ab" to line "cd", when they are
c parallel (ipar = 1).
c Truncated to zero if less than the estimated error
c in its calculation. See tol. Size np.
c
c du, dv Input The u and v coordinates of the second point on line
c segment "cd" in the uv plane. Size np.
c
c eu, ev Output The u and v coordinates of the intersection of lines
c "ab" and "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 "e" along line "cd". Size np.
c Meaningless if ipar = 3 or 4.
c
c ipar Output 0 if lines are not parallel. 1 if they are. 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 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
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---- Coordinate u of point "c".
dimension cu (1)
c---- Coordinate v of point "c".
dimension cv (1)
c---- Distance from line "ab" to line "cd".
dimension dpmin (1)
c---- Coordinate u of point "d".
dimension du (1)
c---- Coordinate v of point "d".
dimension dv (1)
c---- Coordinate u of point "e".
dimension eu (1)
c---- Coordinate v of point "e".
dimension ev (1)
c---- Fractional distance of "e" on "ab".
dimension fracab (1)
c---- Fractional distance of "e" on "cd".
dimension fraccd (1)
c---- 1 if lines parallel, 2 if bad.
dimension ipar (1)
c.... Local variables.
c---- Cross product of "ab" and "cd".
common /laptlnlc/ abcd (64)
c---- Component u of unit vector along "ab".
common /laptlnlc/ abu (64)
c---- Component v of unit vector along "ab".
common /laptlnlc/ abv (64)
c---- Cross product of "ac" and "ab".
common /laptlnlc/ acab (64)
c---- Cross product of "ac" and "cd".
common /laptlnlc/ accd (64)
c---- Component u of vector "ac".
common /laptlnlc/ acu (64)
c---- Component v of vector "ac".
common /laptlnlc/ acv (64)
c---- Component u of unit vector along "cd".
common /laptlnlc/ cdu (64)
c---- Component v of unit vector along "cd".
common /laptlnlc/ cdv (64)
c---- Distance from "a" to "b".
common /laptlnlc/ distab (64)
c---- Distance from "a" to "c".
common /laptlnlc/ distac (64)
c---- Distance from "c" to "d".
common /laptlnlc/ distcd (64)
c---- A very small number.
common /laptlnlc/ fuz
c---- Distance between parallel lines.
common /laptlnlc/ dpar
c---- Distance between "c" and "ab".
common /laptlnlc/ dparab
c---- Distance between "a" and "cd".
common /laptlnlc/ dparcd
c---- Coordinate u of intersection on "ab".
common /laptlnlc/ eabu
c---- Coordinate v of intersection on "ab".
common /laptlnlc/ eabv
c---- Coordinate u of intersection on "cd".
common /laptlnlc/ ecdu
c---- Coordinate v of intersection on "cd".
common /laptlnlc/ ecdv
c---- Index in external array.
common /laptlnlc/ n
c---- First index of subset of data.
common /laptlnlc/ n1
c---- Last index of subset of data.
common /laptlnlc/ n2
c---- Index in external array.
common /laptlnlc/ nn
c---- Size of current subset of data.
common /laptlnlc/ ns
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptlnlc finding intersection of the line through' /
cbug & ' (np=',i3,' tol=',1pe13.5,'):')
cbug 9902 format (i3,' au,av= ',1p2e22.14 /
cbug & ' bu,bv= ',1p2e22.14 /
cbug & ' with the line through:' /
cbug & ' cu,cv= ',1p2e22.14 /
cbug & ' du,dv= ',1p2e22.14)
cbug write ( 3, 9901) np, tol
cbug write ( 3, 9902) (n, au(n), av(n), bu(n), bv(n),
cbug & cu(n), cv(n), du(n), dv(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 2-D vector "ab" from point "a" to point "b".
call aptvdic (au(n1), av(n1), bu(n1), bv(n1), ns, tol,
& abu, abv, distab, nerr)
c.... Find the 2-D vector "cd" from point "c" to point "d".
call aptvdic (cu(n1), cv(n1), du(n1), dv(n1), ns, tol,
& cdu, cdv, distcd, nerr)
c.... Find the 2-D vector "ac" from point "a" to point "c".
call aptvdic (au(n1), av(n1), cu(n1), cv(n1), ns, tol,
& acu, acv, distac, nerr)
c.... Find the vector product of vectors "ab" and "cd".
call aptvaxc (abu, abv, cdu, cdv, ns, tol, abcd, nerr)
c.... Find the vector product of vectors "ac" and "ab".
call aptvaxc (acu, acv, abu, abv, ns, tol, acab, nerr)
c.... Find the vector product of vectors "ac" and "cd".
call aptvaxc (acu, acv, cdu, cdv, ns, tol, accd, nerr)
c.... Find the fractional distances of the intersection "e" along lines
c.... "ab" and "cd". Set dpmin to zero.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
dpmin(nn) = 0.0
fracab(nn) = accd(n) / (abcd(n) + fuz)
fraccd(nn) = acab(n) / (abcd(n) + fuz)
c---- End of loop over subset of data.
120 continue
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
eabu = au(nn) + fracab(nn) * abu(n)
eabv = av(nn) + fracab(nn) * abv(n)
ecdu = cu(nn) + fraccd(nn) * cdu(n)
ecdv = cv(nn) + fraccd(nn) * cdv(n)
c.... See if on line "ab" or on line "cd". Use shorter if on both.
if (abs (fracab(nn)) * distab(n) .lt.
& abs (fraccd(nn)) * distcd(n)) then
eu(nn) = eabu
ev(nn) = eabv
else
eu(nn) = ecdu
ev(nn) = ecdv
endif
c---- End of loop over subset of data.
130 continue
c.... If the lines are parallel, find the distance between them.
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
if (abcd(n) .ne. 0.0) 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
dparab = abs (acab(n) / (distab(n) + fuz))
dparcd = abs (accd(n) / (distcd(n) + fuz))
if (distab(n) .ne. 0.0) then
dpar = dparab
else
dpar = dparcd
endif
if (distcd(n) .eq. 0.0) then
dpar = distac(n)
endif
if (ipar(nn) .ne. 1) then
dpmin(nn) = 0.0
else
dpmin(nn) = dpar
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 (/ 'aptlnlc results:' /
cbug & (i3,' dpmin= ',1pe22.14,' ipar=',i2 /
cbug & ' fracab=',1pe22.14,' fraccd=',1pe22.14 /
cbug & ' eu,ev= ',1p2e22.14))
cbug write ( 3, 9903) (n, dpmin(n), ipar(n), fracab(n), fraccd(n),
cbug & eu(n), ev(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptlnlc. (+1 line.)
end
UCRL-WEB-209832