subroutine aptsphl (r1, r2, r3, r4, tol, x4, y4, z4, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPHL
c
c call aptsphl (r1, r2, r3, r4, tol, x4, y4, z4, nerr)
c
c Version: aptsphl Updated 2001 June 13 13:20.
c aptsphl Originated 1999 January 28 11:30.
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, with the center of the first sphere at the origin, the
c center of the second sphere at coordinate r1 + r2 on the x axis,
c and the center of the third sphere at
c x3 = r1 + r3 * (r1 - r2) / (r1 + r2)
c y3 = 2.0 * sqrt (r1*r2*r3*(r1+r2+r3)) / abs (r1 + r2)
c z3 = 0,
c unless r1 + r2 + r3 = 0, and r4 has a unique required value,
c r4 = 1 / sqrt (1/r1**2 + 1/r2**2 + 1/r3**2)
c in which case:
c x3 = r2
c y3 = 0
c z3 = 0,
c find the center point p4 = (x4, y4, z4) of the fourth sphere.
c Flag nerr indicates any input or result error, if not zero.
c
c Method: The center of the fourth sphere, point p4 = (x4, y4, z4) is at:
c
c x4 = r1 + r4 * (r1 - r2) / (r1 + r2)
c y4 = 2 * (r1*r2*(r1+r2)*(r3+r4)-r3*r4*(r1**2+r2**2)) /
c (y3*(r1+r2)**2)
c y4**2 + z4**2 = 4.0 * r1*r2*r4*(r1+r2+r4) / (r1 + r2)**2
c z4 = sqrt ((y4**2 + z4**2) - y4**2)
c unless r1 + r2 + r3 = 0, and r4 has a unique required value,
c r4 = 1 / sqrt (1/r1**2 + 1/r2**2 + 1/r3**2),
c in which case:
c x4 = r1 + r4 * (r1 - r2) / (r1 + r2).
c y4 = 2.0 * r4
c z4 = 0.
c
c Input: r1, r2, r3, r4, tol.
c
c Output: x4, y4, z4, nerr.
c
c Calls: aptcirl
c
c Glossary:
c
c nerr Output Indicates an input or result error, if not 0.
c -1 if the magnitude of z4 is small but uncertain.
c A zero value of z4 is returned.
c 1 if two or more of the radii are negative, or if
c one radius is negative, but not large enough in
c magnitude for that sphere to enclose the other
c three, or one radius is too small for that sphere
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. A negative
c radius indicates a sphere that surrounds the other
c three spheres. A very large radius, relative to the
c other three, indicates a straight line.
c
c x4,y4,z4 Output The x, y, z coordinates of the center of the fourth
c sphere. Returned as -1.e99 if nerr = 1.
c The sign of z4 is arbitrary. The positive value of
c z4 is returned.
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 (/ 'aptsphl finding coordinates of 4th of four',
cbug & ' tangent 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.
x4 = -1.e99
y4 = -1.e99
z4 = -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 curvature of the first four spheres.
b1 = 1.0 / r1
b2 = 1.0 / r2
b3 = 1.0 / r3
b4 = 1.0 / r4
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.
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
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
go to 210
endif ! Special case.
c.... Find the center of sphere 3.
call aptcirl (r1, r2, r3, tol, x3, y3, nerra)
if (nerra .ne. 0) then
nerr = 1
go to 210
endif
if (y3 .eq. 0.0) then
nerr = 1
go to 210
endif
z3 = 0.0
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 (r3) + abs (r1))
if (abs (d31) .le. err31) d31 = 0.0
d14 = r4 - r1
err14 = tol * (abs (r1) + abs (r4))
if (abs (d14) .le. err14) d14 = 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 the center of sphere 4 (x4, y4, z4).
c.... Find x4, and the value of argyz = sqrt (y4**2 + z4**2).
call aptcirl (r1, r2, r4, tol, x4, argyz, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9501 format (' x4, argyz =',1p2e22.14)
cbugcbug write ( 3, 9501) x4, argyz
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 1
go to 210
endif
c.... Find y4.
if ((d12 .eq. 0.0) .and.
& ((d23 .eq. 0.0) .or. (d24 .eq. 0.0))) then
y4 = r1**2 / y3
elseif (((d31 .eq. 0.0) .and. (d24 .eq. 0.0)) .or.
& ((d14 .eq. 0.0) .and. (d23 .eq. 0.0))) then
y4 = 4.0 / ((b1 + b2)**2 * y3)
elseif (d12 .eq. 0.0) then
arg = b3 + b4 - b1
err = tol * (abs (b3) + abs (b4) + abs (b1))
cbugcbugc***DEBUG begins.
cbugcbug 9701 format (' arg12, err =',1p2e22.14 )
cbugcbug write ( 3, 9701) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9701) arg, err
cbugcbugc***DEBUG ends.
y4 = arg / (b1 * b3 * b4 * y3)
elseif (d23 .eq. 0.0) then
arg = b1 * (b2 + b4 - b1) + b2 * b4
err = tol * (abs (b1 * b2) + abs (b1 * b4) +
& b1**2 + abs (b2 * b4))
cbugcbugc***DEBUG begins.
cbugcbug 9702 format (' arg23, err =',1p2e22.14 )
cbugcbug write ( 3, 9702) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9702) arg, err
cbugcbugc***DEBUG ends.
y4 = 2.0 * arg / (b2 * b4 * (b1 + b2)**2 * y3)
elseif (d31 .eq. 0.0) then
arg = b2 * (b1 + b4 - b2) + b1 * b4
err = tol * (abs (b1 * b2) + abs (b2 * b4) +
& b2**2 + abs (b1 * b4))
cbugcbugc***DEBUG begins.
cbugcbug 9703 format (' arg31, err =',1p2e22.14 )
cbugcbug write ( 3, 9703) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9703) arg, err
cbugcbugc***DEBUG ends.
y4 = 2.0 * arg / (b1 * b4 * (b1 + b2)**2 * y3)
elseif (d14 .eq. 0.0) then
arg = b2 * (b1 + b3 - b2) + b1 * b3
err = tol * (abs (b1 * b2) + abs (b2 * b3) +
& b2**2 + abs (b1 * b3))
cbugcbugc***DEBUG begins.
cbugcbug 9704 format (' arg14, err =',1p2e22.14 )
cbugcbug write ( 3, 9704) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9704) arg, err
cbugcbugc***DEBUG ends.
y4 = 2.0 * arg / (b1 * b3 * (b1 + b2)**2 * y3)
elseif (d24 .eq. 0.0) then
arg = b1 * (b2 + b3 - b1) + b2 * b3
err = tol * (abs (b1 * b2) + abs (b1 * b3) +
& b1**2 + abs (b2 * b3))
cbugcbugc***DEBUG begins.
cbugcbug 9705 format (' arg24, err =',1p2e22.14 )
cbugcbug write ( 3, 9705) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9705) arg, err
cbugcbugc***DEBUG ends.
y4 = 2.0 * arg / (b2 * b3 * (b1 + b2)**2 * y3)
else
arg = (b1 + b2) * (b3 + b4) - (b1**2 + b2**2)
err = tol * (abs (b1 * b3) + abs (b1 * b4) +
& abs (b2 * b3) + abs (b2 * b4) +
& b1**2 + b2**2)
cbugcbugc***DEBUG begins.
cbugcbug 9706 format (' arg , err =',1p2e22.14 )
cbugcbug write ( 3, 9706) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9706) arg, err
cbugcbugc***DEBUG ends.
y4 = 2.0 * arg / (b3 * b4 * (b1 + b2)**2 * y3)
endif
c.... Find z4.
if (d12 .eq. 0.0) then
cbugcbugc***DEBUG begins.
cbugcbug 9301 format ('r1 = r2 =',1p2e22.14 )
cbugcbug 9302 format ('y3, y3**2 =',1p2e22.14 )
cbugcbug 9303 format ('y4, y4**2 =',1p2e22.14 )
cbugcbug 9304 format ('y4**2 + z4**2=',1p2e22.14 )
cbugcbug 9305 format ('z4**2, z4**2 =',1p2e22.14 )
cbugcbug 9306 format ('z4a , z4b =',1p2e22.14 )
cbugcbug write ( 3, 9301) r1, r2
cbugcbug y3a = sqrt (r3 * (2.0 * r1 + r3))
cbugcbug y3a2 = r3 * (2.0 * r1 + r3)
cbugcbug write ( 3, 9302) y3a, y3a2
cbugcbug y4a = (r1 * (r3 + r4) - r3 * r4) / sqrt (r3 * (2.0 * r1 + r3))
cbugcbug y4a2 = (r1 * (r3 + r4) - r3 * r4)**2 / (r3 * (2.0 * r1 + r3))
cbugcbug write ( 3, 9303) y4a, y4a2
cbugcbug ssq = r4 * (2.0 * r1 + r4)
cbugcbug write ( 3, 9304) ssq
cbugcbug z4a2 = ssq - y4a2
cbugcbug z4b2 = r1 * (r3 + r4) * (4.0 * r3 * r4 + r1 * (r4 - r3)) /
cbugcbug & (r3 * (2.0 * r1 + r3))
cbugcbug write ( 3, 9305) z4a2, z4b2
cbugcbug z4a = -1.e99
cbugcbug z4b = -1.e99
cbugcbug if (z4a2 .ge. 0.0) z4a = sqrt (z4a2)
cbugcbug if (z4b2 .ge. 0.0) z4b = sqrt (z4b2)
cbugcbug write ( 3, 9306) z4a, z4b
cbugcbugc***DEBUG ends.
if (d34 .eq. 0.0) then
arg = 8.0 * r1 * r3**2 / (2.0 * r1 + r3)
err = 0.0
else
arg = r1 * (4.0 * r3 * r4 * (r3 + r4) - r1 * (r3 - r4)**2) /
& (r3 * (2.0 * r1 + r3))
a3 = abs (r3)
a4 = abs (r4)
err = tol * r1 * (4.0 * a3 * a4 * (a3 + a4) +
& r1 * (a3 + a4)**2) /
& (a3 * (2.0 * r1 + a3))
cold ssq = r4 * (2.0 * r1 + r4)
cold y4a2 = (r1 * (r3 + r4) - r3 * r4)**2 / (r3 * (2.0 * r1 + r3))
cold arg = ssq - y4a2
cbugcbugc***DEBUG begins.
cbugcbug adiffg = 4.0 * r1 * r3 * r4 * (r3 + r4) - r1**2 * (r3 - r4)**2
cbugcbug aden = r3 * (2.0 * r1 + r3)
cbugcbug arg2 = adiffg / aden
cbugcbug
cbugcbug 9609 format ('arg, arg2 =',1p2e22.14 )
cbugcbug write ( 3, 9609) arg, arg2
cbugcbugc***DEBUG ends.
err = 0.0
endif
else
arg = argyz**2 - y4**2
err = tol**2 * (argyz**2 + y4**2)
endif
cbugcbugc***DEBUG begins.
cbugcbug arg1 = 8.0 * r1 * r3**2 / (2.0 * r1 + r3)
cbugcbug arg2 = r1 * (4.0 * r3 * r4 * (r3 + r4) - r1 * (r3 - r4)**2) /
cbugcbug & (r3 * (2.0 * r1 + r3))
cbugcbug arg3 = argyz**2 - y4**2
cbugcbug 9708 format (' arg1,2,3 =',1p3e22.14 )
cbugcbug write ( 3, 9708) arg1, arg2, arg3
cbugcbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug 9707 format (' argz , err =',1p2e22.14 )
cbugcbug write ( 3, 9707) arg, err
cbugcbugc***DEBUG ends.
if (abs (arg) .le. err) arg = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9707) arg, err
cbugcbugc***DEBUG ends.
if (arg .ge. 0.0) then
z4 = sqrt (arg)
else
nerr = -1
z4 = 0.0
endif
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptsphl results: nerr=',i3 /
cbug & ' x4, y4, z4 =',1p3e22.14 )
cbug write ( 3, 9902) nerr, x4, y4, z4
cbugc***DEBUG ends.
return
c.... End of subroutine aptsphl. (+1 line.)
end
UCRL-WEB-209832