subroutine aptspsp (ra, ax, ay, az, rb, bx, by, bz, np, tol, & rc, cx, cy, cz, abx, aby, abz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTSPSP c c call aptspsp (ra, ax, ay, az, rb, bx, by, bz, np, tol, c & rc, cx, cy, cz, abx, aby, abz, nerr) c c Version: aptspsp Updated 1990 December 3 16:20. c aptspsp Originated 1990 March 20 14:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for each of np sets of input data, the radius rc and c center c = (cx, cy, cz) of the circle of intersection of the c sphere of radius ra at point a = (ax, ay, az) and the sphere of c radius rb at point b = (bx, by, bz), if an intersection occurs. c Vector ab = (abx, aby, abz) is normal to the plane of the c circle. Flag nerr indicates any input error. c c Input: ra, ax, ay, az, rb, bx, by, bz, np, tol. c c Output: rc, cx, cy, cz, abx, aby, abz, nerr. c c Calls: aptvdis, aptvadd, aptvuna c c c Glossary: c c abx,y,z Output The x, y, z components of the unit vector "ab", normal c to any plane containing any circle of intersection of c the two spheres. In the direction from point "a" c to point "b". Zero if points "a" and "b" coincide, c within the limit of precision tol. Size np. c c ax,ay,az Input The x, y, z coordinates of point "a" at the center c of the sphere with radius ra. Size np. c c bx,by,bz Input The x, y, z coordinates of point "b" at the center c of the sphere with radius rb. Size np. c c cx,cy,cz Output The x, y, z coordinates of point "c" at the center c of the circle with radius rc, at the intersection c of the two spheres, if an intersection occurs. c Size np. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c c np Input Size of arrays. c c ra Input The radius of the sphere centered at point "a". c Size np. Absolute value will be used. c c rb Input The radius of the sphere centered at point "b". c Size np. Absolute value will be used. c c rc Output The radius of the circle centered at point "c", at the c intersection of the two spheres, if an intersection c occurs. Size np. Positive if the spheres intersect. c Zero if the spheres are tangent. Negative if the c spheres do not intersect. 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. c---- Component x of vector "ab". dimension abx (1) c---- Component y of vector "ab". dimension aby (1) c---- Component z of vector "ab". dimension abz (1) c---- Coordinate x of point "a". dimension ax (1) c---- Coordinate y of point "a". dimension ay (1) c---- Coordinate z of point "a". dimension az (1) c---- Coordinate x of point "b". dimension bx (1) c---- Coordinate y of point "b". dimension by (1) c---- Coordinate z of point "b". dimension bz (1) c---- Coordinate x of point "c". dimension cx (1) c---- Coordinate y of point "c". dimension cy (1) c---- Coordinate z of point "c". dimension cz (1) c---- Radius of sphere at point "a". dimension ra (1) c---- Radius of sphere at point "b". dimension rb (1) c---- Radius of circle at point "c". dimension rc (1) c.... Local variables. c---- Local value of abx. common /laptspsp/ abxx (64) c---- Factor fact1*fact2*fact3*fact4. common /laptspsp/ arg c---- Magnitude of vector "ab". common /laptspsp/ dab (64) c---- Fractional distance of "c" along "ab". common /laptspsp/ dac (64) c---- Estimated error in factor of dac. common /laptspsp/ dacerr c---- Factor dab + rpr. common /laptspsp/ fact1 c---- Factor dab - rmr. common /laptspsp/ fact2 c---- Factor rpr - dab common /laptspsp/ fact3 c---- Factor rmr + dab. common /laptspsp/ fact4 c---- Estimated error in fact2, 3, 4. common /laptspsp/ ferr c---- A very small number. common /laptspsp/ fuz c---- Index in arrays. common /laptspsp/ n c---- First index of subset of data. common /laptspsp/ n1 c---- Last index of subset of data. common /laptspsp/ n2 c---- Index in external array. common /laptspsp/ nn c---- Size of current subset of data. common /laptspsp/ ns c---- Local value of ra. common /laptspsp/ raa (64) c---- Local value of rb. common /laptspsp/ rbb (64) c---- Local value of rc. common /laptspsp/ rcc (64) c---- Difference abs(rb) - abs(ra). common /laptspsp/ rmr c---- Sum abs(rb) + abs(ra). common /laptspsp/ rpr c---- Length of a vector. common /laptspsp/ vlen (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptspsp finding intersection between spheres:') cbug 9902 format (i3,' ra,rb= ',1p2e22.14 / cbug & ' ax,ay,az=',1p3e22.14 / cbug & ' bx,by,bz=',1p3e22.14) cbug write ( 3, 9901) cbug write ( 3, 9902) (n, ra(n), rb(n), ax(n), ay(n), az(n), cbug & bx(n), by(n), bz(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. c---- A very small number. fuz = 1.e-99 nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif n1 = 1 n2 = min (np, 64) 110 ns = n2 - n1 + 1 c.... Find the unit vector "ab" from point "a" to point "b". c.... The magnitude dab may be truncated to zero. call aptvdis (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), ns, & tol, abx(n1), aby(n1), abz(n1), dab, nerr) call aptvuna (abx(n1), aby(n1), abz(n1), ns, tol, dab, nerr) c.... Gather external arrays into local temporary arrays. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 raa(n) = ra(nn) rbb(n) = rb(nn) abxx(n) = abx(nn) c---- End of loop over subset of data. 120 continue c.... Find the radius of the circle of intersection, and the distance c.... of its center from point "a". c---- No truncation tests. if (tol .le. 0.0) then c---- Loop over subset of data. do 130 n = 1, ns rpr = abs (rbb(n)) + abs (raa(n)) rmr = abs (rbb(n)) - abs (raa(n)) dac(n) = 0.5 * (dab(n)**2 - rpr * rmr) / (dab(n) + fuz) fact1 = dab(n) + rpr fact2 = dab(n) - rmr fact3 = rpr - dab(n) fact4 = rmr + dab(n) arg = fact1 * fact2 * fact3 * fact4 rcc(n) = 0.5 * sign (sqrt (abs (arg)), arg) / (dab(n) + fuz) if ((dab(n) .le. 0.0) .and. (rmr .eq. 0.0)) then rcc(n) = 0.5 * rpr endif if (dab(n) .le. 0.0) then abxx(n) = 0.0 endif c---- End of loop over subset of data. 130 continue c---- Test for truncation. else c---- Loop over subset of data. do 140 n = 1, ns rpr = abs (rbb(n)) + abs (raa(n)) rmr = abs (rbb(n)) - abs (raa(n)) if (abs (rmr) .lt. tol * rpr) then rmr = 0.0 endif dac(n) = dab(n)**2 - rpr * rmr dacerr = tol * (dab(n)**2 + rpr * abs (rmr)) if (abs (dac(n)) .lt. dacerr) then dac(n) = 0.0 endif dac(n) = 0.5 * dac(n) / (dab(n) + fuz) fact1 = dab(n) + rpr ferr = tol * fact1 fact2 = dab(n) - rmr if (abs (fact2) .lt. ferr) then fact2 = 0.0 endif fact3 = rpr - dab(n) if (abs (fact3) .lt. ferr) then fact3 = 0.0 endif fact4 = rmr + dab(n) if (abs (fact4) .lt. ferr) then fact4 = 0.0 endif arg = fact1 * fact2 * fact3 * fact4 rcc(n) = 0.5 * sign (sqrt (abs (arg)), arg) / (dab(n) + fuz) if ((dab(n) .le. tol) .and. & (abs (rmr) .le. tol) ) then rcc(n) = 0.5 * rpr endif if (dab(n) .le. tol) then abxx(n) = 0.0 endif c---- End of loop over subset of data. 140 continue c---- Tested tol. endif c.... Scatter local temporary arrays into external arrays. c---- Loop over subset of data. do 150 n = 1, ns nn = n + n1 - 1 rc(nn) = rcc(n) abx(nn) = abxx(n) c---- End of loop over subset of data. 150 continue c.... Find the coordinates of point "c", the center of the circle c.... of intersection. c = a + dac * ab = b - fbc * ab. call aptvadd (ax(n1), ay(n1), az(n1), 1., dac, & abx(n1), aby(n1), abz(n1), ns, tol, & cx(n1), cy(n1), cz(n1), vlen, nerr) 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 cbugc***DEBUG begins. cbug 9903 format (/ 'aptspsp results:' / cbug & (i3,' rc= ',1pe22.14 / cbug & ' cx,cy,cz=',1p3e22.14 / cbug & ' abx,y,z= ',1p3e22.14)) cbug write ( 3, 9903) (n, rc(n), cx(n), cy(n), cz(n), cbug & abx(n), aby(n), abz(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptspsp. (+1 line.) end UCRL-WEB-209832