subroutine aptesis (ax, ay, az, np, tol, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTESIS
c
c     call aptesis (ax, ay, az, np, tol, qc, qx, qy, qz,
c    &              qxy, qyz, qzx, qxx, qyy, qzz, nerr)
c
c     Version:  aptesis  Updated    1990 December 13 10:30.
c               aptesis  Originated 1990 December 13 10:30.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  For each of np sets of input data, to find the equation of the
c               ellipsoid centered at the origin, with semiaxes of ax, ay
c               and az on the x, y and z axes, respectively.
c
c               The equation for the ellipsoid is as follows:
c
c                 f(u,v,w) = qc + qx * u + qy * v + qz * w +
c                            qxy * u * v + qyz * v * w + qzx * w * u +
c                            qxx * u**2 + qyy * v**2 + qzz * w**2 = 0
c
c               where qx = qy = qz = qxy = qyz = qzx = 0,
c               qc = - ax**2 * ay**2 * az**2, qxx = ay**2 * az**2,
c               qyy = az**2 * ax**2, qzz = ax**2 * ay**2.
c
c               General case:
c
c               If ax, ay and az are all non-zero, the surface will be
c               an ellipsoid with semiaxes ax, ay, az on the x, y and
c               z axes, respectively.
c
c               Indeterminate cases:
c
c               If none, or only one of ax, ay or az is nonzero, as determined
c               by tol, the surface is indeterminate (or a point at the origin),
c               and all coefficients will be zero.
c
c               Other singular cases:
c
c               If ax = ay = az, the surface will be a sphere of radius ax.
c
c               If any two of ax, ay or az are equal, the surface will be an
c               ellipsoid symmetric around the axis of the non-equal semiaxis.
c
c               If just one of the values ax, ay or az is zero, as determined
c               by tol, the surface will be an x, y or z plane through the
c               origin, respectively.
c
c     Note:     Use apttris and aptrois to translate and rotate the ellipsoid
c               into any desired position and orientation.
c               Use aptelis to find the equation of an ellipsoid symmetric
c               around a major axis and passing through two given points.
c
c     Input:    ax, ay, az, np, tol.
c
c     Output:   qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr.
c
c     Glossary:
c
c     ax,ay,az  Input    Semiaxes of the ellipsoid.  Aligned in the directions
c                          of the x, y and z axis, respectively.
c
c     nerr      Output   Indicates an input error, if not zero:
c                          1 if np is not positive.
c
c     np        Input    Size of arrays ax, ay, az, qc, qx, qy, qz,
c                          qxy, qyz, qzx, qxx, qyy, qzz.
c
c     qc        Output   Constant term in the equation of the ellipsoid.
c                          Size np.  qc = -ax**2 * ay**2 * az**2.
c
c     qx,qy,qz  Output   Coefficients of u, v, w in the equation of the
c                          ellipsoid.  Size np.  qx = qy = qz = 0.
c
c     q..       Output   Coefficients of second-order terms in the equation of
c                          the ellipsoid:  qxy, qyz, qzx, qxx, qyy, qzz.
c                          Coefficients of u*v, v*w, w*u, u**2, v**2, w**2,
c                          respectively.  Size np.  qxy = qyz = qzx = 0.
c                          qxx = ay**2 * az**2, qyy = az**2 * ax**2,
c                          qzz = ax**2 * ay**2.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Semiaxis along the x axis.
      dimension ax      (1)
c---- Semiaxis along the y axis.
      dimension ay      (1)
c---- Semiaxis along the z axis.
      dimension az      (1)
c---- Coefficient of term 1.0.
      dimension qc      (1)
c---- Coefficient of term u.
      dimension qx      (1)
c---- Coefficient of term u**2.
      dimension qxx     (1)
c---- Coefficient of term u * v.
      dimension qxy     (1)
c---- Coefficient of term v.
      dimension qy      (1)
c---- Coefficient of term v**2.
      dimension qyy     (1)
c---- Coefficient of term v * w.
      dimension qyz     (1)
c---- Coefficient of term w.
      dimension qz      (1)
c---- Coefficient of term w * u.
      dimension qzx     (1)
c---- Coefficient of term w**2.
      dimension qzz     (1)

c.... Local variables.

c---- Estimated error in ax, ay, az.
      common /laptesis/ aerr
c---- Index in external array.
      common /laptesis/ n

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptesis finding the equation of a ellipsoid.',
cbug     &  '  tol=',1pe13.5)
cbug 9902 format (i3,'  ax,ay,az=',1p3e22.14)
cbug      write ( 3, 9901) tol
cbug      write ( 3, 9902) (n, ax(n), ay(n), az(n), n = 1, np)
cbugc***DEBUG ends.

c.... Initialize.

      nerr = 0

c.... Test for input errors.

      if (np .le. 0) then
        nerr = 1
cbugc***DEBUG begins.
cbug 9903   format (/ "aptesis fatal.  Non-positive np =",i5)
cbug        write ( 3, 9903) np
cbugc***DEBUG ends.
        go to 210
      endif

c.... Find the coefficients of the equation of the ellipsoid.

c---- Loop over data.
      do 120 n = 1, np

        qc(n)  = -ax(n)**2 * ay(n)**2 * az(n)**2
        qx(n)  = 0.0
        qy(n)  = 0.0
        qz(n)  = 0.0
        qxy(n) = 0.0
        qyz(n) = 0.0
        qzx(n) = 0.0
        qxx(n) = ay(n)**2 * az(n)**2
        qyy(n) = az(n)**2 * ax(n)**2
        qzz(n) = ax(n)**2 * ay(n)**2

c---- End of loop over data.
  120 continue

c.... See if the results should be tested for truncation error.

c---- Test for truncation to zero.
      if (tol .gt. 0.0) then

c---- Loop over data.
        do 130 n = 1, np

          aerr = tol * (abs (ax(n)) + abs (ay(n)) + abs (az(n)))

          if (abs (ax(n)) .lt. aerr) then
            qc(n)  = 0.0
            qyy(n) = 0.0
            qzz(n) = 0.0
          endif

          if (abs (ay(n)) .lt. aerr) then
            qc(n)  = 0.0
            qxx(n) = 0.0
            qzz(n) = 0.0
          endif

          if (abs (az(n)) .lt. aerr) then
            qc(n)  = 0.0
            qxx(n) = 0.0
            qyy(n) = 0.0
          endif

c---- End of loop over data.
  130   continue

c---- Tested tol.
      endif
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptesis results:' /
cbug     &  (i3,'  qc=      ',1pe22.14 /
cbug     &  '  qx,qy,qz=   ',1p3e22.14 /
cbug     &  '  qxy,qyz,qzx=',1p3e22.14 /
cbug     &  '  qxx,qyy,qzz=',1p3e22.14))
cbug      write ( 3, 9904) (n, qc(n), qx(n), qy(n), qz(n),
cbug     &  qxy(n), qyz(n), qzx(n), qxx(n), qyy(n), qzz(n),
cbug     &  n = 1, np)
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832