subroutine aptspht (r1, r2, r3, r4, tol, r5, r6, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPHT
c
c call aptspht (r1, r2, r3, r4, tol, r5, r6, nerr)
c
c Version: aptspht Updated 2001 May 29 15:20.
c aptspht Originated 1999 January 5 16:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the four mutually tangent spheres with radii r1, r2, r3
c and r4, find the radii r5 and r6 of two additional spheres,
c each tangent to the first four spheres.
c Flag nerr indicates any input error, if not zero.
c
c Method: Defining the curvature b = 1 / r, negative if the tangency
c is on the inside of the sphere:
c s = b1 + b2 + b3 + b4 => 0
c p = b1*b2 + b2*b3 + b3*b4 + b4*b1 + b1*b3 + b2*b4
c t = b1**2 + b2**2 + b3**2 + b4**2 > 0
c q = 6 * p - 3 * t => 0
c g = t - p = b5 * b6
c If s < 0 (impossible):
c b6 = 0.5 * s - 0.5 * sqrt (q)
c b5 = g / b6
c If s = 0 (impossible?):
c b5 = 0.5 * sqrt (q)
c b6 = -0.5 * sqrt (q)
c If s > 0:
c b5 = 0.5 * s + 0.5 * sqrt (q)
c b6 = g / b5
c (Other more accurate forms are used when any of the
c first four radii are equal.)
c r5 = 1 / b5 (1.e99 if b5 = 0)
c r6 = 1 / b6 (1.e99 if b6 = 0).
c Special case: if r1 + r2 + r3 = 0. then r5 = r6 = r4.
c
c Input: r1, r2, r3, r4, tol.
c
c Output: r5, r6, nerr.
c
c Calls: none
c
c Glossary:
c
c nerr Output Indicates an input error, if not 0.
c 1 if radii r1, r2, r3 and r4 are physically
c impossible.
c 2 if any of the radii are zero.
c 3 if r1 + r2 = 0.
c
c r1 - r4 Input The radii of four mutually tangent spheres.
c A negative radius indicates a sphere that surrounds
c the other three spheres. A very large radius,
c relative to the other three, indicates a plane.
c
c r5 Output The radius of the sphere tangent to the first four
c spheres, in the space between them.
c Returned as -1.e99 if nerr = 1 or 2.
c
c r6 Output The radius of the sphere tangent to the first four
c spheres.
c A negative value of r5 means that sphere surrounds
c the first four spheres.
c A value of r5 near 1.e99 indicates a plane.
c Returned as -1.e99 if nerr = 1 or 2.
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 (/ 'aptspht finding radii of spheres tangent to',
cbug & ' spheres with radii:' /
cbug & ' r1, r2 =',1p2e22.14 /
cbug & ' r3, r4 =',1p2e22.14 /
cbug & ' tol = ',1pe22.14 )
cbug write ( 3, 9901) r1, r2, r3, r4, tol
cbug b5 = -1.e99
cbug b6 = -1.e99
cbug q = -1.e99
cbug s = -1.e99
cbug t = -1.e99
cbug g = -1.e99
cbugc***DEBUG ends.
c.... Initialize.
r5 = -1.e99
r6 = -1.e99
nerr = 0
c.... Test for input errors.
if (r1 * r2 * r3 * r4 .eq. 0.0) then
nerr = 2
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 the first four spheres.
b1 = 1.0 / r1
b2 = 1.0 / r2
b3 = 1.0 / r3
b4 = 1.0 / r4
cbugcbugc***DEBUG begins.
cbugcbug 9801 format (' b1, b2 = ',1p3e22.14 /
cbugcbug & ' b3, b4 = ',1p2e22.14 )
cbugcbug write (3, 9801) b1, b2, b3, b4
cbugcbugc***DEBUG ends.
c.... Find the curvatures of the spheres tangent to spheres 1, 2, 3 and 4.
s = b1 + b2 + b3 + b4
errs = tol * (abs (b1) + abs (b2) + abs (b3) + abs (b4))
cbugcbugc***DEBUG begins.
cbugcbug 9802 format (' s, errs= ',1p2e22.14)
cbugcbug write (3, 9802) s, errs
cbugcbugc***DEBUG ends.
if (abs (s) .le. errs) s = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write (3, 9802) s, errs
cbugcbugc***DEBUG ends.
c.... Test for the special case r1 + r2 + r3 = 0, r1 + r2 /= 0.
scir = r1 + r2 + r3
errsc = tol * (abs (r1) + abs (r2) + abs (r3))
if (abs (scir) .le. errsc) scir = 0.0
if (scir .eq. 0.0) then ! Special case.
db = b1 + b2 + b3 - b4
if (abs (db) .le. errs) then ! Qualifies.
r5 = r4
r6 = r4
go to 210
else
nerr = 1
go to 210
endif
endif
c.... Solve the more general case.
p = b1*b2 + b2*b3 + b3*b1 + b1*b4 + b2*b4 + b3*b4
errp = tol * (abs (b1 * b2) + abs (b2 * b3) +
& abs (b3 * b1) + abs (b1 * b4) +
& abs (b2 * b4) + abs (b3 * b4) )
cbugcbugc***DEBUG begins.
cbugcbug 9803 format (' p, errp= ',1p2e22.14)
cbugcbug write (3, 9803) p, errp
cbugcbugc***DEBUG ends.
if (abs (p) .le. errp) p = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write (3, 9803) p, errp
cbugcbugc***DEBUG ends.
t = b1**2 + b2**2 + b3**2 + b4**2
errt = tol * (b1**2 + b2**2 + b3**2 + b4**2)
cbugcbugc***DEBUG begins.
cbugcbug 9804 format (' t, errt= ',1p2e22.14)
cbugcbug write (3, 9804) t, errt
cbugcbugc***DEBUG ends.
q = 6.0 * p - 3.0 * t
errq = 6.0 * errp + 3.0 * errt
cbugcbugc***DEBUG begins.
cbugcbug 9805 format (' q, errq= ',1p2e22.14)
cbugcbug write (3, 9805) q, errq
cbugcbugc***DEBUG ends.
if (abs (q) .le. errq) q = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write (3, 9805) q, errq
cbugcbugc***DEBUG ends.
if (q .lt. 0.0) then
nerr = 1
go to 210
endif
c.... Find the differences between radii.
d12 = r2 - r1
err12 = tol * (abs (r1) + abs (r2))
if (abs (d12) .le. err12) d12 = 0.0
d13 = r3 - r1
err13 = tol * (abs (r1) + abs (r3))
if (abs (d13) .le. err13) d13 = 0.0
d14 = r4 - r1
err14 = tol * (abs (r1) + abs (r4))
if (abs (d14) .le. err14) d14 = 0.0
d23 = r3 - r2
err23 = tol * (abs (r2) + abs (r3))
if (abs (d23) .le. err23) d23 = 0.0
d24 = r4 - r2
err24 = tol * (abs (r2) + abs (r4))
if (abs (d24) .le. err24) d24 = 0.0
d34 = r4 - r3
err34 = tol * (abs (r3) + abs (r4))
if (abs (d34) .le. err34) d34 = 0.0
c.... Find q and g = b5 * b6 = t - p, by the most accurate equation..
if ((d12 .eq. 0.0) .and. (d23 .eq. 0.0) .and. (d34 .eq. 0.0)) then
q = 24.0 * b1**2
errq = 0.0
g = -2.0 * b1**2
errg = 0.0
elseif ((d12 .eq. 0.0) .and. (d23 .eq. 0.0)) then
q = 9.0 * b1 * (b1 + 2.0 * b4) - 3.0 * b4**2
errq = tol * (9.0 * b1**2 + 18.0 * abs (b1 * b4) + 3.0 * b4**2)
g = b4 * (b4 - 3.0 * b1)
errg = tol * (b4**2 - 3.0 * abs (b1 * b4))
elseif ((d23 .eq. 0.0) .and. (d34 .eq. 0.0)) then
q = 9.0 * b2 * (b2 + 2.0 * b1) - 3.0 * b1**2
errq = tol * (9.0 * b2**2 + 18.0 * abs (b2 * b1) + 3.0 * b1**2)
g = b1 * (b1 - 3.0 * b2)
errg = tol * (b1**2 - 3.0 * abs (b1 * b2))
elseif ((d14 .eq. 0.0) .and. (d34 .eq. 0.0)) then
q = 9.0 * b3 * (b3 + 2.0 * b2) - 3.0 * b2**2
errq = tol * (9.0 * b3**2 + 18.0 * abs (b3 * b2) + 3.0 * b2**2)
g = b2 * (b2 - 3.0 * b3)
errg = tol * (b2**2 - 3.0 * abs (b2 * b3))
elseif ((d14 .eq. 0.0) .and. (d24 .eq. 0.0)) then
q = 9.0 * b4 * (b4 + 2.0 * b3) - 3.0 * b3**2
errq = tol * (9.0 * b4**2 + 18.0 * abs (b4 * b3) + 3.0 * b3**2)
g = b3 * (b3 - 3.0 * b4)
errg = tol * (b3**2 - 3.0 * abs (b3 * b4))
elseif ((d12 .eq. 0.0) .and. (d34 .eq. 0.0)) then
q = 24.0 * b1 * b3
errq = tol * 24.0 * abs (b1 * b3)
g = b1**2 + b3**2 - 4.0 * b1 * b3
errg = tol * (b1**2 + b3**2 + 4.0 * abs (b1 * b3))
elseif ((d13 .eq. 0.0) .and. (d24 .eq. 0.0)) then
q = 24.0 * b1 * b2
errq = tol * 24.0 * abs (b1 * b2)
g = b1**2 + b2**2 - 4.0 * b1 * b2
errg = tol * (b1**2 + b2**2 + 4.0 * abs (b1 * b2))
elseif ((d14 .eq. 0.0) .and. (d23 .eq. 0.0)) then
q = 24.0 * b1 * b2
errq = tol * 24.0 * abs (b1 * b2)
g = b1**2 + b2**2 - 4.0 * b1 * b2
errg = tol * (b1**2 + b2**2 + 4.0 * abs (b1 * b2))
elseif (d12 .eq. 0.0) then
q = 12.0 * b1 * (b3 + b4) - 3.0 * (b3 - b4)**2
errq = tol * (12.0 * (abs (b1 * b3) + abs (b1 * b4)) +
& 6.0 * abs (b3 * b4) + 3.0 * (b3**2 + b4**2))
g = b1**2 + b3**2 + b4**2 - 2.0 * b1 * (b3 + b4) - b3 * b4
errg = tol * (b1**2 + b3**2 + b4**2 + 2.0 * abs (b1 * b3) +
& 2.0 * abs (b1 * b4) + 2.0 * abs (b3 * b4))
elseif (d23 .eq. 0.0) then
q = 12.0 * b2 * (b4 + b1) - 3.0 * (b4 - b1)**2
errq = tol * (12.0 * (abs (b2 * b4) + abs (b2 * b1)) +
& 6.0 * abs (b4 * b1) + 3.0 * (b4**2 + b1**2))
g = b2**2 + b4**2 + b1**2 - 2.0 * b2 * (b4 + b1) - b4 * b1
errg = tol * (b2**2 + b4**2 + b1**2 + 2.0 * abs (b2 * b4) +
& 2.0 * abs (b2 * b1) + 2.0 * abs (b4 * b1))
elseif (d34 .eq. 0.0) then
q = 12.0 * b3 * (b1 + b2) - 3.0 * (b1 - b2)**2
errq = tol * (12.0 * (abs (b3 * b1) + abs (b3 * b2)) +
& 6.0 * abs (b1 * b2) + 3.0 * (b1**2 + b2**2))
g = b3**2 + b1**2 + b2**2 - 2.0 * b3 * (b1 + b2) - b1 * b2
errg = tol * (b3**2 + b1**2 + b2**2 + 2.0 * abs (b3 * b1) +
& 2.0 * abs (b3 * b2) + 2.0 * abs (b1 * b2))
elseif (d14 .eq. 0.0) then
q = 12.0 * b4 * (b2 + b3) - 3.0 * (b2 - b3)**2
errq = tol * (12.0 * (abs (b4 * b2) + abs (b4 * b3)) +
& 6.0 * abs (b2 * b3) + 3.0 * (b2**2 + b3**2))
g = b4**2 + b2**2 + b3**2 - 2.0 * b4 * (b2 + b3) - b2 * b3
errg = tol * (b4**2 + b2**2 + b3**2 + 2.0 * abs (b4 * b2) +
& 2.0 * abs (b4 * b3) + 2.0 * abs (b2 * b3))
elseif (d13 .eq. 0.0) then
q = 12.0 * b1 * (b2 + b4) - 3.0 * (b2 - b4)**2
errq = tol * (12.0 * (abs (b1 * b2) + abs (b1 * b4)) +
& 6.0 * abs (b2 * b4) + 3.0 * (b2**2 + b4**2))
g = b1**2 + b2**2 + b4**2 - 2.0 * b1 * (b2 + b4) - b2 * b4
errg = tol * (b1**2 + b2**2 + b4**2 + 2.0 * abs (b1 * b2) +
& 2.0 * abs (b1 * b4) + 2.0 * abs (b2 * b4))
elseif (d24 .eq. 0.0) then
q = 12.0 * b2 * (b3 + b1) - 3.0 * (b3 - b1)**2
errq = tol * (12.0 * (abs (b2 * b3) + abs (b2 * b1)) +
& 6.0 * abs (b3 * b1) + 3.0 * (b3**2 + b1**2))
g = b2**2 + b3**2 + b1**2 - 2.0 * b2 * (b3 + b1) - b3 * b1
errg = tol * (b2**2 + b3**2 + b1**2 + 2.0 * abs (b2 * b3) +
& 2.0 * abs (b2 * b1) + 2.0 * abs (b3 * b1))
else
q = 6.0 * p - 3.0 * t
errq = 6.0 * errp + 3.0 * errt
g = t - p
errg = errt + errp
endif
cbugcbugc***DEBUG begins.
cbugcbug 9807 format (' g, errg= ',1p2e22.14)
cbugcbug write (3, 9807) g, errg
cbugcbugc***DEBUG ends.
if (abs (q) .le. errq) q = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write (3, 9805) q, errq
cbugcbugc***DEBUG ends.
if (q .lt. 0.0) then
nerr = 1
go to 210
endif
if (abs (g) .le. errg) g = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write (3, 9807) g, errg
cbugcbugc***DEBUG ends.
c.... Find the curvatures of spheres 5 and 6.
if (s .lt. 0.0) then
b6 = 0.5 * s - 0.5 * sqrt (q)
b5 = g / b6
elseif (s .eq. 0.0) then
b5 = 0.5 * sqrt (q)
b6 = -0.5 * sqrt (q)
else
b5 = 0.5 * s + 0.5 * sqrt (q)
b6 = g / b5
endif
c.... Find the radii of spheres 5 and 6.
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
if (b6 .eq. 0.0) then
r6 = 1.e99
elseif (abs (b6) .ge. 1.e99) then
r6 = 0.0
else
r6 = 1.0 / b6
endif
c.... Return the smaller absolute radius as r5.
if (abs (r6) .lt. abs (r5)) then
rx = r5
r5 = r6
r6 = rx
endif
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptspht results: nerr=',i3 /
cbug & ' q, s, p =',1p3e22.14 /
cbug & ' t, g =',1p2e22.14 /
cbug & ' b5, b6 =',1p2e22.14 /
cbug & ' r5, r6 =',1p2e22.14 )
cbug write ( 3, 9902) nerr, q, s, p, t, g, b5, b6, r5, r6
cbugc***DEBUG ends.
return
c.... End of subroutine aptspht. (+1 line.)
end
UCRL-WEB-209832