subroutine aptpfit (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
     &                    qc, qx, qy, qz, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTPFIT
c
c     call aptpfit (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
c    &              qc, qx, qy, qz, nerr)
c
c     Version:  aptpfit  Updated    1999 March 3 10:30.
c               aptpfit  Originated 1999 March 3 10:30.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  For the three points a = (ax, ay, az), b = (bx, by, bz),
c               c = (cx, cy, cz), to find the coefficients qc, qx, qy and qz
c               of the equation of the plane through the three points:
c                 F(x,y,z) = qc + qx * x + qy * y + qz * z = 0
c               The distance of the plane from the origin is qc.  The unit
c               vector normal to the plane is q = (qx, qy, qz).
c               Flag nerr indicates an input error, if not zero.
c
c     Input:    ax, ay, az, bx, by, bz, cx, cy, cz, tol.
c
c     Output:   qc, qx, qy, qz, nerr.
c
c     Calls: aptvdot, aptvpln 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z coordinates of point "a".
c
c     bx,by,bz  Input    The x, y, z coordinates of point "b".
c
c     cx,cy,cz  Input    The x, y, z coordinates of point "c".
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if the points are colinear or any are coincident.
c                            If so, -1.e99 is returned for qc, qx, qy, qz.
c
c     qc        Output   The distance of the plane from the origin.
c
c     qx,qy,qz  Output   The x, y, z components of the unit vector normal to the
c                          plane.  A component may be truncated to zero, if less
c                          than the estimated numerical error in its
c                          calculation, based on  tol.
c
c     tol       Input    Numerical tolerance limit for qc, qx, qy, qz.
c                          On computers with 48-bit floating point numbers,
c                          recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpfit finding the plane through the three points:' /
cbug     &  '  ax,ay,az=   ',1p3e22.14 /
cbug     &  '  bx,by,bz=   ',1p3e22.14 /
cbug     &  '  cx,cy,cz=   ',1p3e22.14 )
cbug      write ( 3, 9901) ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.

c.... Initialize.

      nerr = 0

c.... Find the normal vectors.

      call aptvpln (ax, ay, az, bx, by, bz, cx, cy, cz, 1, tol,
     &              qx, qy, qz, vlen, nerra)

c.... See if points are colinear, or if any are coincident.

      if (vlen .eq. 0.0) then
        nerr = 1
        qc   = -1.e99
        qx   = -1.e99
        qy   = -1.e99
        qz   = -1.e99
        go to 210
      endif

c.... Find the dot product of the vector to point "a" and the normal vector.

      call aptvdot (ax, ay, az, qx, qy, qz, 1, tol, qc, nerra)

c.... Find the coefficients of the quadric equation for the plane.

      qc = -qc / vlen
      qx =  qx / vlen
      qy =  qy / vlen
      qz =  qz / vlen

  210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptpfit results.  nerr=',i3 /
cbug     &  '  qc=         ',1pe22.14 /
cbug     &  '  qx,qy,qz=   ',1p3e22.14)
cbug      write ( 3, 9902) nerr, qc, qx, qy, qz
cbugc***DEBUG ends.

      return

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

UCRL-WEB-209832