subroutine aptsphk (r1, r2, r3, r4, tol, r5, r6,
& x1, y1, z1, x2, y2, z2, x3, y3, z3,
& x4, y4, z4, x5, y5, z5, x6, y6, z6, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPHK
c
c call aptsphk (r1, r2, r3, r4, tol, r5, r6,
c & x1, y1, z1, x2, y2, z2, x3, y3, z3,
c & x4, y4, z4, x5, y5, z5, x6, y6, z6, nerr)
c
c Version: aptsphk Updated 2001 May 29 15:30.
c aptsphk Originated 1999 February 1 16:20.
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 and
c r4, find the radii r5 and r6 of two additional spheres, each
c tangent to those four. Also find the coordinates of the centers
c of each of the 6 spheres, with the first at the origin, the
c second on the x axis, the third in the x-y plane.
c Flag nerr indicates any input or result error, if not zero.
c
c Method: See aptcirk, aptcirl, aptcirt, aptsphl, aptspht.
c
c Input: r1, r2, r3, r4, tol.
c
c Output: r5, r6, x1, y1, z1, x2, y2, z2, x3, y3, z3,
c x4, y4, z4, x5, y5, z5, x6, y6, z6, nerr.
c
c Calls: aptcirl, aptcirt
c
c Glossary:
c
c nerr Output Indicates an input error, if not 0.
c -1 if one or more of the spheres is so small, and its
c y and/or z center coordinate so small, that its
c sign or magnitude can not be determined.
c 1 if two or more of the first four radii are
c negative, or one radius is negative, but not large
c enough in magnitude for that sphere to enclose the
c other three, or one sphere has a radius too small
c to touch each of the other three spheres.
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 a sphere tangent to the first four
c spheres, in the space between them, if the first
c four radii are all positive.
c Returned as -1.e99 if nerr is positive.
c
c r6 Output The radius of another sphere tangent to the first four
c spheres. A negative value of r6 means that sphere 6
c surrounds the first four spheres.
c A value of r6 near 1.e99 indicates the surface of
c that sphere is a plane.
c Returned as -1.e99 if nerr is positive.
c
c x, y, z Output The x, y, z coordinates of the center of a sphere,
c i.e., (xn, yn, zn) for the n'th sphere, n = 1, 6.
c Returned as -1.e99 if nerr is positive.
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 (/ 'aptsphk finding 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
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
r5 = -1.e99
r6 = -1.e99
x1 = -1.e99
y1 = -1.e99
z1 = -1.e99
x2 = -1.e99
y2 = -1.e99
z2 = -1.e99
x3 = -1.e99
y3 = -1.e99
z3 = -1.e99
x4 = -1.e99
y4 = -1.e99
z4 = -1.e99
x5 = -1.e99
y5 = -1.e99
z5 = -1.e99
x6 = -1.e99
y6 = -1.e99
z6 = -1.e99
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.... Test for special case r1 + r2 + r3 = 0, r4 just fits.
s123 = r1 + r2 + r3
errs = tol * (abs (r1) + abs (r2) + abs (r3))
if (abs (s123) .le. errs) then ! Special case.
b1 = 1.0 / r1
b2 = 1.0 / r2
b3 = 1.0 / r3
b4 = sqrt (b1**2 + b2**2 + b3**2)
rx = 1.0 / b4
if (abs (rx - r4) .gt. tol * (rx + r4)) then
nerr = 1
go to 210
endif
r5 = rx
r6 = rx
x1 = 0.0
y1 = 0.0
z1 = 0.0
x2 = r1 + r2
y2 = 0.0
z2 = 0.0
x3 = r2
y3 = 0.0
z3 = 0.0
x4 = r1 + r4 * (r1 - r2) / (r1 + r2)
y4 = 2.0 * r4
z4 = 0.0
x5 = x4
y5 = r4
z5 = r4 * sqrt (3.0)
x6 = x4
y6 = r4
z6 = -z5
go to 210
endif ! Special case.
c.... Find the center of sphere 1. Arbitrarily at origin.
x1 = 0.0
y1 = 0.0
z1 = 0.0
c.... Find the center of sphere 2. Arbitrarily on x axis.
x2 = r1 + r2
y2 = 0.0
z2 = 0.0
c.... Find the center of sphere 3. Arbitrarily at z = 0.
call aptcirl (r1, r2, r3, tol, x3, y3, nerra)
if (nerra .ne. 0) then
nerr = 1
go to 210
endif
z3 = 0.0
c.... Find the radii of the spheres tangent to spheres 1, 2, 3 and 4.
call aptspht (r1, r2, r3, r4, tol, r5, r6, nerra)
if (nerra .gt. 0) then
nerr = 1
go to 210
endif
c.... Find the center of sphere 4.
call aptsphl (r1, r2, r3, r4, tol, x4, y4, z4, nerra)
z4 = abs (z4)
if (nerra .lt. 0) then
nerr = -1
endif
if (nerra .gt. 0) then
nerr = 1
go to 210
endif
c.... Find the center of sphere 5.
call aptsphl (r1, r2, r3, r5, tol, x5, y5, z5, nerra)
if (nerra .gt. 0) then
nerr = 1
go to 210
endif
if (nerra .lt. 0) then
nerr = -1
endif
argt = r1 * (r1 + r4 + r5) - (r4 * r5 + x4 * x5 + y4 * y5)
test1 = argt - z4 * z5
test2 = argt + z4 * z5
terr = 2.0 * tol * (r1**2 + (abs (r1 * r4) + abs (r1 * r5)) +
& abs (r4 * r5) + abs (x4 * x5) +
& abs (y4 * y5) + abs (z4 * z5))
cbugcbugc***DEBUG begins.
cbugcbug 9842 format ('DEBUG z: test1,test2=',1p2e22.14)
cbugcbug write ( 3, 9842) test1, test2
cbugcbugc***DEBUG ends.
if (abs (test1) .le. terr) test1 = 0.0
if (abs (test2) .le. terr) test2 = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9842) test1, test2
cbugcbugc***DEBUG ends.
if (abs (test2) .lt. abs (test1)) then
z5 = -z5
endif
if (abs (test1) .eq. abs (test2)) then
if (z5 .ne. 0.0) then
nerr = -1
endif
endif
c.... Find the center of sphere 6.
call aptsphl (r1, r2, r3, r6, tol, x6, y6, z6, nerra)
if (nerra .gt. 0) then
nerr = 1
go to 210
endif
if (nerra .lt. 0) then
nerr = -1
endif
argt = r1 * (r1 + r4 + r6) - (r4 * r6 + x4 * x6 + y4 * y6)
test1 = argt - z4 * z6
test2 = argt + z4 * z6
terr = 2.0 * tol * (r1**2 + abs (r1 * r4) + abs (r1 * r6) +
& abs (r4 * r6) + abs (x4 * x6) +
& abs (y4 * y6) + abs (z4 * z6))
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9842) test1, test2
cbugcbugc***DEBUG ends.
if (abs (test1) .le. terr) test1 = 0.0
if (abs (test2) .le. terr) test2 = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9842) test1, test2
cbugcbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug 9851 format ('DEBUG 1: z5,z6=',1p2e22.14)
cbugcbug write ( 3, 9851) z5, z6
cbugcbugc***DEBUG ends.
if (abs (test2) .lt. abs (test1)) then
z6 = -z6
endif
cbugcbugc***DEBUG begins.
cbugcbug 9852 format ('DEBUG 2: z5,z6=',1p2e22.14)
cbugcbug write ( 3, 9852) z5, z6
cbugcbugc***DEBUG ends.
if (abs (test1) .eq. abs (test2)) then
if (abs (z6 - z5) .le. tol * (abs (z6) + abs (z5))) then
z6 = -z6
elseif (z6 .ne. 0.0) then
nerr = -1
endif
endif
cbugcbugc***DEBUG begins.
cbugcbug 9853 format ('DEBUG 3: z5,z6=',1p2e22.14)
cbugcbug write ( 3, 9853) z5, z6
cbugcbugc***DEBUG ends.
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptsphk results: nerr=',i3 /
cbug & ' r5, r6 = ',1p2e22.14 /
cbug & ' x1, y1, z1 =',1p3e22.14 /
cbug & ' x2, y2, z2 =',1p3e22.14 /
cbug & ' x3, y3, z3 =',1p3e22.14 /
cbug & ' x4, y4, z4 =',1p3e22.14 /
cbug & ' x5, y5, z5 =',1p3e22.14 /
cbug & ' x6, y6, z6 =',1p3e22.14 )
cbug write ( 3, 9902) nerr, r5, r6, x1, y1, z1, x2, y2, z2,
cbug & x3, y3, z3, x4, y4, z4, x5, y5, z5, x6, y6, z6
cbugc***DEBUG ends.
return
c.... End of subroutine aptsphk. (+1 line)
end
UCRL-WEB-209832