subroutine aptcirt (r1, r2, r3, tol, r4, r5, nerr)
c.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCIRT
c
c call aptcirt (r1, r2, r3, tol, r4, r5, nerr)
c
c Version: aptcirt Updated 2001 May 29 15:00.
c aptcirt Originated 1998 December 4 16:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the three coplanar and mutually tangent circles with radii
c r1, r2 and r3, find the radius r4 of the circle tangent to those
c three, in the space between them, and the radius r5 of an
c additional circle tangent to the first three, and not in the
c space between them.
c Flag nerr indicates any input error, if not zero.
c
c Method: Defining b = 1 / r, negative if the tangency is on the inside
c of the circle:
c s = b1 + b2 + b3 > 0
c q = b1 * b2 + b2 * b3 + b3 * b1 => 0
c sq = b1**2 + b2**2 + b3**2 > 0
c g = sq - 2.0 * q = b4 * b5
c
c If s > 0
c b4 = s + sqrt (q)
c b5 = g / b4
c If s < 0 (probably impossible):
c b5 = s - sqrt (q)
c b4 = g / b5
c
c r4 = 1 / b4 (1.e99 if b4 = 0)
c r5 = 1 / b5 (1.e99 if b5 = 0).
c
c Input: r1, r2, r3, tol.
c
c Output: r4, r5, nerr.
c
c Calls: none
c
c Glossary:
c
c nerr Output Indicates an input error, if not 0.
c 1 if b1*b2 + b2*b3 + b3*b4 is negative, which means
c radii r1, r2 and r3 are physically impossible.
c 2 if any of the radii are zero.
c 3 if r1 + r2 = 0.
c
c r1,r2,r3 Input The radii of three mutually tangent circles, all in
c the same plane. A negative radius indicates a circle
c that surrounds the other two circles. A very large
c radius, relative to the other two, indicates a
c straight line.
c
c r4 Output The radius of the circle coplanar with and tangent to
c the first three circles, in the space between them.
c E.g., r1 = 1, r2 = 2, r3 = 3, r4 = 0.2608695652174.
c Returned as -1.e99 if nerr = 1.
c
c r5 Output The radius of the circle coplanar with and tangent to
c the first three circles, but not in the space between
c them. A negative value of r5 means that circle
c surrounds the first four circles. E.g., r1 = 1,
c r2 = 2, r3 = 3, = 0.2608695652174..., r5 = -6.
c A value of r5 near 1.e99 indicates the circumference
c of that circle is a straight line. E.g., r1 = 1,
c r2 = 1, r3 = 0.25, r4 = 0.08333333333333, r5 = 1.e99.
c Returned as -1.e99 if nerr = 1.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit floating point numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcirt finding radii of circles tangent to',
cbug & ' circles with radii:' /
cbug & ' r1, r2, r3 =',1p3e22.14 /
cbug & ' tol =',1pe22.14 )
cbug write ( 3, 9901) r1, r2, r3, tol
cbug g = -1.e99
cbug q = -1.e99
cbug s = -1.e99
cbug b4 = -1.e99
cbug b5 = -1.e99
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
r4 = -1.e99
r5 = -1.e99
c.... Test for input errors.
if (((r1 .lt. 0.0) .and. (-r1 .lt. (r2 + r3))) .or.
& ((r2 .lt. 0.0) .and. (-r2 .lt. (r3 + r1))) .or.
& ((r3 .lt. 0.0) .and. (-r3 .lt. (r1 + r2)))) then
nerr = 1 ! Impossible fit.
go to 210
endif
if (r1 * r2 * r3 .eq. 0.0) then
nerr = 2 ! A radius is zero.
go to 210
endif
s12 = r1 + r2
if (abs (s12) .le. tol * (abs (r1) + abs (r2))) s12 = 0.0
if (s12 .eq. 0.0) then
nerr = 3
go to 210
endif
c.... Find the curvatures of circles 1, 2 and 3.
b1 = 1.0 / r1
b2 = 1.0 / r2
b3 = 1.0 / r3
cbugcbugc***DEBUG begins.
cbugcbug 9801 format (' b1, b2, b3 =',1p3e22.14)
cbugcbug write (3, 9801) b1, b2, b3
cbugcbugc***DEBUG ends.
c.... See if a solution is possible.
q = b1 * b2 + b2 * b3 + b3 * b1
cbugcbugc***DEBUG begins.
cbugcbug 9701 format (' q =',1pe22.14 )
cbugcbug write ( 3, 9701) q
cbugcbugc***DEBUG ends.
errq = 2.0 * tol * (abs (b1*b2) + abs (b2*b3) + abs (b3*b1))
if (abs (q) .le. errq) q = 0.0
if (q .lt. 0.0) then
nerr = 1
go to 210
endif
c.... See if there are one or two solutions.
s = b1 + b2 + b3
errs = tol * (abs (b1) + abs (b2) + abs (b3))
if (abs (s) .le. errs) s = 0.0
sq = b1**2 + b2**2 + b3**2
cbugcbugc***DEBUG begins.
cbugcbug 9702 format (' s =',1pe22.14 )
cbugcbug write ( 3, 9702) s
cbugcbugc***DEBUG ends.
errs = tol * (abs (b1) + abs (b2) + abs (b3))
if (abs (s) .le. errs) s = 0.0
if (q .eq. 0.0) then
b4 = s
b5 = s
go to 180
endif
c.... Find the differences between radii.
d12 = r2 - r1
err12 = tol * (abs (r1) + abs (r2))
if (abs (d12) .le. err12) d12 = 0.0
d23 = r3 - r2
err23 = tol * (abs (r2) + abs (r3))
if (abs (d23) .le. err23) d23 = 0.0
d31 = r1 - r3
err31 = tol * (abs (r1) + abs (r3))
if (abs (d31) .le. err31) d31 = 0.0
c.... Find g = b4 * b5, by the most accurate equations.
if ((d12 .eq. 0.0) .and. (d23 .eq. 0.0)) then
g = -3.0 * b1**2
errg = 0.0
elseif (d12 .eq. 0.0) then
g = b3 * (b3 - 4.0 * b1)
errg = 2.0 * tol * abs (b3) * (abs (b3) + 4.0 * abs (b1))
elseif (d23 .eq. 0.0) then
g = b1 * (b1 - 4.0 * b2)
errg = 2.0 * tol * abs (b1) * (abs (b1) + 4.0 * abs (b2))
elseif (d31 .eq. 0.0) then
g = b2 * (b2 - 4.0 * b3)
errg = 2.0 * tol * abs (b2) * (abs (b2) + 4.0 * abs (b3))
else ! Radii r1, r2 and r3 different.
g = sq - 2.0 * q
errg = 2.0 * tol * sq + 2.0 * errq
endif ! Tested for equal radii.
cbugcbugc***DEBUG begins.
cbugcbug 9802 format (' g =',1pe22.14)
cbugcbug write (3, 9802) g
cbugcbugc***DEBUG ends.
if (abs (g) .le. errg) g = 0.0
c.... Find the curvatures b4 and b5.
if (s .lt. 0.0) then ! This may be impossible.
b5 = s - 2.0 * sqrt (q)
b4 = g / b5
elseif (s .eq. 0.0) then
b4 = 2.0 * sqrt (q)
b5 = -2.0 * sqrt (q)
elseif (s .gt. 0.0) then
b4 = s + 2.0 * sqrt (q)
b5 = g / b4
endif ! Tested sign of s.
c.... Find the radii r4 and r5.
180 if (b4 .eq. 0.0) then
r4 = 1.e99
elseif (abs (b4) .ge. 1.e99) then
r4 = 0.0
else
r4 = 1.0 / b4
endif
if (b5 .eq. 0.0) then
r5 = 1.e99
elseif (abs (b5) .ge. 1.e99) then
r5 = 0.0
else
r5 = 1.0 / b5
endif
c.... Put radii r4 and r5 in increasing order of absolute value.
if (abs (r5) .lt. abs (r4)) then
rx = r4
r4 = r5
r5 = rx
endif
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptcirt results: nerr=',i3 /
cbug & ' q, s, g =',1p3e22.14 /
cbug & ' b4, b5 =',1p2e22.14 /
cbug & ' r4, r5 =',1p2e22.14 )
cbug write ( 3, 9902) nerr, q, s, g, b4, b5, r4, r5
cbugc***DEBUG ends.
return
c.... End of subroutine aptcirt. (+1 line.)
end
UCRL-WEB-209832