subroutine aptqext (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
     &                    tol, ntype, cenx, ceny, cenz, np, ex, ey, ez,
     &                    adir, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTQEXT
c
c     call aptqext (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
c                   tol, ntype, cenx, ceny, cenz, np, ex, ey, ez,
c                   adir, nerr)
c
c     Version:  aptqext  Updated    1997 December 2 18:05.
c               aptqext  Originated 1994 April 14 17:00.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  For the quadric surface specified by the coefficients qc to qzz,
c               to find the surface type ntype, the center of symmetry
c               (cenx, ceny, cenz), and for any extreme values of x, y or z on
c               the surface, to find the coordinates (ex(n), ey(n), ez(n)), the
c               type adir(n), n = 1, np, and the total number np (0 to 6).
c               The possible types of extrema are minima and maxima (points,
c               lines or planes), saddle points, constant planes, intersection
c               lines, and isolated points.
c 
c               The implicit equation of the quadric surface is:
c
c                 F(x,y,z) = qc + qx  * x     + qy  * y     +
c                                 qxy * x * y + qyz * y * z + qzx * z * x +
c                                 qxx * x**2  + qyy * y**2  + qzz * z**2  = 0  .
c
c               Results may be unpredictable if the coefficients of the quadric
c               surface differ by many orders of magnitude.
c
c     Method:
c
c               The components of the normal vector to the surface, N, are the
c                 partial derivatives of F with respect to x, y and z:
c
c                 N = (Nx, Ny, Nz) = (df/dx, df/dy, df/dz)
c                 Nx = qx + 2.0 * qxx * x + qxy * y + qzx * z
c                 Ny = qy + qxy * x + 2.0 * qyy * y + qyz * z
c                 Nz = qz + qzx * x + qyz * y + 2.0 * qzz * z
c
c               Extreme values of x:  none if Nx is always zero.  Otherwise,
c               from zero to two where Ny and Nz are both zero, and Nx is not.
c               Extreme values of y:  none if Ny is always zero.  Otherwise,
c               from zero to two where Nz and Nx are both zero, and Ny is not.
c               Extreme values of z:  none if Nz is always zero.  Otherwise,
c               from zero to two where Nx and Ny are both zero, and Nz is not.
c
c               General cases:
c
c               Null surfaces (all coefficients zero):  no extremes (nerr = 1).
c               Skew planes (simple, coincident or parallel): no extremes.
c               Simple or coincident planes parallel to a major plane:
c                 one constant plane of x, y or z, with the other two
c                 coordinates arbitrary.
c               Parallel planes parallel to a major plane:
c                 minimum and maximum planes of x, y or z, with the other two
c                 coordinates arbitrary.
c               Intersecting planes, with the intersection perpendicular to a
c                 major axis:  one intersection line, at a coordinate on the
c                 major axis.  The non-extreme coordinates nay be anywhere on
c                 the line.
c               Cylinders:  zero, two or four extreme lines if the axis is
c                 skew, in a major plane, or on a major axis, respectively.
c                 The non-extreme coordinates may be anywhere on each line.
c               Spheres or ellipsoids:  six extreme points (two x, two y,
c                 two z).
c               Cones:  no extreme points, but the vertex resembles one.
c               Hyperboloids of one sheet:  zero, two or four saddle points.
c               Hyperboloids of two sheets:  zero, two or four extreme points.
c                 With two or four extreme points, the minimum of one sheet
c                 will be larger than the maximum of the other sheet.
c               Elliptic paraboloids:  one extreme point if aligned with an
c                 axis, and three extreme points (one x, one y, one z), if not.
c               Hyperbolic paraboloids:  one saddle point, if aligned with an
c                 axis, and three saddle points (one x, one y, one z) if not.
c
c     Input:    qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol.
c
c     Output:   np, ex, ey, ez, adir, nerr.
c
c     Calls: aptaxis, aptcris, aptplpl, aptqexc, aptqrts, aptvadd 
c               
c
c     Glossary:
c
c     adir      Output   The type of extreme point:
c                          Minimum:  "Min x", "Min y" or "Min z".
c                          Maximum:  "Max x", "Max y" or "Max z".
c                          Aligned planes: "Const x", "Const y" or "Const z".
c                          Intersection of two planes, perpendicular to a major
c                          axis:
c                          "Inter x", "Inter y" or "Inter z".
c                          Saddle point:  "Saddle x", Saddle y" or "Saddle z".
c                          Singular isolated point:  "Point x", "Point y" or
c                          "Point z".
c                          If none, return "unknown".  Size = 6.
c
c     cenx,y,z  Output   The x, y and z coordinates of the center of symmetry
c                          of the quadric surface.
c
c     ex,ey,ez  Output   The x, y and z coordinates of an extreme point, a
c                          saddle point, or an intersection of two planes, for
c                          which the x, y or z coordinate is a maximum or a
c                          minimum.
c                          If none, return (-1.e99,-1.e99,-1.e99).  Size = 6.
c
c     nerr      Output   Indicates an input error, if not zero:
c                          1 if all coefficients, except possible qc, are zero,
c                          or the quadric surface type is 17, 19, 20, 23 or 24.
c
c     np        Output   Number of extreme or intersection points found.
c
c     ntype     Output   Indicates the classification of the quadric surface:
c                           0 = Plane.
c                           1 = Coincident planes.
c                           2 = Real parallel planes.
c                           3 = Real intersecting planes.
c                           4 = Parabolic cylinder.
c                           5 = Hyperbolic cylinder.
c                           6 = Real elliptic cylinder.
c                           7 = Real circular cylinder.
c                           8 = Hyperbolic paraboloid.
c                           9 = Elliptic paraboloid.
c                          10 = Circular paraboloid.
c                          11 = Real elliptic cone.
c                          12 = Real circular cone.
c                          13 = Hyperboloid of one sheet.
c                          14 = Hyperboloid of two sheets.
c                          15 = Real ellipsoid.
c                          16 = Real sphere.
c                          17 = Imaginary parallel planes.
c                          18 = Imaginary intersecting planes (line x = y = 0).
c                          19 = Imaginary elliptic cylinder.
c                          20 = Imaginary circular cylinder.
c                          21 = Imaginary elliptic cone (point x = y = z = 0).
c                          22 = Imaginary circular cone (point x = y = z = 0).
c                          23 = Imaginary ellipsoid.
c                          24 = Imaginary sphere.
c
c     qc        In       Constant term in the equation of the quadric surface.
c
c     qx,qy,qz  In       Coefficients of x, y and z, respectively, in the
c                          equation of the quadric surface.
c
c     qxx       In       Coefficient of x**2 in the equation of the quadric
c                          surface.
c
c     qxy       In       Coefficient of x * y in the equation of the quadric
c                          surface.
c
c     qyy       In       Coefficient of y**2 in the equation of the quadric
c                          surface.
c
c     qyz       In       Coefficient of y * z in the equation of the quadric
c                          surface.
c
c     qzx       In       Coefficient of z * x in the equation of the quadric
c                          surface.
c
c     qzz       In       Coefficient of z**2 in the equation of the quadric
c                          surface.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-11.
c                          On 32-bit computers, recommend 1.e-6.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

      dimension ex(6)                 ! Coordinate x of extreme point.
      dimension ey(6)                 ! Coordinate y of extreme point.
      dimension ez(6)                 ! Coordinate z of extreme point.

      dimension adir(6)               ! Normal direction of extreme point.
      character*8 adir                ! Normal direction of extreme point.

c.... Local variables.

      common /laptqext/ rotm (3,3)
      common /laptqext/ ccx(4)        ! Coordinate x of center of curvature.
      common /laptqext/ ccy(4)        ! Coordinate y of center of curvature.
      common /laptqext/ ccz(4)        ! Coordinate z of center of curvature.
      common /laptqext/ ucx(4)        ! Component x of unit vector along curve.
      common /laptqext/ ucy(4)        ! Component y of unit vector along curve.
      common /laptqext/ ucz(4)        ! Component z of unit vector along curve.
      common /laptqext/ rcurve(4)     ! An extreme radius of curvature.

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqext finding the extrema of quadric surface.' /
cbug     &  '  tol=',1pe13.5 )
cbug 9902 format (
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
cbug      write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.

c.... Initialize.

      nerr  = 0
      np    = 0

      do 110 n = 1, 6
        ex(n)   = -1.e99
        ey(n)   = -1.e99
        ez(n)   = -1.e99
        adir(n) = 'unknown'
  110 continue

      tolx = tol
      if (tol .eq. 0.0) tolx = 1.e-11

c.... Test for input errors.

c.... Find the type and center of the quadric surface.

      sc  = qc
      sx  = qx
      sy  = qy
      sz  = qz
      sxy = qxy
      syz = qyz
      szx = qzx
      sxx = qxx
      syy = qyy
      szz = qzz

      call aptaxis (sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz,
     &              tolx, cenx, ceny, cenz, rotm, ntype, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9922 format (
cbugcbug     &  '  sc=    (std)',1pe22.14 /
cbugcbug     &  '  sx,sy,sz=   ',1p3e22.14 /
cbugcbug     &  '  sxy,syz,szx=',1p3e22.14 /
cbugcbug     &  '  sxx,syy,szz=',1p3e22.14)
cbugcbug      write ( 3, 9922) sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz
cbugcbugc***DEBUG ends.

      if ((nerra .ne. 0)   .or.
     &    (ntype .eq. 17)  .or.
     &    (ntype .eq. 19)  .or.
     &    (ntype .eq. 20)  .or.
     &    (ntype .eq. 23)  .or.
     &    (ntype .eq. 24)) then
        nerr = 1
        go to 210
      endif

c.... There are no extreme points on cones, only the vertex point.

      if ((ntype .eq. 11)  .or.
     &    (ntype .eq. 12)) then
        np      = 1
        adir(1) = 'Vertex'
        ex(1)   = cenx
        ey(1)   = ceny
        ez(1)   = cenz
        go to 210
      endif

c=======================================================================********

c.... Test for a cylinder aligned with the x axis.  Only y and z extrema.
c....   Nx = 0 everywhere.

      if ((qx  .eq. 0.0)  .and.
     &    (qxx .eq. 0.0)  .and.
     &    (qxy .eq. 0.0)  .and.
     &    (qzx .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9903 format ('aptqext DEBUG.  X cylinder, y and z extrema.')
cbugcbug      write ( 3, 9903)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qy, qz, qyz, qyy, qzz, tol,
     &                np, ey, ez, adir, nerra)

          do 120 n = 1, np
            ex(n) = 0.0
            if (adir(n) .eq. 'ucon') adir(n) = 'Const y'
            if (adir(n) .eq. 'uint') adir(n) = 'Inter y'
            if (adir(n) .eq. 'umax') adir(n) = 'Max y'
            if (adir(n) .eq. 'umin') adir(n) = 'Min y'
            if (adir(n) .eq. 'vcon') adir(n) = 'Const z'
            if (adir(n) .eq. 'vint') adir(n) = 'Inter z'
            if (adir(n) .eq. 'vmax') adir(n) = 'Max z'
            if (adir(n) .eq. 'vmin') adir(n) = 'Min z'
  120     continue

        go to 210

      endif

c=======================================================================********

c.... Test for a cylinder aligned with the y axis.  Only z and x extrema.
c....   Ny = 0 everywhere.

      if ((qy  .eq. 0.0)  .and.
     &    (qyy .eq. 0.0)  .and.
     &    (qyz .eq. 0.0)  .and.
     &    (qxy .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9904 format ('aptqext DEBUG.  Y cylinder, z and x extrema.')
cbugcbug      write ( 3, 9904)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qz, qx, qzx, qzz, qxx, tol,
     &                np, ez, ex, adir, nerra)

          do 130 n = 1, np
            ey(n) = 0.0
            if (adir(n) .eq. 'ucon') adir(n) = 'Const z'
            if (adir(n) .eq. 'uint') adir(n) = 'Inter z'
            if (adir(n) .eq. 'umax') adir(n) = 'Max z'
            if (adir(n) .eq. 'umin') adir(n) = 'Min z'
            if (adir(n) .eq. 'vcon') adir(n) = 'Const x'
            if (adir(n) .eq. 'vint') adir(n) = 'Inter x'
            if (adir(n) .eq. 'vmax') adir(n) = 'Max x'
            if (adir(n) .eq. 'vmin') adir(n) = 'Min x'
  130     continue

        go to 210

      endif

c=======================================================================********

c.... Test for a cylinder aligned with the z axis.  Only x and y extrema.
c....   Nz = 0 everywhere.

      if ((qz  .eq. 0.0)  .and.
     &    (qzz .eq. 0.0)  .and.
     &    (qzx .eq. 0.0)  .and.
     &    (qyz .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9905 format ('aptqext DEBUG.  Z cylinder, x and y extrema.')
cbugcbug      write ( 3, 9905)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qx, qy, qxy, qxx, qyy, tol,
     &                np, ex, ey, adir, nerra)

          do 140 n = 1, np
            ez(n) = 0.0
            if (adir(n) .eq. 'ucon') adir(n) = 'Const x'
            if (adir(n) .eq. 'uint') adir(n) = 'Inter x'
            if (adir(n) .eq. 'umax') adir(n) = 'Max x'
            if (adir(n) .eq. 'umin') adir(n) = 'Min x'
            if (adir(n) .eq. 'vcon') adir(n) = 'Const y'
            if (adir(n) .eq. 'vint') adir(n) = 'Inter y'
            if (adir(n) .eq. 'vmax') adir(n) = 'Max y'
            if (adir(n) .eq. 'vmin') adir(n) = 'Min y'
  140     continue

        go to 210

      endif

c=======================================================================********

cbugcbugc***DEBUG begins.
cbugcbug 9701 format (/ 'aptqext DEBUG.  Look for x extrema.')
cbugcbug      write ( 3, 9701)
cbugcbugc***DEBUG ends.

c.... Find any x extrema, where Ny = Nz = 0.  Nx is not zero everywhere.

c.... Find the unit normal vector and proximal point of any plane Ny = 0.

      ax = qxy
      ay = 2.0 * qyy
      az = qyz

      vlensq = ax**2 + ay**2 + az**2

      if (vlensq .eq. 0.0) go to 160  ! No Ny = 0 plane exists.  No extreme x.

      vlen = sqrt (vlensq)

      ax = ax / vlen
      ay = ay / vlen
      az = az / vlen

      bx = -ax * qy / vlen
      by = -ay * qy / vlen
      bz = -az * qy / vlen
cbugcbugc***DEBUG begins.
cbugcbug 9501 format ('  A plane exists where Ny=0.' /
cbugcbug     &        '  Ny=0 unit normal',1p3e20.12 /            
cbugcbug     &        '  Ny=0 prox pt    ',1p3e20.12 )
cbugcbug      write ( 3, 9501) ax, ay, az, bx, by, bz
cbugcbugc***DEBUG ends

c.... Find the unit normal vector and proximal point of any plane Nz = 0.

      cx = qzx
      cy = qyz
      cz = 2.0 * qzz

      vlensq = cx**2 + cy**2 + cz**2

      if (vlensq .eq. 0.0) go to 160  ! No Nz = 0 plane exists.  No extreme x.

      vlen = sqrt (vlensq)

      cx = cx / vlen
      cy = cy / vlen
      cz = cz / vlen

      dx = -cx * qz / vlen
      dy = -cy * qz / vlen
      dz = -cz * qz / vlen
cbugcbugc***DEBUG begins.
cbugcbug 9502 format ('  A plane exists where Nz=0.' /
cbugcbug     &        '  Nz=0 unit normal',1p3e20.12 /            
cbugcbug     &        '  Nz=0 prox pt    ',1p3e20.12 )
cbugcbug      write ( 3, 9502) cx, cy, cz, dx, dy, dz
cbugcbug 9906 format ('aptqext DEBUG.  X extrema may exist.')
cbugcbug      write ( 3, 9906)
cbugcbugc***DEBUG ends.

c.... Find any line of intersection of the two planes Ny = 0, Nz = 0.

      call aptplpl (bx, by, bz, ax, ay, az, dx, dy, dz, cx, cy, cz,
     &              1, tol, eex, eey, eez, ffx, ffy, ffz, ux, uy, uz,
     &              ipar, dpmin, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9601 format ('  Intersection of planes Ny=0, Nz=0:' /
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line unit v ',1p3e20.12 /
cbugcbug     &        '  ipar = ',i3 )            
cbugcbug      write ( 3, 9601) eex, eey, eez, ffx, ffy, ffz, ux, uy, uz, ipar
cbugcbugc***DEBUG ends.

      if ((ipar .eq. 1) .and. (dpmin .eq. 0.0)) then  ! Planes coincide.
cbugcbugc***DEBUG begins.
cbugcbug 9602 format ('Planes Ny=0 and Nz=0 coincide.  Find x extrema.' /
cbugcbug     &        '  Main axis (xyz) ',1p3e20.12)            
cbugcbug      write ( 3, 9602) rotm(3,1), rotm(3,2), rotm(3,3)
cbugcbugc***DEBUG ends.

c....   Surface is a cylinder.  Axis has y and z components, no x component.
c....     Set biggest component of axis vector to zero.

        if (abs (rotm(3,2)) .gt. abs (rotm(3,3))) then

c....   Find x extrema at y = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9401 format ('Find x extrema at y = 0.')
cbugcbug      write ( 3, 9401)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qz, qx, qzx, qzz, qxx, tol,
     &                npp, ez, ex, adir, nerra)

        np = 0
        do 165 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = 0.0
            ez(np)   = ez(n)
            adir(np) = 'Const x'
          endif
          if (adir(n) .eq. 'vint') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = 0.0
            ez(np)   = ez(n)
            adir(np) = 'Inter x'
          endif
          if (adir(n) .eq. 'vmin') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = 0.0
            ez(np)   = ez(n)
            adir(np) = 'Min x'
          endif
          if (adir(n) .eq. 'vmax') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = 0.0
            ez(np)   = ez(n)
            adir(np) = 'Max x'
          endif
  165   continue

        else                          ! Z component of axis is bigger.

c....   Find x extrema at z = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9402 format ('Find x extrema at z = 0.')
cbugcbug      write ( 3, 9402)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qy, qx, qxy, qyy, qxx, tol,
     &                npp, ey, ex, adir, nerra)

        np = 0
        do 168 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Const x'
          endif
          if (adir(n) .eq. 'vint') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Inter x'
          endif
          if (adir(n) .eq. 'vmin') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Min x'
          endif
          if (adir(n) .eq. 'vmax') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Max x'
          endif
  168   continue

        endif                         ! Set biggest component of axis to zero.

        go to 210

      endif                           ! Planes coincide.

      if (ipar .ne. 0) go to 160      ! Planes are parallel.  No extreme x.

c....   Find any intersections of the line through points (eex, eey, eez) and
c....   (ffx, ffy, ffz) with unit vector direction (ux, uy, uz) and the
c....   quadric surface.  Such intersections will be x extrema.

c.... Find the coefficients of the equation for the distance to intersection.

        aa = ux * (qxx * ux + qxy * uy) +
     &       uy * (qyy * uy + qyz * uz) +
     &       uz * (qzz * uz + qzx * ux)

        bb = ux * (qx + 2.0 * qxx * eex + qxy * eey + qzx * eez) +
     &       uy * (qy + 2.0 * qyy * eey + qyz * eez + qxy * eex) +
     &       uz * (qz + 2.0 * qzz * eez + qzx * eex + qyz * eey)

        cc = qc +
     &       eex * (qx + qxx * eex + qxy * eey) +
     &       eey * (qy + qyy * eey + qyz * eez) +
     &       eez * (qz + qzz * eez + qzx * eex)


c.... Find the roots of the equation for the distance to intersection.

      call aptqrts (0, aa, bb, cc, qq, tol,
     &              nroots, dint1, dint2, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9301 format ('  Roots = ',i3,5x,1p2e22.14)
cbugcbug      write ( 3, 9301) nroots, dint1, dint2
cbugcbugc***DEBUG ends.

      if ((nroots .eq. 2) .and. (dint2 .lt. dint1)) then
        dintx = dint1
        dint1 = dint2
        dint2 = dintx
      endif

c.... Find the minimum and maximum x points.

      derr = tol**2 * (eex**2 + eey**2 + eez**2)

c-----------------------------------------------------------------------********

      if (nroots .ge. 1) then         ! Test first x root.
        if (dint1**2 .lt. derr) dint1 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint1, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        xnorm = qx + 2.0 * qxx * ex(np) + qxy * ey(np) + qzx * ez(np)
        errxn = tol * (abs (qx) + abs (2.0 * qxx * ex(np)) +
     &                 abs (qxy * ey(np)) + abs ( qzx * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug 9503 format ('  Extreme point   ',1p3e20.12 /             
cbugcbug     &        '  xnorm, error    ',1p2e20.12 )
cbugcbug      write ( 3, 9503) ex(np), ey(np), ez(np), xnorm, errxn
cbugcbugc***DEBUG ends.
        if (abs (xnorm) .le. errxn) then       ! Reject root.
cbugcbugc***DEBUG begins.
cbugcbug 9302 format ('  Rejected point.')
cbugcbug      write ( 3, 9302)
cbugcbugc***DEBUG ends.
          np = np - 1
          go to 150
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point x'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle x'
        else                          ! Quadric surface is real.

          adir(np) = 'Const x'

          do 112 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 112
            if (ccx(nn) .lt. ex(np)) then
               if (adir(np) .eq. 'Min x') then
                 adir(np) = 'Saddle x'
               else
                 adir(np) = 'Max x'
               endif
            elseif (ccx(nn) .gt. ex(np)) then
               if (adir(np) .eq. 'Max x') then
                 adir(np) = 'Saddle x'
               else
                 adir(np) = 'Min x'
               endif
            endif
  112     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9801 format ('aptqext DEBUG  1.  np=  ',i2,' ex=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9801) np, ex(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested first x root.

c-----------------------------------------------------------------------********

  150 if (nroots .eq. 2) then         ! Test second x root.
        if (dint2**2 .lt. derr) dint2 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint2, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        xnorm = qx + 2.0 * qxx * ex(np) + qxy * ey(np) + qzx * ez(np)
        errxn = tol * (abs (qx) + abs (2.0 * qxx * ex(np)) +
     &                 abs (qxy * ey(np)) + abs ( qzx * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9503) ex(np), ey(np), ez(np), xnorm, errxn
cbugcbugc***DEBUG ends.
        if (abs (xnorm) .le. errxn) then       ! Reject root.
          np = np - 1
          go to 160
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point x'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle x'
        else                          ! Quadric surface is real.

          adir(np) = 'Const x'

          do 122 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 122
            if (ccx(nn) .lt. ex(np)) then
               if (adir(np) .eq. 'Min x') then
                 adir(np) = 'Saddle x'
               else
                 adir(np) = 'Max x'
               endif
            elseif (ccx(nn) .gt. ex(np)) then
               if (adir(np) .eq. 'Max x') then
                 adir(np) = 'Saddle x'
               else
                 adir(np) = 'Min x'
               endif
            endif
  122     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9802 format ('aptqext DEBUG  2.  np=  ',i2,' ex=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9802) np, ex(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested second x root.

c=======================================================================********

c.... Find any y extrema, where Nz = Nx = 0.  Ny is not zero everywhere.

c.... Find the unit normal vector and proximal point of any plane Nz = 0.

  160 ax = qzx
      ay = qyz
      az = 2.0 * qzz
cbugcbugc***DEBUG begins.
cbugcbug 9702 format (/ 'aptqext DEBUG.  Look for y extrema.')
cbugcbug      write ( 3, 9702)
cbugcbugc***DEBUG ends.

      vlensq = ax**2 + ay**2 + az**2

      if (vlensq .eq. 0.0) go to 180  ! No Nz = 0 plane exists.  No extreme y.

      vlen = sqrt (vlensq)

      ax = ax / vlen
      ay = ay / vlen
      az = az / vlen

      bx = -ax * qz / vlen
      by = -ay * qz / vlen
      bz = -az * qz / vlen
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9502) ax, ay, az, bx, by, bz
cbugcbugc***DEBUG ends

c.... Find the unit normal vector and proximal point of any plane Nx = 0.

      cx = 2.0 * qxx
      cy = qxy
      cz = qzx

      vlensq = cx**2 + cy**2 + cz**2

      if (vlensq .eq. 0.0) go to 180  ! No Nx = 0 plane exists.  No extreme y.

      vlen = sqrt (vlensq)

      cx = cx / vlen
      cy = cy / vlen
      cz = cz / vlen

      dx = -cx * qx / vlen
      dy = -cy * qx / vlen
      dz = -cz * qx / vlen
cbugcbugc***DEBUG begins.
cbugcbug 9500 format ('  A plane exists where Nx=0.' /
cbugcbug     &        '  Nx=0 unit normal',1p3e20.12 /            
cbugcbug     &        '  Nx=0 prox pt    ',1p3e20.12 )
cbugcbug      write ( 3, 9500) cx, cy, cz, dx, dy, dz
cbugcbug 9907 format ('aptqext DEBUG.  Y extrema may exist.')
cbugcbug      write ( 3, 9907)
cbugcbugc***DEBUG ends.

c.... Find any line of intersection of the two planes Nz = 0, Nx = 0.

      call aptplpl (bx, by, bz, ax, ay, az, dx, dy, dz, cx, cy, cz,
     &              1, tol, eex, eey, eez, ffx, ffy, ffz, ux, uy, uz,
     &              ipar, dpmin, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9603 format ('  Intersection of planes Nz=0, Nx=0:' /
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line unit v ',1p3e20.12 /
cbugcbug     &        '  ipar = ',i3 )            
cbugcbug      write ( 3, 9603) eex, eey, eez, ffx, ffy, ffz, ux, uy, uz, ipar
cbugcbugc***DEBUG ends.

      if ((ipar .eq. 1) .and. (dpmin .eq. 0.0)) then  ! Planes coincide.
cbugcbugc***DEBUG begins.
cbugcbug 9604 format ('Planes Nz=0 and Nx=0 coincide.  Find y extrema.' /
cbugcbug     &        '  Main axis (xyz) ',1p3e20.12)            
cbugcbug      write ( 3, 9604) rotm(3,1), rotm(3,2), rotm(3,3)
cbugcbugc***DEBUG ends.

c....   Surface is a cylinder.  Axis has z and x components, no y component.
c....     Set biggest component of axis vector to zero.

        if (abs (rotm(3,3)) .gt. abs (rotm(3,1))) then

c....   Find y extrema at z = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9403 format ('Find y extrema at z = 0.')
cbugcbug      write ( 3, 9403)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qx, qy, qxy, qxx, qyy, tol,
     &                npp, ex, ey, adir, nerra)

        np = 0
        do 175 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Const y'
          endif
          if (adir(n) .eq. 'vint') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Inter y'
          endif
          if (adir(n) .eq. 'vmin') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Min y'
          endif
          if (adir(n) .eq. 'vmax') then
            np       = np + 1
            ex(np)   = ex(n)
            ey(np)   = ey(n)
            ez(np)   = 0.0
            adir(np) = 'Max y'
          endif
  175   continue

        else                          ! X component of axis is bigger.

c....   Find y extrema at x = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9404 format ('Find y extrema at x = 0.')
cbugcbug      write ( 3, 9404)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qz, qy, qyz, qzz, qyy, tol,
     &                npp, ez, ey, adir, nerra)

        np = 0
        do 178 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np       = np + 1
            ex(np)   = 0.0
            ey(np)   = ey(n)
            ez(np)   = ez(n)
            adir(np) = 'Const y'
          endif
          if (adir(n) .eq. 'vint') then
            np       = np + 1
            ex(np)   = 0.0
            ey(np)   = ey(n)
            ez(np)   = ez(n)
            adir(np) = 'Inter y'
          endif
          if (adir(n) .eq. 'vmin') then
            np       = np + 1
            ex(np)   = 0.0
            ey(np)   = ey(n)
            ez(np)   = ez(n)
            adir(np) = 'Min y'
          endif
          if (adir(n) .eq. 'vmax') then
            np       = np + 1
            ex(np)   = 0.0
            ey(np)   = ey(n)
            ez(np)   = ez(n)
            adir(np) = 'Max y'
          endif
  178   continue

        endif                         ! Set biggest component of axis to zero.

        go to 210

      endif                           ! Planes coincide.

      if (ipar .ne. 0) go to 180      ! Planes are parallel.  No extreme y.

c....   Find any intersections of the line through points (eex, eey, eez) and
c....   (ffx, ffy, ffz) with unit direction vector (ux, uy, uz) and the
c....   quadric surface.  Such intersections will be y extrema.

c.... Find the coefficients of the equation for the distance to intersection.

        aa = uy * (qyy * uy + qyz * uz) +
     &       uz * (qzz * uz + qzx * ux) +
     &       ux * (qxx * ux + qxy * uy)

        bb = uy * (qy + 2.0 * qyy * eey + qyz * eez + qxy * eex) +
     &       uz * (qz + 2.0 * qzz * eez + qzx * eex + qyz * eey) +
     &       ux * (qx + 2.0 * qxx * eex + qxy * eey + qzx * eez)

        cc = qc +
     &       eey * (qy + qyy * eey + qyz * eez) +
     &       eez * (qz + qzz * eez + qzx * eex) +
     &       eex * (qx + qxx * eex + qxy * eey)


c.... Find the roots of the equation for the distance to intersection.

      call aptqrts (0, aa, bb, cc, qq, tol,
     &              nroots, dint1, dint2, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9301) nroots, dint1, dint2
cbugcbugc***DEBUG ends.

c.... Find the minimum and maximum y points.

      derr = tol**2 * (eex**2 + eey**2 + eez**2)

c-----------------------------------------------------------------------********

      if (nroots .ge. 1) then         ! Test first y root.
        if (dint1**2 .lt. derr) dint1 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint1, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        ynorm = qy + qxy * ex(np) + 2.0 * qyy * ey(np) + qyz * ez(np)
        erryn = tol * (abs (qy) + abs (qxy * ex(np)) +
     &                 abs (2.0 * qyy * ey(np)) + abs (qyz * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug 9504 format ('  Extreme point   ',1p3e20.12 /             
cbugcbug     &        '  ynorm, error    ',1p2e20.12 )
cbugcbug      write ( 3, 9504) ex(np), ey(np), ez(np), ynorm, erryn
cbugcbugc***DEBUG ends.
        if (abs (ynorm) .le. erryn) then       ! Reject root.
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9302)
cbugcbugc***DEBUG ends.
          np = np - 1
          go to 170
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point y'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle y'
        else                          ! Quadric surface is real.

          adir(np) = 'Const y'

          do 132 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 132
            if (ccy(nn) .lt. ey(np)) then
               if (adir(np) .eq. 'Min y') then
                 adir(np) = 'Saddle y'
               else
                 adir(np) = 'Max y'
               endif
            elseif (ccy(nn) .gt. ey(np)) then
               if (adir(np) .eq. 'Max y') then
                 adir(np) = 'Saddle y'
               else
                 adir(np) = 'Min y'
               endif
            endif
  132     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9803 format ('aptqext DEBUG  3.  np=  ',i2,' ey=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9803) np, ey(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested first y root.

c-----------------------------------------------------------------------********

  170 if (nroots .eq. 2) then         ! Test second y root.
        if (dint2**2 .lt. derr) dint2 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint2, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        ynorm = qy + qxy * ex(np) + 2.0 * qyy * ey(np) + qyz * ez(np)
        erryn = tol * (abs (qy) + abs (qxy * ex(np)) +
     &                 abs (2.0 * qyy * ey(np)) + abs (qyz * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9504) ex(np), ey(np), ez(np), ynorm, erryn
cbugcbugc***DEBUG ends.
        if (abs (ynorm) .le. erryn) then       ! Reject root.
          np = np - 1
          go to 180
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point y'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle y'
        else                          ! Quadric surface is real.

          adir(np) = 'Const y'

          do 142 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 142
            if (ccy(nn) .lt. ey(np)) then
               if (adir(np) .eq. 'Min y') then
                 adir(np) = 'Saddle y'
               else
                 adir(np) = 'Max y'
               endif
            elseif (ccy(nn) .gt. ey(np)) then
               if (adir(np) .eq. 'Max y') then
                 adir(np) = 'Saddle y'
               else
                 adir(np) = 'Min y'
               endif
            endif
  142     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9804 format ('aptqext DEBUG  4.  np=  ',i2,' ey=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9804) np, ey(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested second y root.

c=======================================================================********

c.... Find any z extrema, where Nx = Ny = 0.  Nz is not zero everywhere.

c.... Find the unit normal vector and proximal point of any plane Nx = 0.

  180 ax = 2.0 * qxx
      ay = qxy
      az = qzx
cbugcbugc***DEBUG begins.
cbugcbug 9703 format (/ 'aptqext DEBUG.  Look for z extrema.')
cbugcbug      write ( 3, 9703)
cbugcbugc***DEBUG ends.

      vlensq = ax**2 + ay**2 + az**2

      if (vlensq .eq. 0.0) go to 210  ! No Nx = 0 plane exists.  No extreme z.

      vlen = sqrt (vlensq)

      ax = ax / vlen
      ay = ay / vlen
      az = az / vlen

      bx = -ax * qx / vlen
      by = -ay * qx / vlen
      bz = -az * qx / vlen
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9500) ax, ay, az, bx, by, bz
cbugcbugc***DEBUG ends

c.... Find the unit normal vector and proximal point of any plane Ny = 0.

      cx = qxy
      cy = 2.0 * qyy
      cz = qyz

      vlensq = cx**2 + cy**2 + cz**2

      if (vlensq .eq. 0.0) go to 210  ! No Ny = 0 plane exists.  No extreme z.

      vlen = sqrt (vlensq)

      cx = cx / vlen
      cy = cy / vlen
      cz = cz / vlen

      dx = -cx * qy / vlen
      dy = -cy * qy / vlen
      dz = -cz * qy / vlen
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9501) cx, cy, cz, dx, dy, dz
cbugcbug 9908 format ('aptqext DEBUG.  Z extrema may exist.')
cbugcbug      write ( 3, 9908)
cbugcbugc***DEBUG ends.

c.... Find any line of intersection of the two planes Nx = 0, Ny = 0.

      call aptplpl (bx, by, bz, ax, ay, az, dx, dy, dz, cx, cy, cz,
     &              1, tol, eex, eey, eez, ffx, ffy, ffz, ux, uy, uz,
     &              ipar, dpmin, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9605 format ('  Intersection of planes Nx=0, Ny=0:' /
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line pt xyz ',1p3e20.12 /            
cbugcbug     &        '  Int line unit v ',1p3e20.12 /
cbugcbug     &        '  ipar = ',i3 )            
cbugcbug      write ( 3, 9605) eex, eey, eez, ffx, ffy, ffz, ux, uy, uz, ipar
cbugcbugc***DEBUG ends.

      if ((ipar .eq. 1) .and. (dpmin .eq. 0.0)) then  ! Planes coincide.
cbugcbugc***DEBUG begins.
cbugcbug 9606 format ('Planes Nx=0 and Ny=0 coincide.  Find z extrema.' /
cbugcbug     &        '  Main axis (xyz) ',1p3e20.12)            
cbugcbug      write ( 3, 9606) rotm(3,1), rotm(3,2), rotm(3,3)
cbugcbugc***DEBUG ends.

c....   Surface is a cylinder.  Axis has x and y components, no z component.
c....     Set biggest component of axis vector to zero.

        if (abs (rotm(3,1)) .gt. abs (rotm(3,2))) then

c....   Find z extrema at x = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9405 format ('Find z extrema at x = 0.')
cbugcbug      write ( 3, 9405)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qy, qz, qyz, qyy, qzz, tol,
     &                npp, ey, ez, adir, nerra)

        np = 0
        do 185 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np      = np + 1
            ex(np)  = 0.0
            ey(np)  = ey(n)
            ez(np)  = ez(n)
            adir(np) = 'Const z'
          endif
          if (adir(n) .eq. 'vint') then
            np      = np + 1
            ex(np)  = 0.0
            ey(np)  = ey(n)
            ez(np)  = ez(n)
            adir(np) = 'Inter z'
          endif
          if (adir(n) .eq. 'vmin') then
            np      = np + 1
            ex(np)  = 0.0
            ey(np)  = ey(n)
            ez(np)  = ez(n)
            adir(np) = 'Min z'
          endif
          if (adir(n) .eq. 'vmax') then
            np      = np + 1
            ex(np)  = 0.0
            ey(np)  = ey(n)
            ez(np)  = ez(n)
            adir(np) = 'Max z'
          endif
  185   continue

        else                          ! Y component of axis is bigger.

c....   Find z extrema at y = 0.
cbugcbugc***DEBUG begins.
cbugcbug 9406 format ('Find z extrema at y = 0.')
cbugcbug      write ( 3, 9406)
cbugcbugc***DEBUG ends.

        call aptqexc (qc, qx, qz, qzx, qxx, qzz, tol,
     &                npp, ex, ez, adir, nerra)

        np = 0
        do 188 n = 1, npp
          if (adir(n) .eq. 'vcon') then
            np      = np + 1
            ex(np)  = ex(n)
            ey(np)  = 0.0
            ez(np)  = ez(n)
            adir(np) = 'Const z'
          endif
          if (adir(n) .eq. 'vint') then
            np      = np + 1
            ex(np)  = ex(n)
            ey(np)  = 0.0
            ez(np)  = ez(n)
            adir(np) = 'Inter z'
          endif
          if (adir(n) .eq. 'vmin') then
            np      = np + 1
            ex(np)  = ex(n)
            ey(np)  = 0.0
            ez(np)  = ez(n)
            adir(np) = 'Min z'
          endif
          if (adir(n) .eq. 'vmax') then
            np      = np + 1
            ex(np)  = ex(n)
            ey(np)  = 0.0
            ez(np)  = ez(n)
            adir(np) = 'Max z'
          endif
  188   continue

        endif                         ! Set biggest component of axis to zero.

        go to 210

      endif                           ! Planes coincide.

      if (ipar .ne. 0) go to 210      ! Planes are parallel.  No extreme z.

c....   Find any intersections of the line through points (eex, eey, eez) and
c....   (ffx, ffy, ffz) with unit direction vector (ux, uy, uz) and the
c....   quadric surface.  Such intersections will be z extrema.

c.... Find the coefficients of the equation for the distance to intersection.

        aa = uz * (qzz * uz + qzx * ux) +
     &       ux * (qxx * ux + qxy * uy) +
     &       uy * (qyy * uy + qyz * uz)

        bb = uz * (qz + 2.0 * qzz * eez + qzx * eex + qyz * eey) +
     &       ux * (qx + 2.0 * qxx * eex + qxy * eey + qzx * eez) +
     &       uy * (qy + 2.0 * qyy * eey + qyz * eez + qxy * eex)

        cc = qc +
     &       eez * (qz + qzz * eez + qzx * eex) +
     &       eex * (qx + qxx * eex + qxy * eey) +
     &       eey * (qy + qyy * eey + qyz * eez)


c.... Find the roots of the equation for the distance to intersection.

      call aptqrts (0, aa, bb, cc, qq, tol,
     &              nroots, dint1, dint2, itrun, nerra)
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9301) nroots, dint1, dint2
cbugcbugc***DEBUG ends.

c.... Find the minimum and maximum z points.

      derr = tol**2 * (eex**2 + eey**2 + eez**2)

c-----------------------------------------------------------------------********

      if (nroots .ge. 1) then         ! Test first z root.
        if (dint1**2 .lt. derr) dint1 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint1, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        znorm = qz + qzx * ex(np) + qyz * ey(np) + 2.0 * qzz * ez(np)
        errzn = tol * (abs (qz) + abs (qzx * ex(np)) +
     &                 abs (qyz * ey(np)) + abs (2.0 * qzz * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug 9505 format ('  Extreme point   ',1p3e20.12 /             
cbugcbug     &        '  znorm, error    ',1p2e20.12 )
cbugcbug      write ( 3, 9505) ex(np), ey(np), ez(np), znorm, errzn
cbugcbugc***DEBUG ends.
        if (abs (znorm) .le. errzn) then       ! Reject root.
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9302)
cbugcbugc***DEBUG ends.
          np = np - 1
          go to 190
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point z'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle z'
        else                          ! Quadric surface is real.

          adir(np) = 'Const z'

          do 152 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 152
            if (ccz(nn) .lt. ez(np)) then
               if (adir(np) .eq. 'Min z') then
                 adir(np) = 'Saddle z'
               else
                 adir(np) = 'Max z'
               endif
            elseif (ccz(nn) .gt. ez(np)) then
               if (adir(np) .eq. 'Max z') then
                 adir(np) = 'Saddle z'
               else
                 adir(np) = 'Min z'
               endif
            endif
  152     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9805 format ('aptqext DEBUG  5.  np=  ',i2,' ez=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9805) np, ez(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested first z root.

c-----------------------------------------------------------------------********

  190 if (nroots .eq. 2) then         ! Test second z root.
        if (dint2**2 .lt. derr) dint2 = 0.0
        np = np + 1

        call aptvadd (eex, eey, eez, 1.0, dint2, ux, uy, uz, 1, tol,
     &                ex(np), ey(np), ez(np), vlen, nerra)

        znorm = qz + qzx * ex(np) + qyz * ey(np) + 2.0 * qzz * ez(np)
        errzn = tol * (abs (qz) + abs (qzx * ex(np)) +
     &                 abs (qyz * ey(np)) + abs (2.0 * qzz * ez(np)))
cbugcbugc***DEBUG begins.
cbugcbug      write ( 3, 9505) ex(np), ey(np), ez(np), znorm, errzn
cbugcbugc***DEBUG ends.
        if (abs (znorm) .le. errzn) then       ! Reject root.
          np = np - 1
          go to 210
        endif

c....   Find the extreme centers of curvature at the extreme point.

        call aptcris (ex(np), ey(np), ez(np), qc, qx, qy, qz,
     &                qxy, qyz, qzx, qxx, qyy, qzz, tolx,
     &                bcx, bcy, bcz, nrad, rcurve, ccx, ccy, ccz,
     &                ucx, ucy, ucz, nerra)

        if ((nerra .ne. 0) .or. (nrad .eq. 0)) then  ! Quadric surface bad.
          adir(np) = 'Point z'
        elseif ((ntype .eq. 8) .or. (ntype .eq. 13)) then
          adir(np) = 'Saddle z'
        else                          ! Quadric surface is real.

          adir(np) = 'Const z'

          do 162 nn = 1, nrad         ! Loop over finite radii of curvature.
            if (rcurve(nn) .eq. 1.e99) go to 162
            if (ccz(nn) .lt. ez(np)) then
               if (adir(np) .eq. 'Min z') then
                 adir(np) = 'Saddle z'
               else
                 adir(np) = 'Max z'
               endif
            elseif (ccz(nn) .gt. ez(np)) then
               if (adir(np) .eq. 'Max z') then
                 adir(np) = 'Saddle z'
               else
                 adir(np) = 'Min z'
               endif
            endif
  162     continue                    ! End of loop over radii of curvature.
        endif                         ! Tested nerra, nrad.

cbugcbugc***DEBUG begins.
cbugcbug 9806 format ('aptqext DEBUG  6.  np=  ',i2,' ez=',1pe20.12,' adir=',a8)
cbugcbug      write ( 3, 9806) np, ez(np), adir(np)
cbugcbugc***DEBUG ends.

      endif                           ! Tested second z root.

c=======================================================================********

  210 continue
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptqext results:  ntype=',i3,' np=',i2,' nerr=',i2)
cbug 9911 format ('Center    xyz=',1p3e22.14)
cbug 9912 format (a8,'  xyz=',1p3e22.14)
cbug      write ( 3, 9910) ntype, np, nerr
cbug      write ( 3, 9911) cenx, ceny, cenz
cbug      do 991 n = 1, np
cbug        write ( 3, 9912) adir(n), ex(n), ey(n), ez(n)
cbug  991 continue
cbugc***DEBUG ends.

      return

c.... End of subroutine aptqext.      (+1 line.)
      end

UCRL-WEB-209832