subroutine aptcirc (au, av, bu, bv, cu, cv, np, tol,
& pu, pv, rad, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCIRC
c
c call aptcirc (au, av, bu, bv, cu, cv, np, tol,
c & pu, pv, rad, nerr)
c
c Version: aptcirc Updated 1993 April 28 10:00.
c aptcirc Originated 1993 April 28 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For each of the np sets of 3 points a = (au, av), b = (bu, bv)
c and c = (cu, cv), all in the uv plane, to find the center
c p = (pu, pv) and radius rad of the circle through points
c (a, b, c).
c Flag nerr indicates any input error, if not zero.
c
c Method: Point "p" is on the lines bisecting vectors "ab" and "ac":
c (p - 0.5 * (b + a)) dot (b - a) = 0
c (p - 0.5 * (c + a)) dot (c - a) = 0
c These 2 equations are solved for pu and pv.
c
c Input: au, av, bu, bv, cu, cv, np, tol.
c
c Output: pu, pv, rad, nerr.
c
c Calls: aptvdic, aptvsuc, aptvdoc
c
c
c Glossary:
c
c au, av Input The u and v coordinates of point "a". Size np.
c
c bu, bv Input The u and v coordinates of point "b". Size np.
c
c cu, cv Input The u and v coordinates of point "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 The number of sets of points (a, b, c) for which the
c circle is to be found. Must be positive.
c
c pu, pv Output The u and v coordinates of the center of the circle.
c Each 1.e99 if the points (a, b, c) are colinear or
c congruent. Size np.
c
c rad Output The radius of the circle, larger than 1.e99 if the
c points (a, b, c) are colinear or congruent. 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 input point "a".
dimension au (1)
c---- Coordinate v of input point "a".
dimension av (1)
c---- Coordinate u of input point "b".
dimension bu (1)
c---- Coordinate v of input point "b".
dimension bv (1)
c---- Coordinate u of input point "c".
dimension cu (1)
c---- Coordinate v of input point "c".
dimension cv (1)
c---- Coordinate u of center of circle "p".
dimension pu (1)
c---- of center of circle "p".
dimension pv (1)
c---- Radius of circle.
dimension rad (1)
c.... Local variables. Floating point.
c---- Dot product of "ab" and "abmid".
common /laptcirc/ abdot (64)
c---- Length of vector "abmid".
common /laptcirc/ abmidl (64)
c---- Component x of vector "abmid".
common /laptcirc/ abmidu (64)
c---- Component y of vector "abmid".
common /laptcirc/ abmidv (64)
c---- Length of vector "ab".
common /laptcirc/ abvlen (64)
c---- Component x of vector "ab".
common /laptcirc/ abu (64)
c---- Component y of vector "ab".
common /laptcirc/ abv (64)
c---- Dot product of "ac" and "acmid".
common /laptcirc/ acdot (64)
c---- Length of vector "acmid".
common /laptcirc/ acmidl (64)
c---- Component x of vector "acmid".
common /laptcirc/ acmidu (64)
c---- Component y of vector "acmid".
common /laptcirc/ acmidv (64)
c---- Length of vector "ac".
common /laptcirc/ acvlen (64)
c---- Component x of vector "ac".
common /laptcirc/ acu (64)
c---- Component y of vector "ac".
common /laptcirc/ acv (64)
c---- A big number.
common /laptcirc/ big
c.... Local variables. Integers.
c---- Index in arrays.
common /laptcirc/ n
c---- First index of subset of data.
common /laptcirc/ n1
c---- Last index of subset of data.
common /laptcirc/ n2
c---- Index in external array.
common /laptcirc/ nn
c---- Size of current subset of data.
common /laptcirc/ ns
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcirc finding circle through points:' /
cbug & (i3,' au,av=',1p2e22.14 /
cbug & ' bu,bv=',1p2e22.14 /
cbug & ' cu,cv=',1p2e22.14))
cbug write ( 3, 9901) (n, au(n), av(n), bu(n), bv(n),
cbug & cu(n), cv(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
big = 1.e99
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 difference vectors "ab" and "ac".
call aptvdic (au(n1), av(n1), bu(n1), bv(n1), ns,
& tol, abu, abv, abvlen, nerr)
call aptvdic (au(n1), av(n1), cu(n1), cv(n1), ns,
& tol, acu, acv, acvlen, nerr)
c.... Find the midpoint vectors "abmid" and "acmid".
call aptvsuc (0, 0.5, au(n1), av(n1),
& 0.5, bu(n1), bv(n1), ns, tol,
& abmidu, abmidv, abmidl, nerr)
call aptvsuc (0, 0.5, au(n1), av(n1),
& 0.5, cu(n1), cv(n1), ns, tol,
& acmidu, acmidv, acmidl, nerr)
c.... Find the vector dot products ab.abmid and ac.acmid.
call aptvdoc (abu, abv, abmidu, abmidv, ns, tol, abdot, nerr)
call aptvdoc (acu, acv, acmidu, acmidv, ns, tol, acdot, nerr)
c.... Find the center of each circle.
do 120 n = 1, ns
nn = n1 + n - 1
det = abu(n) * acv(n) - abv(n) * acu(n)
deterr = tol * (abs (abu(n) * acv(n)) + abs (abv(n) * acu(n)))
if (abs (det) .le. deterr) then
pu(nn) = big
pv(nn) = big
else
pu(nn) = (acv(n) * abdot(n) - abv(n) * acdot(n)) / det
pv(nn) = (abu(n) * acdot(n) - acu(n) * abdot(n)) / det
endif
120 continue
c.... Find the radius of each circle.
call aptvdic (au(n1), av(n1), pu(n1), pv(n1), ns,
& tol, abu, abv, rad(n1), 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
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptcirc results:' /
cbug & (i3,' pu,pv=',1p2e22.14 /
cbug & ' rad= ',1pe22.14))
cbug write ( 3, 9902) (n, pu(n), pv(n), rad(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptcirc. (+1 line.)
end
UCRL-WEB-209832