subroutine aptcinc (ra, au, av, rb, bu, bv, np, tol,
& cu, cv, du, dv, nint, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCINC
c
c call aptcinc (ra, au, av, rb, bu, bv, np, tol,
c & cu, cv, du, dv, nint, nerr)
c
c Version: aptcinc Updated 1990 November 27 14:00.
c aptcinc Originated 1990 March 20 14: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 points of
c intersection c = (cu, cv) and d = (du, dv) of the circle
c of radius ra at point a = (au, av) and the circle of radius
c rb at point b = (bu, bv), if an intersection occurs.
c Flag nint indicates the number of intersection points.
c Flag nerr indicates any input error.
c
c Input: ra, au, av, rb, bu, bv, np, tol.
c
c Output: cu, cv, du, dv, nint, nerr.
c
c Calls: aptvdic, aptvadc, aptvplc, aptvuac
c
c
c Glossary:
c
c au, av Input The u and v coordinates of point "a" at the center
c of the circle with radius ra, in the uv plane (2-D).
c Size np.
c
c bu, bv Input The u and v coordinates of point "b" at the center
c of the circle with radius rb, in the uv plane (2-D).
c Size np.
c
c cu, cv Output The u and v coordinates of point "c" at an intersection
c of the two circles centered at points "a" and "b", if
c an intersection occurs (nint = 1 or 2). Size np.
c Same as "d" if nint = 1.
c Meaningless, but set to "a" if nint = 0 or 3.
c
c du, dv Output The u and v coordinates of point "d" at an intersection
c of the two circles centered at points "a" and "b", if
c an intersection occurs (nint = 1 or 2). Size np.
c Same as "c" if nint = 1.
c Meaningless, but set to "a" if nint = 0 or 3.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c nint Output Indicates the number of intersection points:
c 0 if no intersection occurs. Ignore points "c", "d".
c 1 if the circles are tangent at the single point
c "c" = "d".
c 2 if the circles overlap, intersecting at the two
c points "c" and "d".
c 3 if the circles are congruent. The intersection
c includes each circle. Ignore points "c", "d".
c Size np.
c
c np Input Size of arrays.
c
c ra Input The radius of the circle centered at point "a".
c Size np. The absolute value is used.
c
c rb Input The radius of the circle centered at point "b".
c Size np. The absolute value is used.
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---- Coordinate u of point "d".
dimension du (1)
c---- Coordinate v of point "d".
dimension dv (1)
c---- Number of intersection points.
dimension nint (1)
c---- Radius of circle at point "a".
dimension ra (1)
c---- Radius of circle at point "b".
dimension rb (1)
c.... Local variables.
c---- Component u of unit vector along "ab".
common /laptcinc/ abu (64)
c---- Component v of unit vector along "ab".
common /laptcinc/ abv (64)
c---- Factor fact1*fact2*fact3*fact4.
common /laptcinc/ arg
c---- Distance from point "c" to point "d".
common /laptcinc/ cd (64)
c---- Magnitude of vector "ab".
common /laptcinc/ dab (64)
c---- Fractional distance of "e" along "ab".
common /laptcinc/ dae (64)
c---- Estimated error in factor of dae.
common /laptcinc/ daerr
c---- Coordinate u of midpoint of line "cd".
common /laptcinc/ eu (64)
c---- Coordinate v of midpoint of line "cd".
common /laptcinc/ ev (64)
c---- Factor dab + rpr.
common /laptcinc/ fact1
c---- Factor dab - rmr.
common /laptcinc/ fact2
c---- Factor rpr - dab
common /laptcinc/ fact3
c---- Factor rmr + dab.
common /laptcinc/ fact4
c---- Estimated error in fact2, 3, 4.
common /laptcinc/ ferr
c---- A very small number.
common /laptcinc/ fuz
c---- Component u of vector normal to "ab".
common /laptcinc/ hnu (64)
c---- Component v of vector normal to "ab".
common /laptcinc/ hnv (64)
c---- Index in arrays.
common /laptcinc/ n
c---- First index of subset of data.
common /laptcinc/ n1
c---- Last index of subset of data.
common /laptcinc/ n2
c---- Index in external array.
common /laptcinc/ nn
c---- Size of current subset of data.
common /laptcinc/ ns
c---- Difference abs(rb) - abs(ra).
common /laptcinc/ rmr
c---- Sum abs(rb) + abs(ra).
common /laptcinc/ rpr
c---- Length of a vector.
common /laptcinc/ vlen (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcinc finding the intersections of circles.' /
cbug & ' tol=',1pe13.5)
cbug 9902 format (' ra,au,av=',1p3e22.14 / ' rb,bu,bv=',1p3e22.14)
cbug write (3, 9901) tol
cbug write (3, 9902) (ra(n), au(n), av(n), rb(n), bu(n), bv(n),
cbug & 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
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the unit vector "ab" from point "a" to point "b".
c.... The magnitude dab may be truncated to zero.
call aptvdic (au(n1), av(n1), bu(n1), bv(n1), ns,
& tol, abu, abv, dab, nerr)
call aptvuac (abu, abv, ns, tol, dab, nerr)
c.... Find the unit vector in the uv plane perpendicular to the line "ab".
call aptvplc (au(n1), av(n1), bu(n1), bv(n1), ns, tol,
& hnu, hnv, vlen, nerr)
call aptvuac (hnu, hnv, ns, tol, vlen, nerr)
c.... Find the distance between points "c" and "d", and the distance of the
c.... midpoint of the line "cd" from point "a".
c---- No truncation tests.
if (tol .le. 0.0) then
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
rpr = abs (rb(nn)) + abs (ra(nn))
rmr = abs (rb(nn)) - abs (ra(nn))
fact1 = dab(n) + rpr
fact2 = dab(n) - rmr
fact3 = rpr - dab(n)
fact4 = rmr + dab(n)
arg = fact1 * fact2 * fact3 * fact4
cd(n) = sqrt (amax1 (0.0, arg)) / (dab(n) + fuz)
dae(n) = 0.5 * (dab(n)**2 - rpr * rmr) / (dab(n) + fuz)
if (arg .gt. 0.0) then
nint(nn) = 2
else
nint(nn) = 0
endif
if (arg .eq. 0.0) then
nint(nn) = 1
endif
if ((dab(n) .eq. 0.0) .and.
& (rmr .le. 0.0) ) then
nint(nn) = 3
endif
if (nint(nn) .eq. 0) then
dae(n) = 0.0
endif
c---- End of loop over subset of data.
120 continue
c---- Test for truncation.
else
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
rpr = abs (rb(nn)) + abs (ra(nn))
rmr = abs (rb(nn)) - abs (ra(nn))
if (abs (rmr) .lt. tol * rpr) then
rmr = 0.0
endif
fact1 = dab(n) + rpr
ferr = tol * fact1
fact2 = dab(n) - rmr
if (abs (fact2) .lt. ferr) then
fact2 = 0.0
endif
fact3 = rpr - dab(n)
if (abs (fact3) .lt. ferr) then
fact3 = 0.0
endif
fact4 = rmr + dab(n)
if (abs (fact4) .lt. ferr) then
fact4 = 0.0
endif
arg = fact1 * fact2 * fact3 * fact4
cd(n) = sqrt (amax1 (0.0, arg)) / (dab(n) + fuz)
if (arg .gt. 0.0) then
nint(nn) = 2
else
nint(nn) = 0
endif
if (arg .eq. 0.0) then
nint(nn) = 1
endif
if ((dab(n) .le. tol) .and.
& (abs (rmr) .le. tol) ) then
nint(nn) = 3
endif
dae(n) = dab(n)**2 - rpr * rmr
daerr = tol * (dab(n)**2 + rpr * abs (rmr))
if (abs (dae(n)) .lt. daerr) then
dae(n) = 0.0
endif
dae(n) = 0.5 * dae(n) / (dab(n) + fuz)
if (nint(nn) .eq. 0) then
dae(n) = 0.0
endif
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
endif
c.... Find the coordinates of point "e", the center of the "line"
c.... of intersection. e = a + dae * ab.
call aptvadc (au(n1), av(n1), 1., dae, abu, abv, ns, tol,
& eu, ev, vlen, nerr)
c.... Find the coordinates of the two points of intersection, "c" and "d".
call aptvadc (eu, ev, 0.5, cd, hnu, hnv, ns, tol,
& cu(n1), cv(n1), vlen, nerr)
call aptvadc (eu, ev, -0.5, cd, hnu, hnv, ns, tol,
& du(n1), dv(n1), vlen, nerr)
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
210 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptcinc results. nerr =',i3)
cbug 9904 format (' n =',i5,' nint =',i2 /
cbug & ' cu,cv=',1p2e22.14 / ' du,dv=',1p2e22.14)
cbug write (3, 9903) nerr
cbug write (3, 9904) (n, nint(n), cu(n), cv(n), du(n), dv(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
return
c.... End of subroutine aptcinc. (+1 line.)
end
UCRL-WEB-209832