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