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