subroutine aptwhis (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, & qxx, qyy, qzz, tol, dside, bx, by, bz, & nprox, cx, cy, cz, dist, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTWHIS c c call aptwhis (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz, tol, c dside, bx, by, bz, nprox, cx, cy, cz, dist, nerr) c c Version: aptwhis Updated 1998 June 23 17:30. c aptwhis Originated 1994 June 14 15:10. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for point "a" = (ax, ay, az), and the quadric surface c represented by f(x,y,z) = 0, which side of the quadric surface c point "a" is on, by finding dside = f(ax,ay,az), where: c c f(x,y,z) = qc + qx * x + qy * y + qz * z + c qxy * x * y + qyz * y * z + qzx * z * x + c qxx * x**2 + qyy * y**2 + qzz * z**2, c c and to find the value at point "a" of the normal vector c "b" = (bx,by,bz) = grad f: c c bx = qx + 2.0 * qxx * x + qxy * y + qzx * z c by = qy + qxy * x + 2.0 * qyy * y + qyz * z c bz = qz + qzx * x + qyz * y + 2.0 * qzz * z, c c and to find the estimated (nprox = -1 or 0, rarely) or accurate c (nprox = 1, 2, 3 or 4) point or points "c" on the quadric c surface nearest to point "a", and the distance dist of point c or points "c" from point "a". For a circle of proximal c points (nprox = 3), three representative points are returned. c For a sphere of proximal points (nprox = 4), four representative c points are returned. c c One estimate of dist is -dside / |b|. c Another estimate (nprox = 0) is any intersection "c" of the c straight line through point "a" with direction vector "b" c with the quadric surface, at distance dist. c c Any result less than the estimated error in its calculation, c based on tol, will be truncated to zero. c c Input: ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz, tol. c c Output: dside, bx, by, bz, nprox, cx, cy, cz, dist, nerr. c c Calls: aptaxis, aptcubs, aptmopv, aptqrts, aptqnor, aptquar, aptqval, c apttran, aptvunz c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a". c c bx,by,bz Output The x, y, z components of the normal vector "b", the c gradient of function f(x,y,z), at point "a". c c cx,cy,cz Output If nprox = 1, 2, 3, or 4, the x, y, z coordinates of c the nprox points "c" on the quadric surface nearest c point "a". If nprox = 0, point "c(1)" is the c intersection of the line through point "a" with c direction vector "b" with the quadric surface, at c distance dist(1) from "a". Array size = 4. c c dist Output If nprox = 1, 2, 3, or 4, the nprox distances of the c nprox proximal points "c" from point "a". c If nprox = 0 (rarely), dist(1) is the distance from c "a" of the intersection of the line through point "a" c with direction vector "b" with the quadric surface. c Positive if in the same hemisphere from point "a" as c vector "b". Array size = 4. c c dside Output The value of function f(x,y,z) at point "a", c f(ax,ay,az). Positive if on the side of the quadric c surface where the normal vector points away from the c surface. c c nerr Output Indicates an input or method error, if not 0. c 1 if all coefficients are zero except possibly qc. c 2 if the quadric surface is imaginary. c 3 if the method fails. coldc 4 if the calculated proximal point is not on the coldc quadric surface. c c nprox Output The number of accurate proximal points "c" and proximal c distances dist found. c -1 if none (rarely), and the line from "a" in c direction "b" does not intersect the quadric c surface. An arbitrary point "c(1)" on the quadric c surface is returned, at distance dist(1) from point c "a". If nerr = 4, the method failed. c 0 if none (rarely), and the line from "a" in the c normal direction "b" intersects the quadric surface c at point "c(1)", distance dist(1). c 1 if a single accurate proximal point is found. c 2 if two equally distant proximal points are found. c 3 if point "a" is on the axis of a circular cylinder c or circular cone, so the proximal points form a c circle. Three arbitrary points "c", at 120 degree c intervals around the circle, will be returned. c 4 if point "a" is at the center of a sphere, so the c proximal points are the entire sphere. Four c arbitrary points "c", in the -x, x, y and z c directions from the center of the sphere, will be c returned. c c q.. Input Coefficients of the implicit equation of a second-order c surface in xyz space (qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz). c c tol Input Numerical tolerance limit. c On computers with 64-bit floating point words, c recommend 1.e-5 to 1.e-15. c c History: 1997 June 12 14:00. Replaced coding for hyperbolic c cylinder, to correct errors. c c 1997 September 23 14:00. Finished adding coding to c find accurate proximal points for all quadric c surfaces. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. dimension cx (4) dimension cy (4) dimension cz (4) dimension dist (4) c.... Local variables. common /laptwhis/ a(0:6) ! Coefficients of equation for f. common /laptwhis/ atype(3) character*8 atype common /laptwhis/ crtest(3) common /laptwhis/ cxtest(3) common /laptwhis/ cytest(3) common /laptwhis/ cztest(3) common /laptwhis/ distest(3) common /laptwhis/ fuz ! A very small number. common /laptwhis/ rreal(4) common /laptwhis/ rimag(4) common /laptwhis/ rotm (3,3) common /laptwhis/ xreal(3) common /laptwhis/ ximag(3) cbugc***DEBUG begins. cbug 9901 format (/ 'aptwhis finding the side of point "a" relative' / cbug & ' to the quadric surface with the coefficients qc-qzz:' / cbug & ' tol =',1pe22.14 / cbug & ' ax, ay, az =',1p3e22.14 / cbug & ' qc =',1pe22.14 / cbug & ' qx, qy, qz =',1p3e22.14 / cbug & ' qxy,qyz,qzx=',1p3e22.14 / cbug & ' qxx,qyy,qzz=',1p3e22.14 ) cbug write ( 3, 9901) tol, ax, ay, az, cbug & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz cbugc***DEBUG ends. c.... Initialize. nerr = 0 nprox = 0 fuz = 1.e-99 do 10 n = 1, 4 cx(n) = -1.e99 cy(n) = -1.e99 cz(n) = -1.e99 dist(n) = -1.e99 10 continue c.... Find the standard form of the quadric surface. sc = qc sx = qx sy = qy sz = qz sxy = qxy syz = qyz szx = qzx sxx = qxx syy = qyy szz = qzz tolx = tol if (tolx .le. 0.0) tolx = 1.e-11 call aptaxis (sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz, & tolx, cenx, ceny, cenz, rotm, ntype, nerr) if (nerr .ne. 0) go to 710 if (ntype .gt. 16) then nerr = 2 go to 710 endif c.... Find the side of the quadric surface point "a" is on. call aptqval (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, & qxx, qyy, qzz, tol, dside) cold dside = qc + qx * ax + qy * ay + qz * az + cold & qxy * ax * ay + cold & qyz * ay * az + cold & qzx * az * ax + cold & qxx * ax**2 + qyy * ay**2 + qzz * az**2 cold cold derr = tol * (abs (qc) + cold & 2.0 * (abs (qx * ax) + abs (qy * ay) + abs (qz * az)) + cold & 3.0 * (abs (qxy * ax * ay) + cold & abs (qyz * ay * az) + cold & abs (qzx * az * ax) + cold & abs (qxx) * ax**2 + cold & abs (qyy) * ay**2 + cold & abs (qzz) * az**2 )) cbugcbugcoldc***DEBUG begins. cbugcbugcold 9941 format ('DEBUG dside, derr=',1p2e22.14) cbugcbugcold write ( 3, 9941) dside, derr cbugcbugcoldc***DEBUG ends. c.... Find the normal vector at point "a". Truncate small components. call aptqnor (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, & qxx, qyy, qzz, tol, bc, bx, by, bz, vlen, nerr) cold bx = qx + 2.0 * qxx * ax + qxy * ay + qzx * az cold by = qy + qxy * ax + 2.0 * qyy * ay + qyz * az cold bz = qz + qzx * ax + qyz * ay + 2.0 * qzz * az cold cold call aptvunz (bx, by, bz, tol, vlen) cold cold bx = vlen * bx cold by = vlen * by cold bz = vlen * bz if (dside .eq. 0.0) then nprox = 1 ! Exact. cx(1) = ax cy(1) = ay cz(1) = az dist(1) = 0.0 dside = 0.0 go to 710 endif c.... Find the estimated proximal distance from point "a" to the surface. nprox = -1 c.... Find any intersection "c" of the track through point "a" with direction c.... vector "b" and the quadric surface, and its distance from point "a". aa = qxx * bx**2 + qyy * by**2 + qzz * bz**2 + & qxy * bx * by + qyz * by * bz + qzx * bz * bx bb = bx**2 + by**2 + bz**2 cc = dside call aptqrts (0, aa, bb, cc, qq, tol, nroots, root1, root2, & itrun) if (nroots .ge. 1) then ! Normal track intersects surface. nprox = 0 root = root1 if ((nroots .eq. 2) .and. & (abs (root2) .lt. abs (root1))) then root = root2 endif cx(1) = ax + root * bx cy(1) = ay + root * by cz(1) = az + root * bz dist(1) = root * vlen endif ! Normal track intersects surface. c.... Find the proximal point and distance from point "a" to the c.... quadric surface. c.... Find point "aa", the image of point "a" in standard coordinates. aax = ax aay = ay aaz = az call apttran (cenx, ceny, cenz, aax, aay, aaz, 1, tol, nerra) call aptmopv (rotm, 0, 0., 0., 0., aax, aay, aaz, 1, tol, nerra) c.... Move point "aa" onto axis if distance from axis negligible. rsq = aax**2 + aay**2 + aaz**2 if (aax**2 .le. tol * rsq) aax = 0.0 if (aay**2 .le. tol * rsq) aay = 0.0 if (aaz**2 .le. tol * rsq) aaz = 0.0 cbugcbugc***DEBUG begins. cbugcbug cbugcbug 9902 format (/ ' quadric ntype =',i3 / cbugcbug & ' cenx,y,z =',1p3e22.14 ) cbugcbug 9903 format ( cbugcbug & ' x'' axis ',1p3e22.14 / cbugcbug & ' y'' axis ',1p3e22.14 / cbugcbug & ' z'' axis ',1p3e22.14) cbugcbug 9904 format (/ 'aptwhis: standard coordinates:' / cbugcbug & ' aax,aay,aaz=',1p3e22.14 / cbugcbug & ' sc =',1pe22.14 / cbugcbug & ' sx, sy, sz =',1p3e22.14 / cbugcbug & ' sxy,syz,szx=',1p3e22.14 / cbugcbug & ' sxx,syy,szz=',1p3e22.14 ) cbugcbug write ( 3, 9902) ntype, cenx, ceny, cenz cbugcbug write ( 3, 9903) ((rotm(i,j), j = 1, 3), i = 1, 3) cbugcbug write ( 3, 9904) aax, aay, aaz, sc, sx, sy, sz, sxy, syz, szx, cbugcbug & sxx, syy, szz cbugcbugc***DEBUG ends. c.... Use separate coding for each type of quadric surface. nntype = ntype + 1 go to (100, 100, 102, 103, 104, 105, 106, 107, 108, 109, & 110, 111, 112, 113, 114, 115, 116), nntype c=======================================================================******** c.... Type 0. Simple plane: x = 0. c.... Type 1. Coincident planes: x**2 = 0. 100 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'SIMPLE PLANE OR COINCIDENT PLANES.')") cbugcbugc***DEBUG ends. cx(1) = 0.0 cy(1) = aay cz(1) = aaz dist(1) = - aax go to 510 c.... Type 0. Simple plane: x = 0. c.... Type 1. Coincident planes: x**2 = 0. c=======================================================================******** c.... Type 2. Real parallel planes: -|sc| + sxx * x**2 = 0. 102 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'PARALLEL PLANES.')") cbugcbugc***DEBUG ends. xp = sqrt (-sc / sxx) if (aax .lt. 0.0) then ! Nearer negative plane. cx(1) = -xp cy(1) = aay cz(1) = aaz dist(1) = -xp - aax elseif (aax .eq. 0.0) then ! Equidistant between planes. nprox = 2 ! Exact. cx(1) = -xp cy(1) = aay cz(1) = aaz dist(1) = -xp - aax cx(2) = xp cy(2) = aay cz(2) = aaz dist(2) = xp - aax else ! Nearer positive plane. cx(1) = xp cy(1) = aay cz(1) = aaz dist(1) = xp - aax endif go to 510 c.... Type 2. Real parallel planes: -|sc| + sxx * x**2 = 0. c=======================================================================******** c.... Type 3. Real intersecting planes: sxx * x**2 - |syy| * y**2 = 0. 103 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'INTERSECTING PLANES.')") cbugcbugc***DEBUG ends. crat = -syy / sxx cratp = 1.0 + crat cratr = sqrt (crat) if (aax * aay .lt. 0.0) then ! Point "a" in quadrant 2 or 4. cy(1) = (aay - cratr * aax) / cratp cx(1) = -cratr * cy(1) cz(1) = aaz dist(1) = -(aax + cratr * aay) / sqrt (cratp) elseif (aax * aay .gt. 0.0) then ! Point "a" in quadrant 1 or 3. cy(1) = (aay + cratr * aax) / cratp cx(1) = cratr * cy(1) cz(1) = aaz dist(1) = (aax - cratr * aay) / sqrt (cratp) else ! Midway between planes. nprox = 2 ! Exact. cy(1) = (aay - cratr * aax) / cratp cx(1) = -cratr * cy(1) cz(1) = aaz dist(1) = (aax + cratr * aay) / sqrt (cratp) cy(2) = (aay + cratr * aax) / cratp cx(2) = cratr * cy(2) cz(2) = aaz dist(2) = (aax - cratr * aay) / sqrt (cratp) endif ! Tested symmetry condition. go to 510 c.... Type 3. Real intersecting planes: sxx * x**2 - |syy| * y**2 = 0. c=======================================================================******** c.... Type 4. Parabolic cylinder: (+/-)su * u + sxx * x**2 = 0 (u = y or z). 104 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'PARABOLIC CYLINDER.')") cbugcbugc***DEBUG ends. if (aax .eq. 0.0) then ! One or two proximal points. arg = -(sy * (aay + 0.5 * sy / sxx) + & sz * (aaz + 0.5 * sz / sxx)) / sxx if (arg .le. 0.0) then ! One proximal point at vertex. cx(1) = 0.0 if (sy .eq. 0.0) then cy(1) = aay cz(1) = 0.0 dist(1) = aaz else cy(1) = 0.0 cz(1) = aaz dist(1) = aay endif else ! Two proximal points. nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = aay + 0.5 * sy / sxx cz(1) = aaz + 0.5 * sz / sxx cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(1) = sqrt (arg + (0.5 * sy / sxx)**2 + & (0.5 * sz / sxx)**2) dist(2) = dist(1) endif ! One or two proximal points. else ! Possible three proximal points. c.... Solve cubic for one to three values of cx. aacub = 0.0 bbcub = (sy * (aay + 0.5 * sy / sxx) + & sz * (aaz + 0.5 * sz / sxx)) / sxx cccub = -0.5 * (sy**2 + sz**2) * aax / sxx**2 call aptcubs (cccub, bbcub, aacub, tol, nroots, xreal, ximag, & atype, itrun) if (nroots .ge. 1) then cxtest(1) = xreal(1) if (sy .eq. 0.0) then cytest(1) = aay cztest(1) = -sxx * cxtest(1)**2 / sz distest(1) = sqrt ((cxtest(1) - aax)**2 + & (cztest(1) - aaz)**2) else cytest(1) = -sxx * cxtest(1)**2 / sy cztest(1) = aaz distest(1) = sqrt ((cxtest(1) - aax)**2 + & (cytest(1) - aay)**2) endif nuse = 1 endif if (nroots .ge. 2) then cxtest(2) = xreal(2) if (sy .eq. 0.0) then cytest(2) = aay cztest(2) = -sxx * cxtest(2)**2 / sz distest(2) = sqrt ((cxtest(2) - aax)**2 + & (cztest(2) - aaz)**2) else cytest(2) = -sxx * cxtest(2)**2 / sy cztest(2) = aaz distest(2) = sqrt ((cxtest(2) - aax)**2 + & (cytest(2) - aay)**2) endif if (distest(2) .lt. distest(1)) nuse = 2 endif if (nroots .eq. 3) then cxtest(3) = xreal(3) if (sy .eq. 0.0) then cytest(3) = aay cztest(3) = -sxx * cxtest(3)**2 / sz distest(3) = sqrt ((cxtest(3) - aax)**2 + & (cztest(3) - aaz)**2) else cytest(3) = -sxx * cxtest(3)**2 / sy cztest(3) = aaz distest(3) = sqrt ((cxtest(3) - aax)**2 + & (cytest(3) - aay)**2) endif if (distest(3) .lt. distest(nuse)) nuse = 3 endif cx(1) = cxtest(nuse) cy(1) = cytest(nuse) cz(1) = cztest(nuse) dist(1) = distest(nuse) endif ! Tested for number of proximal points. go to 510 c.... Type 4. Parabolic cylinder: (+/-)su * u + sxx * x**2 = 0 (u = y or z). c=======================================================================******** c.... Type 5. Hyperbolic cylinder: (+/-)sc + sxx * x**2 - |syy| * y**2 = 0. 105 if (sc .lt. 0.0) then ! Intercepts are on x axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'HYPERBOLIC CYLINDER INTERCEPTING X AXIS.')") cbugcbugc***DEBUG ends. xint = sqrt (-sc / sxx) ! Positive intercept on x axis. xlim = xint * (1.0 - sxx / syy) ! Positive center of curvature. if ((aax .eq. 0.0) .and. (aay .eq. 0.0)) then nprox = 2 ! Exact. cx(1) = xint cy(1) = 0.0 cz(1) = aaz dist(1) = cx(1) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 elseif (aax .eq. 0.0) then ! In plane x = 0, not at y = 0. nprox = 2 ! Exact. cy(1) = aay / (1.0 - syy / sxx) cx(1) = sqrt (-(sc + syy * cy(1)**2) / sxx) cz(1) = aaz dist(1) = sqrt (cx(1)**2 + (cy(1) - aay)**2) cy(2) = cy(1) cx(2) = -cx(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 elseif (aay .eq. 0.0) then ! In plane y = 0, not at x = 0. if (aax**2 .le. xlim**2) then nprox = 1 ! Exact. At an x-intercept. cx(1) = xint * aax / abs (aax) cy(1) = 0.0 cz(1) = aaz dist(1) = abs (cx(1) - aax) else nprox = 2 ! Exact. cx(1) = aax / (1.0 - sxx / syy) cy(1) = sqrt (-(sc + sxx * cx(1)**2) / syy) cz(1) = aaz dist(1) = sqrt ((cx(1) - aax)**2 + cy(1)**2) cy(2) = -cy(1) cx(2) = cx(1) cz(2) = cz(1) dist(2) = dist(1) endif go to 510 endif ! Point "a" not at x = 0 or y = 0. else ! Intercepts are on y axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'HYPERBOLIC CYLINDER INTERCEPTING Y AXIS.')") cbugcbugc***DEBUG ends. yint = sqrt (-sc / syy) ! Positive intercept on y axis. ylim = yint * (1.0 - syy / sxx) ! Positive center of curvature. if ((aax .eq. 0.0) .and. (aay .eq. 0.0)) then nprox = 2 ! Exact. cx(1) = 0.0 cy(1) = yint cz(1) = aaz dist(1) = cy(1) cx(2) = cx(1) cy(2) = -cy(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 elseif (aay .eq. 0.0) then ! In plane y = 0, not at x = 0. nprox = 2 ! Exact. cx(1) = aax / (1.0 - sxx / syy) cy(1) = sqrt (-(sc + sxx * cx(1)**2) / syy) cz(1) = aaz dist(1) = sqrt (cy(1)**2 + (cx(1) - aax)**2) cx(2) = cx(1) cy(2) = -cy(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 elseif (aax .eq. 0.0) then ! In plane x = 0, not at y = 0. if (aay**2 .le. ylim**2) then nprox = 1 ! Exact. At a y-intercept. cx(1) = 0.0 cy(1) = yint * aay / abs (aay) cz(1) = aaz dist(1) = abs (cy(1) - aay) else nprox = 2 ! Exact. cy(1) = aay / (1.0 - syy / sxx) cx(1) = sqrt (-(sc + syy * cy(1)**2) / sxx) cz(1) = aaz dist(1) = sqrt ((cy(1) - aay)**2 + cx(1)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) endif go to 510 endif ! Point "a" not at x = 0 or y = 0. endif ! Tested sign of sc. c.... Point "a" is not at origin or at x = 0 or y = 0. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aax**2 + syy * aay**2 c1 = -4.0 * (sc * (sxx + syy) + sxx * syy * (aax**2 + aay**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * syy + syy**2) + & sxx * syy * (syy * aax**2 + sxx * aay**2)) c3 = -16.0 * sc * sxx * syy * (sxx + syy) c4 = 16.0 * sc * sxx**2 * syy**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 210 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny distsq = (cxx - aax)**2 + (cyy - aay)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = aaz dist(1) = sqrt (distsq) endif 210 continue go to 510 c.... Type 5. Hyperbolic cylinder: sc + sxx * x**2 - |syy| * y**2 = 0. c=======================================================================******** c.... Type 6. Real elliptic cylinder: -|sc| + sxx * x**2 + |syy| * y**2 = 0. c.... Use two proximal points if point "a" is on the z axis. 106 aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'ELLIPTIC CYLINDER.')") cbugcbugc***DEBUG ends. if (aar .eq. 0.0) then nprox = 2 ! Exact. cz(1) = aaz cz(2) = aaz if (sxx .gt. syy) then ! ALWAYS TRUE. cx(1) = sqrt (-sc / sxx) cy(1) = 0.0 cx(2) = -cx(1) cy(2) = cy(1) dist(1) = cx(1) else ! ALWAYS FALSE. cx(1) = 0.0 cy(1) = sqrt (-sc / syy) cx(2) = cx(1) cy(2) = -cy(1) dist(1) = cy(1) endif dist(2) = dist(1) elseif (aax .eq. 0.0) then cy(1) = aay / (1.0 - syy / sxx) arg = -(sc + syy * cy(1)**2) / sxx if ((syy .gt. sxx) .or. (arg .lt. 0.0)) then nprox = 1 ! Exact. cx(1) = 0.0 cy(1) = sqrt (-sc / syy) if (aay .lt. 0.0) cy(1) = -cy(1) cz(1) = aaz dist(1) = cy(1) - aay else nprox = 2 ! Exact. cx(1) = sqrt (arg) cz(1) = aaz dist(1) = sqrt (arg + (cy(1) - aay)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = aaz endif dist(2) = dist(1) elseif (aay .eq. 0.0) then cx(1) = aax / (1.0 - sxx / syy) arg = -(sc + sxx * cx(1)**2) / syy if ((sxx .gt. syy) .or. (arg .lt. 0.0)) then ! ALWAYS TRUE. nprox = 1 ! Exact. cx(1) = sqrt (-sc / sxx) cy(1) = 0.0 if (aax .lt. 0.0) cx(1) = -cx(1) cz(1) = aaz dist(1) = cx(1) - aax else ! ALWAYS FALSE. nprox = 2 ! Exact. cy(1) = sqrt (arg) cz(1) = aaz dist(1) = sqrt (arg + (cx(1) - aax)**2) cx(2) = cx(1) cy(2) = -cy(1) cz(2) = aaz endif dist(2) = dist(1) else ! Find general solution. c.... Point "a" is not at origin or at x = 0 or y = 0. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aax**2 + syy * aay**2 c1 = -4.0 * (sc * (sxx + syy) + sxx * syy * (aax**2 + aay**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * syy + syy**2) + & sxx * syy * (syy * aax**2 + sxx * aay**2)) c3 = -16.0 * sc * sxx * syy * (sxx + syy) c4 = 16.0 * sc * sxx**2 * syy**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 220 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny distsq = (cxx - aax)**2 + (cyy - aay)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = aaz dist(1) = sqrt (distsq) endif 220 continue endif ! Tested for special or general case. go to 510 c.... Type 6. Real elliptic cylinder: -|sc| + sxx * x**2 + |syy| * y**2 = 0. c=======================================================================******** c.... Type 7. Real circular cylinder: -|sc| + sxx * (x**2 + y**2) = 0. 107 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'CIRCULAR CYLINDER.')") cbugcbugc***DEBUG ends. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. rc = sqrt (-sc / sxx) ! Radius of cylinder. if (aar .eq. 0.0) then ! Point "a" is on the z axis. nprox = 3 ! Exact. 3 points on the proximal ring. cx(1) = rc cy(1) = 0.0 cz(1) = aaz cx(2) = -0.5 * rc cy(2) = 0.5 * sqrt (3.0) * rc cz(2) = cz(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(1) = rc dist(2) = dist(1) dist(3) = dist(1) else ! Point "a" is not on the z axis. cx(1) = aax * rc / aar cy(1) = aay * rc / aar cz(1) = aaz dist(1) = rc - aar endif go to 510 c.... Type 7. Real circular cylinder: -|sc| + sxx * (x**2 + y**2) = 0. c=======================================================================******** c.... Type 8. Hyperbolic paraboloid: c.... (+/-)sz * z + sxx * x**2 - |syy| * y**2 = 0. 108 aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'HYPERBOLIC PARABOLOID.')") cbugcbugc***DEBUG ends. if (aar .eq. 0.0) then ! Point "a" is on the z axis. if (sz * aaz .lt. -0.5 * sz**2 / sxx) then nprox = 2 ! Exact. cy(1) = 0.0 cz(1) = aaz + 0.5 * sz / sxx cx(1) = sqrt (-sz * cz(1) / sxx) dist(1) = sqrt (cx(1)**2 + (cz(1) - aaz)**2) cy(2) = cy(1) cz(2) = cz(1) cx(2) = -cx(1) dist(2) = dist(1) elseif (sz * aaz .gt. -0.5 * sz**2 / syy) then nprox = 2 ! Exact. cx(1) = 0.0 cz(1) = aaz + 0.5 * sz / syy cy(1) = sqrt (-sz * cz(1) / syy) dist(1) = sqrt (cy(1)**2 + (cz(1) - aaz)**2) cy(2) = -cy(1) cz(2) = cz(1) cx(2) = cx(1) dist(2) = dist(1) else nprox = 1 ! Exact. At the origin. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = aaz endif ! Tested location of point "a" on z axis. elseif (aax .eq. 0.0) then ! Point "a" is in plane x = 0. cyy = aay / (1.0 - syy / sxx) czz = aaz + 0.5 * sz / sxx arg = -(sz * czz + syy * cyy**2) / sxx if (arg .ge. 0.0) then ! Proximal points not in plane x = 0. nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = czz dist(1) = sqrt (cx(1)**2 + (cy(1) - aay)**2 + & (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Proximal point in plane x = 0. c0 = sz * aaz + syy * aay**2 c1 = sz * (sz - 4.0 * syy * aaz) c2 = 4.0 * syy * sz * (syy * aaz - sz) c3 = 4.0 * syy**2 * sz**2 c0 = c0 / c3 c1 = c1 / c3 c2 = c2 / c3 call aptcubs (c0, c1, c2, tol, nroot, rreal, rimag, & atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. At x = 0. distsqm = 1.e99 do 260 n = 1, nroot cyy = aay / (1.0 - 2.0 * syy * rreal(n) + fuz) czz = aaz + sz * rreal(n) distsq = (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = 0.0 cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 260 continue endif ! Tested sign of arg. elseif (aay .eq. 0.0) then ! Point "a" is in plane y = 0. cxx = aax / (1.0 - sxx / syy) czz = aaz + 0.5 * sz / syy arg = -(sz * czz + sxx * cxx**2) / syy if (arg .ge. 0.0) then ! Proximal points not in plane y = 0. nprox = 2 cx(1) = cxx cy(1) = sqrt (arg) cz(1) = czz dist(1) = sqrt (cy(1)**2 + (cx(1) - aax)**2 + & (cz(1) - aaz)**2) cx(2) = cx(1) cy(2) = -cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Proximal points in plane y = 0. c0 = sz * aaz + sxx * aax**2 c1 = sz * (sz - 4.0 * sxx * aaz) c2 = 4.0 * sxx * sz * (sxx * aaz - sz) c3 = 4.0 * sxx**2 * sz**2 c0 = c0 / c3 c1 = c1 / c3 c2 = c2 / c3 call aptcubs (c0, c1, c2, tol, nroot, rreal, rimag, & atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. At y = 0. distsqm = 1.e99 do 270 n = 1, nroot cxx = aax / (1.0 - 2.0 * sxx * rreal(n) + fuz) czz = aaz + sz * rreal(n) distsq = (cxx - aax)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cy(1) = 0.0 cx(1) = cxx cz(1) = czz dist(1) = sqrt (distsq) endif 270 continue endif ! Tested sign of arg. else ! General 5th-order elliptic paraboloid. c.... Find coeffients of 5th-order equation for hyperbolic paraboloid. c.... f = (cx - aax) / (2.0 * sxx * cx) (aax not zero) c.... f = (cy - aay) / (2.0 * syy * cy) (aay not zero) c.... f = (cz - aaz) / sz (sz < 0) fsum1 = sxx + syy fsum2 = sxx * syy fsum3 = fsum1**2 + 2.0 * fsum2 sza = abs (sz) abz = abs (aaz) a(0) = sxx * aax**2 + syy * aay**2 + sz * aaz a(1) = sz**2 - 4.0 * (fsum1 * sz * aaz + & fsum2 * (aax**2 + aay**2)) a(2) = 4.0 * (fsum3 * sz * aaz - fsum1 * sz**2 + & fsum2 * (syy * aax**2 + sxx * aay**2)) a(3) = 4.0 * sz * (fsum3 * sz - 4.0 * fsum1 * fsum2 * aaz) a(4) = 16.0 * fsum2 * sz * (fsum2 * aaz - fsum1 * sz) a(5) = 16.0 * fsum2**2 * sz**2 cbugcbugc***DEBUG begins. cbugcbug 9928 format (/ 'a(n) =',1p6e12.4 ) cbugcbug write ( 3, 9928) (a(nn), nn = 0, 5) cbugcbugc***DEBUG ends. c.... Find limits on f for an hyperbolic paraboloid. c.... Find out if external point is inside or outside. funq = sxx * aax**2 + syy * aay**2 side = funq + sz * aaz serr = tol * (3.0 * abs (sxx) * aax**2 + & 3.0 * abs (syy) * aay**2 + & 2.0 * abs (sz * aaz)) cbugcbugc***DEBUG begins. cbugcbug 9931 format (/ 'side, err =',1p2e20.12) cbugcbug write ( 3, 9931) side, serr cbugcbugc***DEBUG ends. argx = -(syy * aay**2 + sz * aaz) / sxx argy = -(sz * aaz + sxx * aax**2) / syy argz = -funq / sz zint = argz if (side .le. 0.0) then ! Inside. xint = sqrt (argx) fmaxx = 0.5 * (1.0 - abs (aax) / xint) / sxx if (argy .gt. 0.0) then yint = sqrt (argy) fmaxy = 0.5 * (1.0 - abs (aay) / yint) / syy else fmaxy = 1.e99 endif fmaxz = (zint - aaz) / sz fmin = 0.0 - 2.0 * tol fmax = min (fmaxx, fmaxy, fmaxz) + 2.0 * tol finit = fmax else ! Outside. if (argx .gt. 0.0) then xint = sqrt (argx) fminx = 0.5 * (1.0 - abs (aax) / xint) / sxx else fminx = -1.e99 endif yint = sqrt (argy) fminy = 0.5 * (1.0 - abs (aay) / yint) / syy fminz = (zint - aaz) / sz fmin = max (fminx, fminy, fminz) - 2.0 * tol fmax = 2.0 * tol finit = fmin endif ! Tested inside or outside. cbugcbugc***DEBUG begins. cbugcbug 9936 format (/ 'fmin, fmax =',1p2e20.12) cbugcbug write ( 3, 9936) fmin, fmax cbugcbugc***DEBUG ends. c.... Find froot, value of f at proximal point. call aptnewt (a, 5, finit, fmin, fmax, 1000, tol, & froot, nerra) if (nerra .eq. 0) then nprox = 1 denx = 1.0 - 2.0 * sxx * froot + fuz deny = 1.0 - 2.0 * syy * froot + fuz cx(1) = aax / (denx + fuz) cy(1) = aay / (deny + fuz) cz(1) = aaz + sz * froot cbugcbugc***DEBUG begins. cbugcbug 9972 format (' cx, cy, cz =',1p3e22.14) cbugcbug write ( 3, 9972) cx(1), cy(1), cz(1) cbugcbugc***DEBUG ends. c.... Make a correction for an initial coordinate very close to an axis. errx = 1.0 + 4.0 * abs (sxx * froot) erry = 1.0 + 4.0 * abs (syy * froot) ratx = abs (denx) / errx raty = abs (deny) / erry ratr = min (ratx, raty) if (ratx .eq. ratr) then argx = -(syy * cy(1)**2 + sz * cz(1)) / sxx cx(1) = sqrt (max (0.0, argx)) if (aax .lt. 0.0) cx(1) = -cx(1) endif if (raty .eq. ratr) then argy = -(sxx * cx(1)**2 + sz * cz(1)) / syy cy(1) = sqrt (max (0.0, argy)) if (aay .lt. 0.0) cy(1) = -cy(1) endif cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9972) cx(1), cy(1), cz(1) cbugcbugc***DEBUG ends. dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) else ! General solution failed. if (nprox .eq. 0) go to 710 ! Found normal intercept. c.... Use an approximate proximal point at the origin. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = sqrt (aax**2 + aay**2 + aaz**2) endif ! Tested for general solution. endif ! Tested for solution type. go to 510 c.... Type 8. Hyperbolic paraboloid: c.... (+/-)sz * z + sxx * x**2 + |syy| * y**2 = 0. c=======================================================================******** c.... Type 9. Elliptic paraboloid: c.... (+/-)sz * z + sxx * x**2 + |syy| * y**2 = 0. 109 aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'ELLIPTIC PARABOLOID.')") cbugcbugc***DEBUG ends. if (aar .eq. 0.0) then ! Point "a" is on the z axis. if (sz * aaz .ge. -0.5 * sz**2 / sxx) then nprox = 1 ! Exact. At the origin. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = aaz else nprox = 2 ! Exact. cy(1) = 0.0 cz(1) = aaz + 0.5 * sz / sxx cx(1) = sqrt (-sz * cz(1) / sxx) dist(1) = sqrt (cx(1)**2 + (cz(1) - aaz)**2) cy(2) = cy(1) cz(2) = cz(1) cx(2) = -cx(1) dist(2) = dist(1) endif ! Tested location of point "a" on z axis. elseif (aax .eq. 0.0) then ! Point "a" is in plane x = 0. cyy = aay / (1.0 - syy / sxx) czz = aaz + 0.5 * sz / sxx arg = -(sz * czz + syy * cyy**2) / sxx if (arg .ge. 0.0) then ! Proximal points not in plane x = 0. nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = czz dist(1) = sqrt (cx(1)**2 + (cy(1) - aay)**2 + & (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Proximal point in plane x = 0. c0 = sz * aaz + syy * aay**2 c1 = sz * (sz - 4.0 * syy * aaz) c2 = 4.0 * syy * sz * (syy * aaz - sz) c3 = 4.0 * syy**2 * sz**2 c0 = c0 / c3 c1 = c1 / c3 c2 = c2 / c3 call aptcubs (c0, c1, c2, tol, nroot, rreal, rimag, & atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. At x = 0. distsqm = 1.e99 do 280 n = 1, nroot cyy = aay / (1.0 - 2.0 * syy * rreal(n) + fuz) czz = aaz + sz * rreal(n) distsq = (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = 0.0 cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 280 continue endif ! Tested sign of arg. elseif (aay .eq. 0.0) then ! Point "a" is in plane y = 0. cxx = aax / (1.0 - sxx / syy) czz = aaz + 0.5 * sz / syy arg = -(sz * czz + sxx * cxx**2) / syy if (arg .ge. 0.0) then ! Proximal points not in plane y = 0. nprox = 2 cx(1) = cxx cy(1) = sqrt (arg) cz(1) = czz dist(1) = sqrt (cy(1)**2 + (cx(1) - aax)**2 + & (cz(1) - aaz)**2) cx(2) = cx(1) cy(2) = -cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Proximal points in plane y = 0. c0 = sz * aaz + sxx * aax**2 c1 = sz * (sz - 4.0 * sxx * aaz) c2 = 4.0 * sxx * sz * (sxx * aaz - sz) c3 = 4.0 * sxx**2 * sz**2 c0 = c0 / c3 c1 = c1 / c3 c2 = c2 / c3 call aptcubs (c0, c1, c2, tol, nroot, rreal, rimag, & atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. At y = 0. distsqm = 1.e99 do 290 n = 1, nroot cxx = aax / (1.0 - 2.0 * sxx * rreal(n) + fuz) czz = aaz + sz * rreal(n) distsq = (cxx - aax)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cy(1) = 0.0 cx(1) = cxx cz(1) = czz dist(1) = sqrt (distsq) endif 290 continue endif ! Tested sign of arg. else ! General 5th-order elliptic paraboloid. c.... Find coeffients of 5th-order equation for elliptic paraboloid. c.... f = (cx - aax) / (2.0 * sxx * cx) (aax not zero) c.... f = (cy - aay) / (2.0 * syy * cy) (aay not zero) c.... f = (cz - aaz) / sz (sz < 0) fsum1 = sxx + syy fsum2 = sxx * syy fsum3 = fsum1**2 + 2.0 * fsum2 sza = abs (sz) abz = abs (aaz) a(0) = sxx * aax**2 + syy * aay**2 + sz * aaz a(1) = sz**2 - 4.0 * (fsum1 * sz * aaz + & fsum2 * (aax**2 + aay**2)) a(2) = 4.0 * (fsum3 * sz * aaz - fsum1 * sz**2 + & fsum2 * (syy * aax**2 + sxx * aay**2)) a(3) = 4.0 * sz * (fsum3 * sz - 4.0 * fsum1 * fsum2 * aaz) a(4) = 16.0 * fsum2 * sz * (fsum2 * aaz - fsum1 * sz) a(5) = 16.0 * fsum2**2 * sz**2 cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9928) (a(nn), nn = 0, 5) cbugcbugc***DEBUG ends. c.... Find limits on f for an elliptic paraboloid. c.... Find out if external point is inside or outside. funq = sxx * aax**2 + syy * aay**2 side = funq + sz * aaz serr = tol * (3.0 * abs (sxx) * aax**2 + & 3.0 * abs (syy) * aay**2 + & 2.0 * abs (sz * aaz)) cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9931) side, serr cbugcbugc***DEBUG ends. if (side .le. 0.0) then ! Inside. argx = -(syy * aay**2 + sz * aaz) / sxx if (argx .gt. 0.0) then xint = sqrt (argx) fmaxx = 0.5 * (1.0 - abs (aax) / xint) / sxx else fmaxx = 1.e99 endif argy = -(sz * aaz + sxx * aax**2) / syy if (argy .gt. 0.0) then yint = sqrt (argy) fmaxy = 0.5 * (1.0 - abs (aay) / yint) / syy else fmaxy = 1.e99 endif argz = -funq / sz zint = argz fmaxz = (zint - aaz) / sz if (funq .gt. 0.0) then rat = sqrt (-sz * aaz / funq) fmaxa = 0.5 * (1.0 - 1.0 / rat) / syy else fmaxa = 0.5 / max (sxx, syy) endif fmin = 0.0 - 2.0 * tolx fmax = min (fmaxx, fmaxy, fmaxz, fmaxa) + 2.0 * tolx finit = fmax else ! Outside. argx = -(syy * aay**2 + sz * aaz) / sxx if (argx .gt. 0.0) then xint = sqrt (argx) fminx = 0.5 * (1.0 - abs (aax) / xint) / sxx else fminx = -1.e99 endif argy = -(sz * aaz + sxx * aax**2) / syy if (argy .gt. 0.0) then yint = sqrt (argy) fminy = 0.5 * (1.0 - abs (aay) / yint) / syy else fminy = -1.e99 endif argz = -funq / sz zint = argz fminz = (zint - aaz) / sz zfoc = -0.5 * sz / sxx ! Use lower focal point. if (aaz / sz .le. zfoc / sz) then if (funq .gt. 0.0) then rat = sqrt (-sz * aaz / funq) fmina = 0.5 * (1.0 - 1.0 / rat) / syy endif else rat = 1.e99 fmina = -1.e99 endif fmin = max (fminx, fminy, fminz, fmina) - 2.0 * tolx fmax = min (-aaz / sz, 0.0) + 2.0 * tolx finit = fmin endif ! Tested inside or outside. cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9936) fmin, fmax cbugcbugc***DEBUG ends. c.... Find froot, value of f at proximal point. call aptnewt (a, 5, finit, fmin, fmax, 1000, tol, & froot, nerra) if (nerra .eq. 0) then nprox = 1 denx = 1.0 - 2.0 * sxx * froot + fuz deny = 1.0 - 2.0 * syy * froot + fuz cx(1) = aax / (denx + fuz) cy(1) = aay / (deny + fuz) cz(1) = aaz + sz * froot cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9972) cx(1), cy(1), cz(1) cbugcbugc***DEBUG ends. c.... Make a correction for an initial coordinate very close to an axis. errx = 1.0 + 4.0 * abs (sxx * froot) erry = 1.0 + 4.0 * abs (syy * froot) ratx = abs (denx) / errx raty = abs (deny) / erry ratr = min (ratx, raty) if (ratx .eq. ratr) then argx = -(syy * cy(1)**2 + sz * cz(1)) / sxx cx(1) = sqrt (max (0.0, argx)) if (aax .lt. 0.0) cx(1) = -cx(1) endif if (raty .eq. ratr) then argy = -(sxx * cx(1)**2 + sz * cz(1)) / syy cy(1) = sqrt (max (0.0, argy)) if (aay .lt. 0.0) cy(1) = -cy(1) endif cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9972) cx(1), cy(1), cz(1) cbugcbugc***DEBUG ends. dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) else ! General solution failed. if (nprox .eq. 0) go to 710 ! Found normal intercept. arg = -sz * aaz if (arg .le. 0.0) then ! Use the vertex point at the origin. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = sqrt (aax**2 + aay**2 + aaz**2) else ! Use a point on the ellipse at z = aaz. cx(1) = aax * sqrt (arg / sxx) / aar cy(1) = aay * sqrt (arg / syy) / aar cz(1) = aaz dist(1) = sqrt((cx(1) - aax)**2 + (cy(1) - aay)**2) endif ! Tested z axis intercept. endif ! Tested for general solution. endif ! Tested for solution type. go to 510 c.... Type 9. Elliptic paraboloid: c.... (+/-)sz * z + sxx * x**2 + |syy| * y**2 = 0. c=======================================================================******** c.... Type 10. Circular paraboloid: (+/-)sz * z + sxx * (x**2 + y**2) = 0. 110 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'CIRCULAR PARABOLOID.')") cbugcbugc***DEBUG ends. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if (aar .eq. 0.0) then ! Point "a" is on the z axis. arg = -sz * (aaz + 0.5 * sz / sxx) / sxx if (arg .le. 0.0) then ! Proximal point is at vertex. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = aaz else ! Proximal points are on a circle. nprox = 3 ! Exact. 3 points on proximal circle. rc = sqrt (arg) ! Radius of proximal circle. cx(1) = rc cy(1) = 0.0 cz(1) = aaz + 0.5 * sz / sxx cx(2) = -0.5 * rc cy(2) = 0.5 * sqrt (3.0) * rc cz(2) = cz(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(1) = sqrt (arg + (0.5 * sz / sxx)**2) dist(2) = dist(1) dist(3) = dist(1) endif else ! Point "a" is not on the z axis. c.... Solve cubic for one to three values of the radius of the proximal point. aacub = 0.0 bbcub = sz * (aaz + 0.5 * sz / sxx) / sxx cccub = -0.5 * sz**2 * aar / sxx**2 call aptcubs (cccub, bbcub, aacub, tol, nroots, rreal, rimag, & atype, itrun) if (nroots .ge. 1) then crtest(1) = rreal(1) cztest(1) = -sxx * crtest(1)**2 / sz distest(1) = sqrt ((crtest(1) - aar)**2 + & (cztest(1) - aaz)**2) nuse = 1 endif if (nroots .ge. 2) then crtest(2) = rreal(2) cztest(2) = -sxx * crtest(2)**2 / sz distest(2) = sqrt ((crtest(2) - aar)**2 + & (cztest(2) - aaz)**2) if (distest(2) .lt. distest(1)) nuse = 2 endif if (nroots .eq. 3) then crtest(3) = rreal(3) cztest(3) = -sxx * crtest(3)**2 / sz distest(3) = sqrt ((crtest(3) - aar)**2 + & (cztest(3) - aaz)**2) if (distest(3) .lt. distest(nuse)) nuse = 3 endif cx(1) = aax * crtest(nuse) / aar cy(1) = aay * crtest(nuse) / aar cz(1) = cztest(nuse) dist(1) = distest(nuse) endif ! Tested for point "a" on the z axis. go to 510 c.... Type 10. Circular paraboloid: (+/-)sz * z + sxx * (x**2 + y**2) = 0. c=======================================================================******** c.... Type 11. Real elliptic cone: c.... sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. 111 aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'ELLIPTIC CONE.')") cbugcbugc***DEBUG ends. c.... Test for point "a" at origin. if ((aar .eq. 0.0) .and. (aaz .eq. 0.0)) then ! At origin. nprox = 1 ! Exact. At point "a" at origin. cx(1) = aax cy(1) = aay cz(1) = aaz dist(1) = 0.0 c.... Test for point "a" on x axis. elseif (aar .eq. 0.0) then ! On z axis. nprox = 2 ! Exact. cz(1) = aaz / (1.0 - szz / sxx) cx(1) = sqrt (-szz * cz(1)**2 / sxx) cy(1) = 0.0 cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(1) = sqrt (cx(1)**2 + (cz(1) - aaz)**2) dist(2) = dist(1) c.... Test for point "a" in plane x = 0. elseif ((aax .eq. 0.0) .and. & (aay .ne. 0.0) .and. & (aaz .ne. 0.0)) then yp = aay / (1.0 - syy / sxx) zp = aaz / (1.0 - szz / sxx) arg = -(syy * yp**2 + szz * zp**2) / sxx if (arg .lt. 0.0) then ! Near extreme y. c0 = syy * aay**2 + szz * aaz**2 c1 = -4.0 * syy * szz * (aay**2 + aaz**2) c2 = 4.0 * syy * szz * (szz * aay**2 + syy * aaz**2) call aptqrts (0, c2, c1, c0, qq, tol, nroots, root1, root2, & itrun) if (nroots .eq. 0) then nerr = 3 go to 710 endif nprox = 1 ! Exact. At x = 0, at extreme y. if (nroots .eq. 1) root2 = root1 xp1 = 0.0 yp1 = aay / (1.0 - 2.0 * syy * root1 + fuz) zp1 = aaz / (1.0 - 2.0 * szz * root1 + fuz) dp1 = sqrt ((yp1 - aay)**2 + (zp1 - aaz)**2) xp2 = 0.0 yp2 = aay / (1.0 - 2.0 * syy * root2 + fuz) zp2 = aaz / (1.0 - 2.0 * szz * root2 + fuz) dp2 = sqrt ((yp2 - aay)**2 + (zp2 - aaz)**2) if (dp1 .le. dp2) then cx(1) = xp1 cy(1) = yp1 cz(1) = zp1 dist(1) = dp1 else cx(1) = xp2 cy(1) = yp2 cz(1) = zp2 dist(1) = dp2 endif else ! Not near extreme y. nprox = 2 ! Exact. Not at x = 0. cx(1) = sqrt (arg) cy(1) = yp cz(1) = zp dist(1) = sqrt (cx(1)**2 + (yp - aay)**2 + (zp - aaz)**2) cx(2) = -sqrt (arg) cy(2) = yp cz(2) = zp dist(2) = sqrt (cx(1)**2 + (yp - aay)**2 + (zp - aaz)**2) endif ! Tested point "a" nearness to ylim. c.... Test for point "a" in plane y = 0. elseif ((aay .eq. 0.0) .and. & (aaz .ne. 0.0) .and. & (aax .ne. 0.0)) then c0 = sxx * aax**2 + szz * aaz**2 c1 = -4.0 * sxx * szz * (aax**2 + aaz**2) c2 = 4.0 * sxx * szz * (szz * aax**2 + sxx * aaz**2) call aptqrts (0, c2, c1, c0, qq, tol, nroots, root1, root2, & itrun) if (nroots .eq. 0) then nerr = 3 go to 710 endif nprox = 1 ! Exact. At y = 0. if (nroots .eq. 1) root2 = root1 xp1 = aax / (1.0 - 2.0 * sxx * root1 + fuz) yp1 = 0.0 zp1 = aaz / (1.0 - 2.0 * szz * root1 + fuz) dp1 = sqrt ((xp1 - aax)**2 + (zp1 - aaz)**2) xp2 = aax / (1.0 - 2.0 * sxx * root2 + fuz) yp2 = 0.0 zp2 = aaz / (1.0 - 2.0 * szz * root2 + fuz) dp2 = sqrt ((xp2 - aax)**2 + (zp2 - aaz)**2) if (dp1 .le. dp2) then cx(1) = xp1 cy(1) = yp1 cz(1) = zp1 dist(1) = dp1 else cx(1) = xp2 cy(1) = yp2 cz(1) = zp2 dist(1) = dp2 endif c.... Test for point "a" in plane z = 0. elseif ((aaz .eq. 0.0) .and. & (aax .ne. 0.0) .and. & (aay .ne. 0.0)) then nprox = 2 ! Exact. cx(1) = aax / (1.0 - qxx / qzz) cy(1) = aay / (1.0 - qyy / qzz) cz(1) = sqrt (-(sxx * cx(1)**2 + syy * cy(1)**2) / szz) dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + cz(1)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) c.... Use general quartic solution. else ! Point "a" not on main axes or planes. c0 = sxx * aax**2 + syy * aay**2 + szz * aaz**2 c1 = -4.0 * (sxx * (syy + szz) * aax**2 + & syy * (szz + sxx) * aay**2 + & szz * (sxx + syy) * aaz**2 ) c2 = 4.0 * & (sxx * (syy**2 + 4.0 * syy * szz + szz**2) * aax**2 + & syy * (szz**2 + 4.0 * szz * sxx + sxx**2) * aay**2 + & szz * (sxx**2 + 4.0 * sxx * syy + syy**2) * aaz**2 ) c3 = -16.0 * sxx * syy * szz * & ((syy + szz) * aax**2 + & (szz + sxx) * aay**2 + & (sxx + syy) * aaz**2 ) c4 = 16.0 * sxx * syy * szz * & (syy * szz * aax**2 + & szz * sxx * aay**2 + & sxx * syy * aaz**2) c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 250 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cxx - aax)**2 + (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 250 continue endif ! Tested for position of point "a". go to 510 c.... Type 11. Real elliptic cone: c.... sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. c=======================================================================******** c.... Type 12. Real circular cone: sxx * (x**2 + y**2) - |szz| * z**2 = 0. 112 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'CIRCULAR CONE.')") cbugcbugc***DEBUG ends. crat = -szz / sxx cratp = 1.0 + crat cratr = sqrt (crat) aar = sqrt (aax**2 + aay**2) rp = cratr * (abs (aaz) + cratr * aar) / cratp if (aar .eq. 0.0) then ! Point "a" is on the axis of the cone. if (aaz .eq. 0.0) then ! Point "a" is at the vertex of the cone. nprox = 1 ! Exact. One point at the vertex. cx(1) = 0.0 cy(1) = 0.0 cz(1) = 0.0 dist(1) = 0.0 else ! Point "a" is not at the vertex. nprox = 3 ! Exact. 3 points on the proximal ring. cx(1) = rp cy(1) = 0.0 cz(1) = aaz / cratp cx(2) = -0.5 * rp cy(2) = 0.5 * sqrt (3.0) * rp cz(2) = cz(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(1) = cratr * abs (aaz) / sqrt (cratp) dist(2) = dist(1) dist(3) = dist(1) endif else ! Point "a" is not on the axis of the cone. if (aaz .eq. 0.0) then ! Point "a" is in the vertex plane. nprox = 2 ! Exact. cx(1) = aax * rp / aar cy(1) = aay * rp / aar cz(1) = cratr * aar / cratp dist(1) = aar / sqrt (cratp) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) else ! Point "a" is not in the vertex plane. nprox = 1 ! Exact. cx(1) = aax * rp / aar cy(1) = aay * rp / aar cz(1) = sign ((abs (aaz) + cratr * aar) / cratp, aaz) dist(1) = abs (cratr * abs (aaz) - aar) / sqrt (cratp) endif endif ! Tested for point on the cone axis. go to 510 c.... Type 12. Real circular cone: sxx * (x**2 + y**2) - |szz| * z**2 = 0. c=======================================================================******** c.... Type 13. Hyperboloid of one sheet: c.... -|sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. 113 if (abs (sxx - syy) .le. tol * (sxx + syy)) then ! Circular. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'CIRCULAR HYPERBOLOID OF ONE SHEET.')") cbugcbugc***DEBUG ends. c=======================================================================******** c.... Type 13a. Circular hyperboloid of one sheet: c.... -|sc| + sxx * (x**2 + y**2) - |szz| * z**2 = 0. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if ((aar .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Proximal ring at z = 0. nprox = 3 ! Exact. 3 points on the proximal circle. rc = sqrt (-sc / sxx) ! Radius of proximal circle. cx(1) = rc cy(1) = 0.0 cz(1) = 0.0 dist(1) = rc cx(2) = -0.5 * rc cy(2) = 0.5 * sqrt (3.0) * rc cz(2) = cz(1) dist(2) = dist(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(3) = dist(1) elseif (aar .eq. 0.0) then ! Proximal ring, not at z = 0. nprox = 3 ! Exact. 3 points on the proximal circle. cz(1) = aaz / (1.0 - szz / sxx) rc = sqrt (-(sc + szz * cz(1)**2) / sxx) ! Radius of prox circle. cx(1) = rc cy(1) = 0.0 dist(1) = sqrt (rc**2 + (cz(1) - aaz)**2) cx(2) = -0.5 * rc cy(2) = 0.5 * sqrt (3.0) * rc cz(2) = cz(1) dist(2) = dist(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(3) = dist(1) elseif (aaz .eq. 0.0) then rc = aar / (1.0 - sxx / szz) ! Radius of proximal circle. arg = -(sc + sxx * rc**2) / szz if (arg .le. 0.0) then ! Proximal point on circle at z = 0. nprox = 1 ! Exact. rc = sqrt (-sc / sxx) ! Radius of waist. cx(1) = aax * rc / aar cy(1) = aay * rc / aar cz(1) = 0.0 dist(1) = rc - aar else ! Two proximal points, not at z = 0. nprox = 2 ! Exact. cx(1) = aax * rc / aar cy(1) = aay * rc / aar cz(1) = sqrt (arg) dist(1) = sqrt ((rc - aar)**2 + arg) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) endif ! Tested arg. else ! Solve general quartic equation. c.... Find the roots of the equation for |d| / |N|. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. c0 = sc + sxx * aar**2 + szz * aaz**2 c1 = -4.0 * (sc * (sxx + szz) + sxx * szz * (aar**2 + aaz**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aar**2 + sxx * aaz**2)) c3 = -16.0 * sc * sxx * szz * (sxx + szz) c4 = 16.0 * sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 310 n = 1, nroot denr = abs (1.0 - 2.0 * sxx * rreal(n) + fuz) cxx = aax / denr cyy = aay / denr czz = aaz / (1.0 - 2.0 * szz * rreal(n) + fuz) distsq = (cxx - aax)**2 + (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 310 continue endif ! Tested for symmetry condition. c.... Type 13a. Circular hyperboloid of one sheet: c.... -|sc| + sxx * (x**2 + y**2) - |szz| * z**2 = 0. c=======================================================================******** c.... Middle of logic for hyperboloid of one sheet. else ! Elliptic hyperboloid of one sheet. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'ELLIPTIC HYPERBOLOID OF ONE SHEET.')") cbugcbugc***DEBUG ends. c=======================================================================******** c.... Type 13b. Elliptic hyperboloid of one sheet: c.... -|sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. cbugcbugc***DEBUG begins. cbugcbug 9981 format ('DEBUG. Doing type 13b.') cbugcbug write ( 3, 9981) cbugcbugc***DEBUG ends. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if ((aar .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" at origin. nprox = 2 ! Exact, at waist, on x semiaxis. cx(1) = sqrt (-sc / sxx) cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) elseif ((aax .eq. 0.0) .and. (aay .eq. 0.0)) then ! Point "a" on z axis. cbugcbugc***DEBUG begins. cbugcbug 9982 format ('DEBUG. Point "a" on z axis.') cbugcbug write ( 3, 9982) cbugcbugc***DEBUG ends. nprox = 2 ! Exact. cz(1) = aaz / (1.0 - szz / sxx) cy(1) = 0.0 cx(1) = sqrt (-(sc + szz * cz(1)**2) / sxx) dist(1) = sqrt (cx(1)**2 + (cz(1) - aaz)**2) cz(2) = cz(1) cy(2) = cy(1) cx(2) = -cx(1) dist(2) = dist(1) elseif ((aay .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" on x axis. cxx = aax / (1.0 - sxx / szz) arg = -(sc + sxx * cxx**2) / szz if (arg .le. 0.0) then nprox = 1 ! Exact. cx(1) = (aax / abs (aax)) * sqrt (-sc / sxx) cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) - aax else nprox = 2 ! Exact cx(1) = cxx cy(1) = 0.0 cz(1) = sqrt (arg) dist(1) = sqrt (arg + (cx(1) - aax)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) endif ! Tested sign of arg. elseif ((aax .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" on y axis. cyyx = aay / (1.0 - syy / sxx) argx = -(sc + syy * cyyx**2) / sxx cyyz = aay / (1.0 - syy / szz) argz = -(sc + syy * cyyz**2) / szz if (argx .gt. 0.0) then nprox = 2 ! Exact, at non-zero cx. cx(1) = sqrt (argx) cy(1) = cyyx cz(1) = 0.0 dist(1) = sqrt (argx + (cy(1) - aay)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) elseif (argz .le. 0.0) then ! Proximal point at y vertex. nprox = 1 ! Exact, at y vertex. cx(1) = 0.0 cy(1) = (aay / abs (aay)) * sqrt (-sc / syy) cz(1) = 0.0 dist(1) = cy(1) - aay else ! Sign of argz is positive. nprox = 2 ! Exact, at non-zero cz. cx(1) = 0.0 cy(1) = cyyz cz(1) = sqrt (argz) dist(1) = sqrt (argz + (cy(1) - aay)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) endif ! Tested signs of argx and argz. elseif (aax .eq. 0.0) then ! In plane x = 0. cyy = aay / (1.0 - syy / sxx) czz = aaz / (1.0 - szz / sxx) arg = -(sc + syy * cyy**2 + szz * czz**2) / sxx if (arg .ge. 0.0) then nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = czz dist(1) = sqrt (arg + (cy(1) - aay)**2 + (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Sign of arg is negative. c0 = sc + syy * aay**2 + szz * aaz**2 c1 = -4.0 * (sc * (syy + szz) + syy * szz * (aay**2 + aaz**2)) c2 = 4.0 * (sc * (syy**2 + 4.0 * syy * szz + szz**2) + & syy * szz * (szz * aay**2 + syy * aaz**2)) c3 = -16.0 * sc * syy * szz * (syy + szz) c4 = 16 * sc * syy**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 330 n = 1, nroot deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = 0.0 cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 330 continue endif ! Tested sign of arg. elseif (aay .eq. 0.0) then ! In plane y = 0. c0 = sc + sxx * aax**2 + szz * aaz**2 c1 = -4.0 * (sc * (sxx + szz) + sxx * szz * (aax**2 + aaz**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aax**2 + sxx * aaz**2)) c3 = -16.0 * sc * sxx * szz * (sxx + szz) c4 = 16 * sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 340 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cxx - aax)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = 0.0 cz(1) = czz dist(1) = sqrt (distsq) endif 340 continue elseif (aaz .eq. 0.0) then ! In plane z = 0. c0 = sc + sxx * aax**2 + syy * aay**2 c1 = -4.0 * (sc * (sxx + syy) + sxx * syy * (aax**2 + aay**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * syy + syy**2) + & sxx * syy * (syy * aax**2 + sxx * aay**2)) c3 = -16.0 * sc * sxx * syy * (sxx + syy) c4 = 16 * sc * sxx**2 * syy**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 350 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny distsq = (cxx - aax)**2 + (cyy - aay)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = 0.0 dist(1) = sqrt (distsq) endif 350 continue else ! General 6th order equation. c.... General elliptic hyperboloid of one sheet. c.... Find the coeffients of the 6th-order equation for f. c.... f = (cx - aax) / (2.0 * sxx * cx) (aax not zero) c.... f = (cy - aay) / (2.0 * syy * cy) (aay not zero) c.... f = (cz - aaz) / (2.0 * szz * cz) (aaz not zero) fsum1 = sxx + syy + szz fsum2 = sxx * syy + syy * szz + szz * sxx fsum3 = sxx * syy * szz a(0) = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 a(1) = -4.0 * (sc * fsum1 + & sxx * (syy + szz) * aax**2 + & syy * (szz + sxx) * aay**2 + & szz * (sxx + syy) * aaz**2 ) a(2) = 4.0 * (sc * (fsum1**2 + 2.0 * fsum2) + & sxx * (syy**2 + 4.0 * syy * szz + szz**2) * aax**2 + & syy * (szz**2 + 4.0 * szz * sxx + sxx**2) * aay**2 + & szz * (sxx**2 + 4.0 * sxx * syy + syy**2) * aaz**2 ) a(3) = -16.0 * (sc * (fsum1 * fsum2 + fsum3) + & fsum3 * ((syy + szz) * aax**2 + & (szz + sxx) * aay**2 + & (sxx + syy) * aaz**2 )) a(4) = 16.0 * (sc * (fsum2**2 + 2.0 * fsum1 * fsum3) + & fsum3 * syy * szz * aax**2 + & fsum3 * szz * sxx * aay**2 + & fsum3 * sxx * syy * aaz**2 ) a(5) = -64.0 * sc * fsum2 * fsum3 a(6) = 64.0 * sc * fsum3**2 c.... Find limits on f for an elliptic hyperboloid of one sheet. c.... Find out if external point is inside or outside. side = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 if (side .le. 0.0) then ! Inside. xint = sqrt (-(sc + syy * aay**2 + szz * aaz**2) / sxx) fmaxx = 0.5 * (1.0 - abs (aax) / xint) / sxx yint = sqrt (-(sc + szz * aaz**2 + sxx * aax**2) / syy) fmaxy = 0.5 * (1.0 - abs (aay) / yint) / syy argz = -(sc + sxx * aax**2 + syy * aay**2) / szz if (argz .gt. 0.0) then zint = sqrt (argz) fmaxz = 0.5 * (1.0 - abs (aaz) / zint) / szz else fmaxz = 1.e99 endif fmin = 0.0 - 2.0 * tolx fmax = min (fmaxx, fmaxy, fmaxz) + 2.0 * tolx finit = fmax else ! Outside. argx = -(sc + syy * aay**2 + szz * aaz**2) / sxx if (argx .gt. 0.0) then xint = sqrt (argx) fminx = 0.5 * (1.0 - abs (aax) / xint) / sxx else fminx = -1.e99 endif argy = -(sc + szz * aaz**2 + sxx * aax**2) / syy if (argy .gt. 0.0) then yint = sqrt (argy) fminy = 0.5 * (1.0 - abs (aay) / yint) / syy else fminy = -1.e99 endif argz = -(sc + sxx * aax**2 + syy * aay**2) / szz zint = sqrt (argz) fminz = 0.5 * (1.0 - abs (aaz) / zint) / szz fmin = max (fminx, fminy, fminz) - 2.0 * tolx fmax = 0.0 + 2.0 * tolx finit = fmin endif ! Tested inside or outside. c.... Find froot, value of f at proximal point. call aptnewt (a, 6, finit, fmin, fmax, 1000, tol, & froot, nerra) if (nerra .eq. 0) then nprox = 1 denx = 1.0 - 2.0 * sxx * froot + fuz deny = 1.0 - 2.0 * syy * froot + fuz denz = 1.0 - 2.0 * szz * froot + fuz cx(1) = aax / (denx + fuz) cy(1) = aay / (deny + fuz) cz(1) = aaz / (denz + fuz) c.... Make a correction for an initial coordinate very close to an axis. errx = 1.0 + 4.0 * abs (sxx * froot) erry = 1.0 + 4.0 * abs (syy * froot) errz = 1.0 + 4.0 * abs (szz * froot) ratx = abs (denx) / errx raty = abs (deny) / erry ratz = abs (denz) / errz ratr = min (ratx, raty, ratz) if (ratx .eq. ratr) then argx = -(sc + syy * cy(1)**2 + szz * cz(1)**2) / sxx cx(1) = sqrt (max (0.0, argx)) if (aax .lt. 0.0) cx(1) = -cx(1) endif if (raty .eq. ratr) then argy = -(sc + szz * cz(1)**2 + sxx * cx(1)**2) / syy cy(1) = sqrt (max (0.0, argy)) if (aay .lt. 0.0) cy(1) = -cy(1) endif if (ratz .eq. ratr) then argz = -(sc + sxx * cx(1)**2 + syy * cy(1)**2) / szz cz(1) = sqrt (max (0.0, argz)) if (aaz .lt. 0.0) cz(1) = -cz(1) endif c.... Find the distance from the proximal point to the surface. dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) else ! Failed to find root. if (nprox .eq. 0) go to 710 ! Found normal intercept. c.... Use approximation for proximal point. arg = -sc - szz * aaz**2 cx(1) = aax * sqrt (arg / sxx) / aar cy(1) = aay * sqrt (arg / syy) / aar cz(1) = aaz dist(1) = sqrt ((cx(1) - aax)**2 + (cy(1) - aay)**2) endif ! Tested for real root. endif ! Tested for special or general case. c.... Type 13b. Elliptic hyperboloid of one sheet: c.... -|sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. c=======================================================================******** c.... End of logic for hyperboloid of one sheet. endif ! Tested for circular or elliptic. go to 510 c.... Type 13. Hyperboloid of one sheet: c.... -|sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. c=======================================================================******** c.... Type 14. Hyperboloid of two sheets: c.... |sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. 114 aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. c=======================================================================******** c.... Type 14a. Circular hyperboloid of two sheets. c.... |sc| + sxx * (x**2 + y**2) - |szz| * z**2 = 0. if (abs (sxx - syy) .le. tol * (sxx + syy)) then ! Circular. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'CIRCULAR HYPERBOLOID OF TWO SHEETS.')") cbugcbugc***DEBUG ends. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if ((aar .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" is at origin. nprox = 2 ! Exact. At both vertices. cx(1) = 0.0 cy(1) = 0.0 cz(1) = sqrt (-sc / szz) dist(1) = cz(1) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) elseif (aar .eq. 0.0) then ! Point "a" is on z axis. cz(1) = aaz / (1.0 - szz / sxx) arg = -(sc + szz * cz(1)**2) / sxx if (arg .le. 0.0) then ! Proximal point at nearest vertex. nprox = 1 ! Exact. cx(1) = 0.0 cy(1) = 0.0 cz(1) = sqrt (-sc / szz) if (aaz .lt. 0.0) cz(1) = -cz(1) dist(1) = cz(1) - aaz else ! Proximal ring. Pick 3 points. nprox = 3 ! Exact. 3 points on the proximal circle. rc = sqrt (arg) ! Radius of proximal circle. cx(1) = rc cy(1) = 0.0 dist(1) = sqrt (arg + (cz(1) - aaz)**2) cx(2) = -0.5 * rc cy(2) = 0.5 * sqrt (3.0) * rc cz(2) = cz(1) dist(2) = dist(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(3) = dist(1) endif elseif (aaz .eq. 0.0) then ! Two proximal points, off-axis. nprox = 2 ! Exact. rc = aar / (1.0 - sxx / szz) ! Radius of proximal points. arg = -(sc + sxx * rc**2) / szz cx(1) = aax * rc / aar cy(1) = aay * rc / aar cz(1) = sqrt (arg) dist(1) = sqrt ((rc - aar)**2 + arg) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) else ! General fourth order solution. c.... Find the roots of the equation for |d| / |N|. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. c0 = sc + sxx * aar**2 + szz * aaz**2 c1 = -4.0 * (sc * (sxx + szz) + sxx * szz * (aar**2 + aaz**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aar**2 + sxx * aaz**2)) c3 = -16.0 * sc * sxx * szz * (sxx + szz) c4 = 16.0 * sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 320 n = 1, nroot denr = abs (1.0 - 2.0 * sxx * rreal(n) + fuz) cxx = aax / denr cyy = aay / denr denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cxx - aax)**2 + (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 320 continue endif ! Tested for symmetry condition. c.... Type 14a. Circular hyperboloid of two sheets. c.... |sc| + sxx * (x**2 + y**2) - |szz| * z**2 = 0. c=======================================================================******** c.... Middle of logic for hyperboloid of two sheets. else ! Elliptic hyperboloid of two sheets. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'ELLIPTIC HYPERBOLOID OF TWO SHEETS.')") cbugcbugc***DEBUG ends. c=======================================================================******** c.... Type 14b. Elliptic hyperboloid of two sheets: c.... |sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if ((aar .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" at origin. nprox = 2 ! Exact, at z vertices. cx(1) = 0.0 cy(1) = 0.0 cz(1) = sqrt (-sc / szz) dist(1) = cz(1) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) elseif ((aax .eq. 0.0) .and. (aay .eq. 0.0)) then ! Point "a" on z axis. czz = aaz / (1.0 - szz / sxx) arg = -(sc + szz * czz**2) / sxx if (arg .le. 0.0) then nprox = 1 ! Exact, at nearest z vertex. cx(1) = 0.0 cy(1) = 0.0 cz(1) = (aaz / abs (aaz)) * sqrt (-sc / szz) dist(1) = cz(1) - aaz else nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = 0.0 cz(1) = czz dist(1) = sqrt (arg + (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) endif ! Tested sign of arg. elseif ((aay .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" on x axis. cxx = aax / (1.0 - sxx / szz) arg = -(sc + sxx * cxx**2) / szz nprox = 2 ! Exact. cx(1) = cxx cy(1) = 0.0 cz(1) = sqrt (arg) dist(1) = sqrt (arg + (cx(1) - aax)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) elseif ((aax .eq. 0.0) .and. (aaz .eq. 0.0)) then ! Point "a" on y axis. cyy = aay / (1.0 - syy / szz) arg = -(sc + syy * cyy**2) / szz nprox = 2 ! Exact, at non-zero cz. cx(1) = 0.0 cy(1) = cyy cz(1) = sqrt (arg) dist(1) = sqrt (arg + (cy(1) - aay)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) elseif (aax .eq. 0.0) then ! In plane x = 0. cyy = aay / (1.0 - syy / sxx) czz = aaz / (1.0 - szz / sxx) arg = -(sc + syy * cyy**2 + szz * czz**2) / sxx if (arg .ge. 0.0) then nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = czz dist(1) = sqrt (arg + (cy(1) - aay)**2 + (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Sign of arg is negative. c0 = sc + syy * aay**2 + szz * aaz**2 c1 = -4.0 * (sc * (syy + szz) + syy * szz * (aay**2 + aaz**2)) c2 = 4.0 * (sc * (syy**2 + 4.0 * syy * szz + szz**2) + & syy * szz * (szz * aay**2 + syy * aaz**2)) c3 = -16.0 * sc * syy * szz * (syy + szz) c4 = 16 * sc * syy**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 360 n = 1, nroot deny = 1.0 - 2.0 * syy * rreal(n) + fuz cyy = aay / deny denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = 0.0 cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 360 continue endif ! Tested sign of arg. elseif (aay .eq. 0.0) then ! In plane y = 0. c0 = sc + sxx * aax**2 + szz * aaz**2 c1 = -4.0 * (sc * (sxx + szz) + sxx * szz * (aax**2 + aaz**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aax**2 + sxx * aaz**2)) c3 = -16.0 * sc * sxx * szz * (sxx + szz) c4 = 16 * sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 370 n = 1, nroot denx = 1.0 - 2.0 * sxx * rreal(n) + fuz cxx = aax / denx denz = 1.0 - 2.0 * szz * rreal(n) + fuz czz = aaz / denz distsq = (cxx - aax)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = 0.0 cz(1) = czz dist(1) = sqrt (distsq) endif 370 continue elseif (aaz .eq. 0.0) then ! In plane z = 0. cxx = aax / (1.0 - sxx / szz) cyy = aay / (1.0 - syy / szz) arg = -(sc + sxx * cxx**2 + syy * cyy**2) / szz nprox = 2 ! Exact. cx(1) = cxx cy(1) = cyy cz(1) = sqrt (arg) dist(1) = sqrt (arg + (cx(1) - aax)**2 + (cy(1) - aay)**2) cx(2) = cx(1) cy(2) = cy(1) cz(2) = -cz(1) dist(2) = dist(1) else ! General 6th order equation. c.... General elliptic hyperboloid of two sheets. c.... Find the coeffients of the 6th-order equation for f. c.... f = (cx - aax) / (2.0 * sxx * cx) (aax not zero) c.... f = (cy - aay) / (2.0 * syy * cy) (aay not zero) c.... f = (cz - aaz) / (2.0 * szz * cz) (aaz not zero) fsum1 = sxx + syy + szz fsum2 = sxx * syy + syy * szz + szz * sxx fsum3 = sxx * syy * szz a(0) = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 a(1) = -4.0 * (sc * fsum1 + & sxx * (syy + szz) * aax**2 + & syy * (szz + sxx) * aay**2 + & szz * (sxx + syy) * aaz**2 ) a(2) = 4.0 * (sc * (fsum1**2 + 2.0 * fsum2) + & sxx * (syy**2 + 4.0 * syy * szz + szz**2) * aax**2 + & syy * (szz**2 + 4.0 * szz * sxx + sxx**2) * aay**2 + & szz * (sxx**2 + 4.0 * sxx * syy + syy**2) * aaz**2 ) a(3) = -16.0 * (sc * (fsum1 * fsum2 + fsum3) + & fsum3 * ((syy + szz) * aax**2 + & (szz + sxx) * aay**2 + & (sxx + syy) * aaz**2 )) a(4) = 16.0 * (sc * (fsum2**2 + 2.0 * fsum1 * fsum3) + & fsum3 * syy * szz * aax**2 + & fsum3 * szz * sxx * aay**2 + & fsum3 * sxx * syy * aaz**2 ) a(5) = -64.0 * sc * fsum2 * fsum3 a(6) = 64.0 * sc * fsum3**2 c.... Find limits on f for an elliptic hyperboloid of two sheets. c.... Find out if external point is inside or outside. side = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 if (side .le. 0.0) then ! Inside. argx = -(sc + syy * aay**2 + szz * aaz**2) / sxx xint = sqrt (argx) fmaxx = 0.5 * (1.0 - abs (aax) / xint) / sxx argy = -(sc + szz * aaz**2 + sxx * aax**2) / syy yint = sqrt (argy) fmaxy = 0.5 * (1.0 - abs (aay) / yint) / syy argz = -(sc + sxx * aax**2 + syy * aay**2) / szz zint = sqrt (argz) fmaxz = 0.5 * (1.0 - abs (aaz) / zint) / szz fmin = 0.0 - 2.0 * tolx fmax = min (fmaxx, fmaxy, fmaxz) + 2.0 * tolx finit = fmax else ! Outside. argx = -(sc + syy * aay**2 + szz * aaz**2) / sxx if (argx .gt. 0.0) then xint = sqrt (argx) fminx = 0.5 * (1.0 - abs (aax) / xint) / sxx else fminx = -1.e99 endif argy = -(sc + szz * aaz**2 + sxx * aax**2) / syy if (argy .gt. 0.0) then yint = sqrt (argy) fminy = 0.5 * (1.0 - abs (aay) / yint) / syy else fminy = -1.e99 endif argz = -(sc + sxx * aax**2 + syy * aay**2) / szz zint = sqrt (argz) fminz = 0.5 * (1.0 - abs (aaz) / zint) / szz fmin = max (fminx, fminy, fminz) - 2.0 * tolx fmax = 0.0 + 2.0 * tolx finit = fmin endif ! Tested inside or outside. c.... Find froot, value of f at proximal point. call aptnewt (a, 6, finit, fmin, fmax, 1000, tol, & froot, nerra) if (nerra .eq. 0) then nprox = 1 denx = 1.0 - 2.0 * sxx * froot + fuz deny = 1.0 - 2.0 * syy * froot + fuz denz = 1.0 - 2.0 * szz * froot + fuz cx(1) = aax / (denx + fuz) cy(1) = aay / (deny + fuz) cz(1) = aaz / (denz + fuz) c.... Make a correction for an initial coordinate very close to an axis. errx = 1.0 + 4.0 * abs (sxx * froot) erry = 1.0 + 4.0 * abs (syy * froot) errz = 1.0 + 4.0 * abs (szz * froot) ratx = abs (denx) / errx raty = abs (deny) / erry ratz = abs (denz) / errz ratr = min (ratx, raty, ratz) if (ratx .eq. ratr) then argx = -(sc + syy * cy(1)**2 + szz * cz(1)**2) / sxx cx(1) = sqrt (max (0.0, argx)) if (aax .lt. 0.0) cx(1) = -cx(1) endif if (raty .eq. ratr) then argy = -(sc + szz * cz(1)**2 + sxx * cx(1)**2) / syy cy(1) = sqrt (max (0.0, argy)) if (aay .lt. 0.0) cy(1) = -cy(1) endif if (ratz .eq. ratr) then argz = -(sc + sxx * cx(1)**2 + syy * cy(1)**2) / szz cz(1) = sqrt (max (0.0, argz)) if (aaz .lt. 0.0) cz(1) = -cz(1) endif c.... Find the distance from the proximal point to the surface. dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) else ! Failed to find root. if (nprox .eq. 0) go to 710 ! Found normal intercept. cz(1) = sqrt (-sc / szz) if (abs (aaz) .le. cz(1)) then ! Pick nearest vertex. cx(1) = 0.0 cy(1) = 0.0 if (aaz .lt. 0.0) cz(1) = -cz(1) dist(1) = sqrt (aax**2 + aay**2 + (cz(1) - aaz)**2) else ! Use a point on ellipse at z = aaz. rcx = sqrt (-(sc + szz * aaz**2) / sxx) rcy = sqrt (-(sc + szz * aaz**2) / syy) cx(1) = aax * rcx / aar cy(1) = aay * rcy / aar cz(1) = aaz dist(1) = sqrt ((cx(1) - aax)**2 + (cy(1) - aay)**2) endif ! Tested aaz. endif ! Tested for real root. endif ! Tested for special and general cases. c.... Type 14b. Elliptic hyperboloid of two sheets: c.... |sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. c=======================================================================******** c.... End of logic for hyperboloid of two sheets. endif ! Tested for circular or elliptic. go to 510 c.... Type 14. Hyperboloid of two sheets: c.... |sc| + sxx * x**2 + |syy| * y**2 - |szz| * z**2 = 0. c=======================================================================******** c.... Type 15. Real ellipsoid: c.... -|sc| + sxx * x**2 + |syy| * y**2 + |szz| * z**2 = 0. c.... sxx => syy => szz. c.... Distance of point "a" from center of ellipsoid. 115 aar = sqrt (aax**2 + aay**2 + aaz**2) c=======================================================================******** c.... Beginning of logic for prolate spheroid. Symmetric around z axis. if (abs (sxx - syy) .le. tol * (sxx + syy)) then ! Prolate spheroid. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'PROLATE SPHEROID.')") cbugcbugc***DEBUG ends. zext = sqrt (-sc / szz) ! Limiting z values (+ and -). rext = sqrt (-sc / sxx) ! Extreme radius at z = 0. aar = sqrt (aax**2 + aay**2) ! Distance of point "a" from the z axis. if ((aaz .eq. 0.0) .and. (aar .eq. 0.0)) then ! Point "a" at center. c.... Use three points on waist, at angles of 0, 120, 240 degrees. nprox = 3 ! Exact. On circular waist. cx(1) = rext cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) cx(2) = -0.5 * cx(1) cy(2) = 0.5 * sqrt (3.0) * cx(1) cz(2) = cz(1) dist(2) = dist(1) cx(3) = cx(2) cy(3) = -cy(2) cz(3) = cz(1) dist(3) = dist(1) go to 510 elseif (aaz .eq. 0.0) then ! Point "a" is off center in plane z = 0. nprox = 1 ! Exact. cx(1) = aax * rext / aar cy(1) = aay * rext / aar cz(1) = 0.0 dist(1) = rext - aar go to 510 elseif (aar .eq. 0.0) then ! Point "a" on the z axis. zcen = (1.0 - szz / sxx) * zext if (abs (aaz) .ge. zcen) then ! Point "a" is near vertex. nprox = 1 ! Exact. At vertex of ellipsoid. cx(1) = 0.0 cy(1) = 0.0 cz(1) = zext * aaz / abs (aaz) dist(1) = cz(1) - aaz go to 510 else ! Point "a" is not near vertex. nprox = 3 ! Exact. On circular ring. zc = aaz / (1.0 - szz / sxx) rc = sqrt (-(sc + szz * zc**2) / sxx) ! Radius of circle. cx(1) = rc cy(1) = 0.0 cz(1) = zc dist(1) = sqrt (cx(1)**2 + (cz(1) - aaz)**2) cx(2) = -0.5 * cx(1) cy(2) = 0.5 * sqrt (3.0) * cx(1) cz(2) = zc dist(2) = dist(1) cx(3) = -0.5 * cx(1) cy(3) = -0.5 * sqrt (3.0) * cx(1) cz(3) = zc dist(3) = dist(1) go to 510 endif ! Tested point "a" nearness to vertex. else ! Point "a" is not in symmetric position. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aar**2 + szz * aaz**2 c1 = -2.0 * (sc * (sxx + szz) + sxx * szz * (aar**2 + aaz**2)) c2 = sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aar**2 + sxx * aaz**2) c3 = -2.0 * sc * sxx * szz * (sxx + szz) c4 = sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the limiting values of the roots. fmaxr = (1.0 - abs (aar) / rext) / sxx fmaxz = (1.0 - abs (aaz) / zext) / szz fmax = min (fmaxr, fmaxz) cbugcbugc***DEBUG begins. cbugcbug 9971 format ('fmaxr, z, rz=',1p3e22.14) cbugcbug write ( 3, 9971) fmaxr, fmaxz, fmax cbugcbugc***DEBUG ends. c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 230 n = 1, nroot denz = 1.0 - szz * rreal(n) + fuz czz = aaz / denz denr = abs (1.0 - sxx * rreal(n) + fuz) crrs = aar / denr arg = -(sc + szz * czz**2) / sxx if (arg .ge. 0.0) then crr = sqrt (arg) else crr = crrs endif cbugcbugc***DEBUG begins. cbugcbug 9967 format ('crrs, crr =',1p2e22.14) cbugcbug write ( 3, 9966) crrs, crr cbugcbugc***DEBUG ends. cxx = aax / denr cyy = aay / denr distsq = (cxx - aax)**2 + (cyy - aay)**2 + (czz - aaz)**2 if ((distsq .lt. distsqm) .and. (rreal(n) .le. fmax)) then cbugcbugc***DEBUG begins. cbugcbug write ( 3, "('DEBUG: ACCEPTED')") cbugcbugc***DEBUG ends. distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif cbugcbugc***DEBUG begins. cbugcbug 9961 format (/,'n=',i2,' rreal =',1pe22.14,' fmax =',1pe22.14 / cbugcbug & 'denr, denz =',1p2e22.14 / cbugcbug & 'cxx,cyy,czz =',1p3e22.14 / cbugcbug & 'distsq =',1pe22.14 ) cbugcbug write ( 3, 9961) n, rreal(n), fmax, denr, denz, cbugcbug & cxx, cyy, czz, distsq cbugcbug rrealx = (cxx - aax) / (sxx * cxx + fuz) cbugcbug rrealy = (cyy - aay) / (syy * cyy + fuz) cbugcbug rrealz = (czz - aaz) / (szz * czz + fuz) cbugcbug qside = sc + sxx * cxx**2 + syy * cyy**2 + szz * czz**2 cbugcbug 9962 format ( 'rrealx,y,z =',1p3e22.14 / cbugcbug & 'qside =',1pe22.14 ) cbugcbug write ( 3, 9962) rrealx, rrealy, rrealz, qside cbugcbugc***DEBUG ends. 230 continue go to 510 endif ! Tested for symmetric point "a". c.... End of logic for prolate spheroid. c=======================================================================******** c.... Beginning of logic for oblate spheroid. Symmetric around x axis. elseif (abs (syy - szz) .le. tol * (syy + szz)) then ! Oblate spheroid. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'OBLATE SPHEROID.')") cbugcbugc***DEBUG ends. xext = sqrt (-sc / sxx) ! Limiting x values (+ and -). rext = sqrt (-sc / syy) ! Extreme radius at x = 0. aar = sqrt (aay**2 + aaz**2) ! Distance of point "a" from the x axis. if ((aax .eq. 0.0) .and. (aar .eq. 0.0)) then ! Point "a" is at origin. nprox = 2 ! Exact. At extreme x values. cx(1) = xext cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) cx(2) = -xext cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 elseif (aar .eq. 0.0) then ! On x axis. nprox = 1 cx(1) = sign (xext, aax) cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) - aax go to 510 elseif (aax .eq. 0.0) then ! Point "a" is off center in plane x = 0. rcen = (1.0 - syy / sxx) * rext if (aar .ge. rcen) then ! Near sharpest tip. nprox = 1 ! Exact. cx(1) = 0 cy(1) = aay * rext / aar cz(1) = aaz * rext / aar dist(1) = rext - aar go to 510 elseif (aar .lt. rcen) then ! Not near sharpest tip. nprox = 2 ! Exact. cr = aar / (1.0 - syy / sxx) cx(1) = sqrt (-(sc + syy * cr**2) / sxx) cy(1) = aay * cr / aar cz(1) = aaz * cr / aar dist(1) = sqrt (cx(1)**2 + (cr - aar)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) go to 510 endif ! Tested nearness to waist. else ! Point "a" not in symmetric position. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aax**2 + syy * aar**2 c1 = -2.0 * (sc * (sxx + syy) + sxx * syy * (aax**2 + aar**2)) c2 = sc * (sxx**2 + 4.0 * sxx * syy + syy**2) + & sxx * syy * (syy * aax**2 + sxx * aar**2) c3 = -2.0 * sc * sxx * syy * (sxx + syy) c4 = sc * sxx**2 * syy**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the limiting values of the roots. fmaxx = (1.0 - abs (aax) / xext) / sxx fmaxr = (1.0 - abs (aar) / rext) / syy fmax = min (fmaxx, fmaxr) cbugcbugc***DEBUG begins. cbugcbug 9965 format ('fmaxx, r, xr=',1p3e22.14) cbugcbug write ( 3, 9965) fmaxx, fmaxr, fmax cbugcbugc***DEBUG ends. c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 240 n = 1, nroot denr = abs (1.0 - syy * rreal(n) + fuz) crr = aar / denr cyy = aay / denr czz = aaz / denr denx = 1.0 - sxx * rreal(n) + fuz cxxs = aax / denx arg = -(sc + syy * crr**2) / sxx if (arg .ge. 0.0) then cxx = sign (sqrt (arg), aax) else cxx = cxxs endif cbugcbugc***DEBUG begins. cbugcbug 9966 format ('cxxs, cxx =',1p2e22.14) cbugcbug write ( 3, 9966) cxxs, cxx cbugcbugc***DEBUG ends. distsq = (cxx - aax)**2 + (cyy - aay)**2 + (czz - aaz)**2 if ((distsq .lt. distsqm) .and. (rreal(n) .le. fmax)) then cbugcbugc***DEBUG begins. cbugcbug write ( 3, "('DEBUG: ACCEPTED')") cbugcbugc***DEBUG ends. distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif cbugcbugc***DEBUG begins. cbugcbug write ( 3, 9961) n, rreal(n), fmax, denr, denz, cbugcbug & cxx, cyy, czz, distsq cbugcbug rrealx = (cxx - aax) / (sxx * cxx + fuz) cbugcbug rrealy = (cyy - aay) / (syy * cyy + fuz) cbugcbug rrealz = (czz - aaz) / (szz * czz + fuz) cbugcbug qside = sc + sxx * cxx**2 + syy * cyy**2 + szz * czz**2 cbugcbug write ( 3, 9962) rrealx, rrealy, rrealz, qside cbugcbugc***DEBUG ends. 240 continue endif ! Tested for symmetric point "a". c.... End of logic for oblate spheroid. c=======================================================================******** c.... Beginning of logic for general ellipsoid. else ! General ellipsoid. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'GENERAL ELLIPSOID.')") cbugcbugc***DEBUG ends. if (aar .eq. 0.0) then ! Point "a" is at center of ellipsoid. nprox = 2 ! Exact. cx(1) = sqrt (-sc / sxx) cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) elseif ((aax .eq. 0.0) .and. (aay .eq. 0.0)) then ! On z axis. czz = aaz / (1.0 - szz / sxx) arg = -(sc + szz * czz**2) / sxx if (arg .gt. 0.0) then nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = 0.0 cz(1) = czz dist(1) = sqrt (arg + (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Arg is not positive. nprox = 1 ! Exact, at z vertex. cx(1) = 0.0 cy(1) = 0.0 cz(1) = (aaz / abs (aaz)) * sqrt (-sc / szz) dist(1) = cz(1) - aaz endif ! Tested sign of arg. elseif ((aay .eq. 0.0) .and. (aaz .eq. 0.0)) then ! On x axis. nprox = 1 ! Exact, at nearest x vertex.j cx(1) = (aax / abs (aax)) * sqrt (-sc / sxx) cy(1) = 0.0 cz(1) = 0.0 dist(1) = cx(1) - aax elseif ((aax .eq. 0.0) .and. (aaz .eq. 0.0)) then ! On y axis. cyy = aay / (1.0 - syy / sxx) arg = -(sc + syy * cyy**2) / sxx if (arg .gt. 0.0) then nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = 0.0 dist(1) = sqrt (arg + (cy(1) - aay)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Arg is not positive. nprox = 1 ! Exact, at nearest y vertex. cx(1) = 0.0 cy(1) = (aay / abs (aay)) * sqrt (-sc / syy) cz(1) = 0.0 dist(1) = cy(1) - aay endif ! Tested sign of arg. elseif (aaz .eq. 0.0) then ! In plane z = 0. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aax**2 + syy * aay**2 c1 = -4.0 * (sc * (sxx + syy) + sxx * syy * (aax**2 + aay**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * syy + syy**2) + & sxx * syy * (syy * aax**2 + sxx * aay**2)) c3 = -16.0 * sc * sxx * syy * (sxx + syy) c4 = 16.0 * sc * sxx**2 * syy**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 410 n = 1, nroot cxx = aax / (1.0 - 2.0 * sxx * rreal(n) + fuz) cyy = aay / (1.0 - 2.0 * syy * rreal(n) + fuz) distsq = (cxx - aax)**2 + (cyy - aay)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = cyy cz(1) = 0.0 dist(1) = sqrt (distsq) endif 410 continue elseif (aay .eq. 0.0) then ! In plane y = 0. c.... Find the roots of the equation for |d| / |N|. c0 = sc + sxx * aax**2 + szz * aaz**2 c1 = -4.0 * (sc * (sxx + szz) + sxx * szz * (aax**2 + aaz**2)) c2 = 4.0 * (sc * (sxx**2 + 4.0 * sxx * szz + szz**2) + & sxx * szz * (szz * aax**2 + sxx * aaz**2)) c3 = -16.0 * sc * sxx * szz * (sxx + szz) c4 = 16 * sc * sxx**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 420 n = 1, nroot cxx = aax / (1.0 - 2.0 * sxx * rreal(n) + fuz) czz = aaz / (1.0 - 2.0 * szz * rreal(n) + fuz) distsq = (cxx - aax)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = cxx cy(1) = 0.0 cz(1) = czz dist(1) = sqrt (distsq) endif 420 continue elseif (aax .eq. 0.0) then ! In plane x = 0. cyy = aay / (1.0 - syy / sxx) czz = aaz / (1.0 - szz / sxx) arg = -(sc + syy * cyy**2 + szz * czz**2) / sxx if (arg .gt. 0.0) then nprox = 2 ! Exact. cx(1) = sqrt (arg) cy(1) = cyy cz(1) = czz dist(1) = sqrt (arg + (cy(1) - aay)**2 + (cz(1) - aaz)**2) cx(2) = -cx(1) cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) else ! Arg not positive. c.... Find the roots of the equation for |d| / |N|. c0 = sc + syy * aay**2 + szz * aaz**2 c1 = -4.0 * (sc * (syy + szz) + & syy * szz * (aay**2 + aaz**2)) c2 = 4.0 * (sc * (syy**2 + 4.0 * syy * szz + szz**2) + & syy * szz * (szz * aay**2 + syy * aaz**2)) c3 = -16.0 * sc * syy * szz * (syy + szz) c4 = 16 * sc * syy**2 * szz**2 c0 = c0 / c4 c1 = c1 / c4 c2 = c2 / c4 c3 = c3 / c4 call aptquar (c0, c1, c2, c3, tol, & nroot, nrooti, rreal, rimag, atype, itrun) if (nroot .eq. 0) then ! Logical error. nerr = 3 go to 710 endif c.... Find the proximal point at the minimum distance from point "a". nprox = 1 ! Exact. distsqm = 1.e99 do 430 n = 1, nroot cyy = aay / (1.0 - 2.0 * syy * rreal(n) + fuz) czz = aaz / (1.0 - 2.0 * szz * rreal(n) + fuz) distsq = (cyy - aay)**2 + (czz - aaz)**2 if (distsq .lt. distsqm) then distsqm = distsq cx(1) = 0.0 cy(1) = cyy cz(1) = czz dist(1) = sqrt (distsq) endif 430 continue endif ! Tested sign of arg. else ! Point "a" not on symmetry plane. c.... General ellipsoid. c.... Find the coeffients of the 6th-order equation for f. c.... f = (cx - aax) / (2.0 * sxx * cx) (aax not zero) c.... f = (cy - aay) / (2.0 * syy * cy) (aay not zero) c.... f = (cz - aaz) / (2.0 * szz * cz) (aaz not zero) fsum1 = sxx + syy + szz fsum2 = sxx * syy + syy * szz + szz * sxx fsum3 = sxx * syy * szz a(0) = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 a(1) = -4.0 * (sc * fsum1 + & sxx * (syy + szz) * aax**2 + & syy * (szz + sxx) * aay**2 + & szz * (sxx + syy) * aaz**2 ) a(2) = 4.0 * (sc * (fsum1**2 + 2.0 * fsum2) + & sxx * (syy**2 + 4.0 * syy * szz + szz**2) * aax**2 + & syy * (szz**2 + 4.0 * szz * sxx + sxx**2) * aay**2 + & szz * (sxx**2 + 4.0 * sxx * syy + syy**2) * aaz**2 ) a(3) = -16.0 * (sc * (fsum1 * fsum2 + fsum3) + & fsum3 * ((syy + szz) * aax**2 + & (szz + sxx) * aay**2 + & (sxx + syy) * aaz**2 )) a(4) = 16.0 * (sc * (fsum2**2 + 2.0 * fsum1 * fsum3) + & fsum3 * syy * szz * aax**2 + & fsum3 * szz * sxx * aay**2 + & fsum3 * sxx * syy * aaz**2 ) a(5) = -64.0 * sc * fsum2 * fsum3 a(6) = 64.0 * sc * fsum3**2 c.... Find limits on f for an ellipsoid. c.... Find out if external point is inside or outside. side = sc + sxx * aax**2 + syy * aay**2 + szz * aaz**2 if (side .le. 0.0) then ! Inside. xint = sqrt (-(sc + syy * aay**2 + szz * aaz**2) / sxx) yint = sqrt (-(sc + szz * aaz**2 + sxx * aax**2) / syy) zint = sqrt (-(sc + sxx * aax**2 + syy * aay**2) / szz) fmaxx = 0.5 * (1.0 - abs (aax) / xint) / sxx fmaxy = 0.5 * (1.0 - abs (aay) / yint) / syy fmaxz = 0.5 * (1.0 - abs (aaz) / zint) / szz fmin = 0.0 - 2.0 * tolx fmax = min (fmaxx, fmaxy, fmaxz) + 2.0 * tolx finit = fmax else ! Outside. xlim = sqrt (-sc / sxx) ylim = sqrt (-sc / syy) zlim = sqrt (-sc / szz) fmaxx = 0.5 * (1.0 - abs (aax) / xlim) / sxx fmaxy = 0.5 * (1.0 - abs (aay) / ylim) / syy fmaxz = 0.5 * (1.0 - abs (aaz) / zlim) / szz fmax = min (0.0, fmaxx, fmaxy, fmaxz) + 2.0 * tolx fminx = 0.5 * (1.0 - aar / xlim) / sxx fminy = 0.5 * (1.0 - aar / ylim) / syy fminz = 0.5 * (1.0 - aar / zlim) / szz fmin = min (0.0, fminx, fminy, fminz) - 2.0 * tolx finit = fmin endif ! Tested inside or outside. c.... Find froot, value of f at proximal point. call aptnewt (a, 6, finit, fmin, fmax, 1000, tol, & froot, nerra) if (nerra .eq. 0) then nprox = 1 denx = 1.0 - 2.0 * sxx * froot + fuz deny = 1.0 - 2.0 * syy * froot + fuz denz = 1.0 - 2.0 * szz * froot + fuz cx(1) = aax / (denx + fuz) cy(1) = aay / (deny + fuz) cz(1) = aaz / (denz + fuz) c.... Make a correction for an initial coordinate very close to an axis. errx = 1.0 + 4.0 * abs (sxx * froot) erry = 1.0 + 4.0 * abs (syy * froot) errz = 1.0 + 4.0 * abs (szz * froot) ratx = abs (denx) / errx raty = abs (deny) / erry ratz = abs (denz) / errz ratr = min (ratx, raty, ratz) if (ratx .eq. ratr) then argx = -(sc + syy * cy(1)**2 + szz * cz(1)**2) / sxx cx(1) = sqrt (max (0.0, argx)) if (aax .lt. 0.0) cx(1) = -cx(1) endif if (raty .eq. ratr) then argy = -(sc + szz * cz(1)**2 + sxx * cx(1)**2) / syy cy(1) = sqrt (max (0.0, argy)) if (aay .lt. 0.0) cy(1) = -cy(1) endif if (ratz .eq. ratr) then argz = -(sc + sxx * cx(1)**2 + syy * cy(1)**2) / szz cz(1) = sqrt (max (0.0, argz)) if (aaz .lt. 0.0) cz(1) = -cz(1) endif c.... Find the distance from the proximal point to the surface. dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) else ! Failed to find root. if (nprox .eq. 0) go to 710 ! Found normal intercept. c.... Use approximation for proximal point. cx(1) = aax * sqrt (-sc / sxx) / aar cy(1) = aay * sqrt (-sc / syy) / aar cz(1) = aaz * sqrt (-sc / szz) / aar dist(1) = sqrt ((cx(1) - aax)**2 + & (cy(1) - aay)**2 + & (cz(1) - aaz)**2 ) endif ! Tested for real root. endif ! Tested for special case or general. endif ! Tested type of ellipsoid. go to 510 c.... Type 15. Real ellipsoid: c.... -|sc| + sxx * x**2 + |syy| * y**2 + |szz| * z**2 = 0. c.... sxx => syy => szz. c=======================================================================******** c.... Type 16. Real sphere: -|sc| + sxx * (x**2 + y**2 + z**2) = 0. 116 nprox = 1 ! Exact. cbugcbugc***DEBUG begins. cbugcbug write ( 3, "(/ 'SPHERE.')") cbugcbugc***DEBUG ends. ra = sqrt (aax**2 + aay**2 + aaz**2) rc = sqrt (-sc / sxx) ! Radius of sphere. if (ra .ne. 0.0) then ! Point "a" is not at center of sphere. cx(1) = aax * rc / ra cy(1) = aay * rc / ra cz(1) = aaz * rc / ra dist(1) = rc - ra else ! Point "a" is at center of sphere. nprox = 4 ! Exact. 4 points on the proximal sphere. cx(1) = -rc cy(1) = 0.0 cz(1) = 0.0 dist(1) = rc cx(2) = rc cy(2) = cy(1) cz(2) = cz(1) dist(2) = dist(1) cx(3) = 0.0 cy(3) = rc cz(3) = cz(1) dist(3) = dist(1) cx(4) = cx(3) cy(4) = cy(1) cz(4) = rc dist(4) = dist(1) endif ! Tested for central or noncentral "a". go to 510 c.... Type 16. Real sphere: -|sc| + sxx * (x**2 + y**2 + z**2) = 0. c=======================================================================******** c.... Rotate and translate proximal points back to original coordinates. 510 continue nproxs = nprox if (nprox .eq. -1) nproxs = 1 cbugcbugc***DEBUG begins. cbugcbug do 7710 n = 1, nproxs cbugcbug write ( 3, 9906) cx(n), cy(n), cz(n), dist(n) cbugcbug 7710 continue cbugcbugc***DEBUG ends. call aptmopv (rotm, 1, 0., 0., 0., cx, cy, cz, nproxs, tol, nerra) cenx = -cenx ceny = -ceny cenz = -cenz call apttran (cenx, ceny, cenz, cx, cy, cz, nproxs, tol, nerra) c.... Impose the correct sign on dist. 710 continue nproxs = nprox if (nprox .eq. -1) nproxs = 1 do 720 n = 1, nproxs ! Loop over proximal points. vdot = (cx(n) - ax) * bx + & (cy(n) - ay) * by + & (cz(n) - az) * bz dist(n) = sign (dist(n), vdot) coldc.... See if point "c" is on quadric surface. cold cold cside = qc + qx * cx(n) + qy * cy(n) + qz * cz(n) + cold & qxy * cx(n) * cy(n) + cold & qyz * cy(n) * cz(n) + cold & qzx * cz(n) * cx(n) + cold & qxx * cx(n)**2 + qyy * cy(n)**2 + qzz * cz(n)**2 cold cold cerr = tol * (abs (qc) + cold & 2.0 * (abs (qx * cx(n)) + cold & abs (qy * cy(n)) + cold & abs (qz * cz(n))) + cold & 3.0 * (abs (qxy * cx(n) * cy(n)) + cold & abs (qyz * cy(n) * cz(n)) + cold & abs (qzx * cz(n) * cx(n)) + cold & abs (qxx) * cx(n)**2 + cold & abs (qyy) * cy(n)**2 + cold & abs (qzz) * cz(n)**2 )) cold cold if (abs (cside) .gt. cerr) then cbugcbugcoldc***DEBUG begins. cbugcbugcold 9951 format ('DEBUG cside, cerr=',1p2e22.14) cbugcbugcold write ( 3, 9951) cside, cerr cbugcbugcoldc***DEBUG ends. cold nerr = 4 cold endif 720 continue ! End of loop over proximal points. cbugc***DEBUG begins. cbug 9905 format (/ 'aptwhis results: nprox =',i3,' nerr=',i3 / cbug & ' dside =',1pe22.14 / cbug & ' bx, by, bz =',1p3e22.14 / cbug & ' cx, cy, cz =',1p3e22.14 / cbug & ' dist =',1pe22.14 ) cbug 9906 format ( cbug & ' cx, cy, cz =',1p3e22.14 / cbug & ' dist =',1pe22.14 ) cbug write ( 3, 9905) nprox, nerr, dside, bx, by, bz, cbug & cx(1), cy(1), cz(1), dist(1) cbug do 7720 n = 2, nproxs cbug write ( 3, 9906) cx(n), cy(n), cz(n), dist(n) cbug 7720 continue cbugc***DEBUG ends. return c.... End of subroutine aptwhis. (+1 line.) end UCRL-WEB-209832