subroutine aptspha (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, np, tol, px, py, pz, rad, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPHA
c
c call aptspha (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, px, py, pz, rad, nerr)
c
c Version: aptspha Updated 1995 June 12 17:20.
c aptspha Originated 1995 June 12 17:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For each of the np sets of input data, find the sphere whose
c center is on the axis through point a = (ax, ay, az) with
c direction vector b = (bx, by, bz), and whose surface passes
c through points c = (cx, cy, cz) and d = (dx, dy,dz), and
c return the center p = (px, py, pz) and radius rad of the sphere.
c
c There is no solution if the axis vector b is perpendicular to
c the vector from point c to point d, within the estimated error
c of the calculation, based on tol.
c
c Flag nerr indicates any input error, if not zero.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, tol.
c
c Output: px, py, pz, rad, nerr.
c
c Calls: aptlnpl, aptvdis, aptvsum
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a". Size np.
c
c bx,by,bz Input The x, y, z components of vector "b". Size np.
c
c cx,cy,cz Input The x, y, z coordinates of point "c". Size np.
c
c dx,dy,dz Input The x, y, z coordinates of point "d". 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 input data (a, b, c, d) for which
c the sphere is to be found. Must be positive.
c
c px,py,pz Output The x, y, z coordinates of the center of the sphere.
c Each 1.e99 if the vector from point c to point d is
c perpendicular to the axis vector b. Size np.
c
c rad Output The radius of the sphere. Equal to 1.e99 * sqrt (3) if
c the vector from point c to point d is perpendicular
c to the axis vector b. 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.
dimension ax (1) ! Coordinate x of input point "a".
dimension ay (1) ! Coordinate y of input point "a".
dimension az (1) ! Coordinate z of input point "a".
dimension bx (1) ! Component x of input vector "b".
dimension by (1) ! Component y of input vector "b".
dimension bz (1) ! Component z of input vector "b".
dimension cx (1) ! Coordinate x of input point "c".
dimension cy (1) ! Coordinate y of input point "c".
dimension cz (1) ! Coordinate z of input point "c".
dimension dx (1) ! Coordinate x of input point "d".
dimension dy (1) ! Coordinate y of input point "d".
dimension dz (1) ! Coordinate z of input point "d".
dimension px (1) ! Coordinate x of center of sphere "p".
dimension py (1) ! Coordinate y of center of sphere "p".
dimension pz (1) ! Coordinate z of center of sphere "p".
dimension rad (1) ! Radius of sphere.
c.... Local variables. Floating point.
common /laptspha/ ex (64) ! Coordinate x of point e = a + b.
common /laptspha/ ey (64) ! Coordinate y of point e = a + b.
common /laptspha/ ez (64) ! Coordinate z of point e = a + b.
common /laptspha/ elen (64) ! Distance of point e from origin.
common /laptspha/ fx (64) ! Coordinate x of point f = .5 * (c + d).
common /laptspha/ fy (64) ! Coordinate y of point f = .5 * (c + d).
common /laptspha/ fz (64) ! Coordinate z of point f = .5 * (c + d).
common /laptspha/ flen (64) ! Distance of point f from origin.
common /laptspha/ cdx (64) ! Component x of vector cd = d - c.
common /laptspha/ cdy (64) ! Component y of vector cd = d - c.
common /laptspha/ cdz (64) ! Component z of vector cd = d - c.
common /laptspha/ cdlen (64) ! Length of vector cd.
common /laptspha/ cpx (64) ! Component x of vector cp = p - c.
common /laptspha/ cpy (64) ! Component y of vector cp = p - c.
common /laptspha/ cpz (64) ! Component z of vector cp = p - c.
common /laptspha/ dpmin (64) ! Unused array.
common /laptspha/ dint (64) ! Unused array.
common /laptspha/ fracap (64) ! Unused array.
common /laptspha/ ipar (64) ! If not zero, no center point was found.
common /laptspha/ big ! A big number.
c.... Local variables. Integers.
common /laptspha/ n ! Index in arrays
common /laptspha/ n1 ! First index of subset of data.
common /laptspha/ n2 ! Last index of subset of data.
common /laptspha/ nn ! Index in external array.
common /laptspha/ ns ! Size of current subset of data.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptspha finding sphere with its center on the',
cbug & ' axis through point a' /
cbug & ' in direction b, with its surface passing through points',
cbug & ' c and d.' /
cbug & (i4,' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' dx,dy,dz=',1p3e22.14))
cbug write ( 3, 9901) (n, ax(n), ay(n), az(n), bx(n), by(n), bz(n),
cbug & cx(n), cy(n), cz(n), dx(n), dy(n), dz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
big = 1.e99 ! A big number.
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 each axis point "e" (e = a + b).
call aptvsum (0, 1.0, ax(n1), ay(n1), az(n1),
& 1.0, bx(n1), by(n1), bz(n1), ns, tol,
& ex, ey, ez, elen, nerr)
cbugcbugc***DEBUG begins.
cbugcbug 9801 format ('ex,ey,ez= ',1p3e22.14)
cbugcbug write ( 3, 9801) (ex(n), ey(n), ez(n), n = 1, ns)
cbugcbugc***DEBUG ends.
c.... Find each midpoint vector "f".
call aptvsum (0, 0.5, cx(n1), cy(n1), cz(n1),
& 0.5, dx(n1), dy(n1), dz(n1), ns, tol,
& fx, fy, fz, flen, nerr)
cbugcbugc***DEBUG begins.
cbugcbug 9802 format ('fx,fy,fz= ',1p3e22.14)
cbugcbug write ( 3, 9802) (fx(n), fy(n), fz(n), n = 1, ns)
cbugcbugc***DEBUG ends.
c.... Find each difference vector "cd".
call aptvdis (cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1), ns,
& tol, cdx, cdy, cdz, cdlen, nerr)
cbugcbugc***DEBUG begins.
cbugcbug 9803 format ('cdx,cdy,cdz= ',1p3e22.14)
cbugcbug write ( 3, 9803) (cdx(n), cdy(n), cdz(n), n = 1, ns)
cbugcbugc***DEBUG ends.
c.... Find the center of each sphere.
call aptlnpl (ax(n1), ay(n1), az(n1), ex, ey, ez,
& fx, fy, fz, cdx, cdy, cdz, ns, tol,
& dpmin, dint, fracap, px(n1), py(n1), pz(n1),
& ipar, nerr)
c.... Find the radius of each sphere.
call aptvdis (cx(n1), cy(n1), cz(n1), px(n1), py(n1), pz(n1), ns,
& tol, cpx, cpy, cpz, rad(n1), nerr)
c.... Set the radius to a big number if no unique solution is found.
do 120 n = 1, ns
if (ipar(n) .ne. 0) then
nn = n + n1 - 1
px(nn) = big
py(nn) = big
pz(nn) = big
rad(nn) = big * sqrt (3.0)
endif
120 continue
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
c.... Done.
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptspha results. nerr=',i3 /
cbug & (i4,' px,py,pz=',1p3e22.14 /
cbug & ' rad= ',1pe22.14))
cbug write ( 3, 9902) nerr, (n, px(n), py(n), pz(n),
cbug & rad(n), n = 1, np)
cbugc***DEBUG ends.
return
c.... End of subroutine aptspha. (+1 line.)
end
UCRL-WEB-209832