subroutine aptcnis (aw, ar, bw, br, np, tol, qc, qu, qv, qw,
     &                    quv, qvw, qwu, quu, qvv, qww, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTCNIS
c
c     call aptcnis (aw, ar, bw, br, np, tol, qc, qu, qv, qw,
c    &              quv, qvw, qwu, quu, qvv, qww, nerr)
c
c     Version:  aptcnis  Updated    1990 December 12 14:30.
c               aptcnis  Originated 1990 December 7 11:50.
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
c               a cone symmetric around the w axis, that passes through the
c               two distinct points a = (aw, ar) and b = (bw, br),
c               where ar = sqrt (au**2 + av**2) and br = sqrt (bu**2 + bv**2).
c               Points "a" and "b" are on the same side of the vertex of the
c               cone if ar and br have the same sign, otherwise are on opposite
c               sides.  The coordinates (u,v,w) are along the three major
c               orthogonal axes, and may represent (x,y,z), (y,z,x) or (z,x,y).
c
c               The equation for the conical surface is as follows:
c
c                 f(u,v,w) = qc + qu * u + qv * v + qw * w +
c                            quv * u * v + qvw * v * w + qwu * w * u +
c                            quu * u**2 + qvv * v**2 + qww * w**2 = 0
c
c               where qu = qv = quv = qvw = qwu = 0, quu = qvv, and
c               qww has the opposite sign of quu.
c
c               Degenerate case:
c
c               If points "a" and "b" coincide, within the limits determined
c               by tol, all coefficients will be zero.
c
c               Other special cases:
c
c               If aw = bw, the surface will be a w plane.
c
c               If ar = br, the surface will be a right circular cylinder.
c
c               If ar = br = 0, the surface will be a line on the w axis.
c
c     Note:     Use apttris and aptrois to translate and rotate the cone
c               into any desired position and orientation.
c               Use aptcois to find the equation of a cone on a given axis,
c               with a given vertex angle and vertex axial coordinate.
c               Use aptcyis to find the equation of a cylinder symmetric
c               Around a major axis.
c               Use aptelis to find the equation of an ellipsoid symmetric
c               around a major axis.
c
c     Input:    aw, ar, bw, br, np, tol.
c
c     Output:   qc, qu, qv, qw, quv, qvw, qwu, quu, qvv, qww, nerr.
c
c     Calls: aptvdic, aptvaxc 
c               
c
c     Glossary:
c
c     ar, aw    Input    Coordinates r and w of first point in conical surface.
c                          Size np.
c                          If points "a" and "b" are on the same side of the
c                          vertex of the cone, the signs of ar and br must be
c                          the same.  Otherwise, the signs must be different.
c
c     br, bw    Input    Coordinates r and w of second point in conical surface.
c                          Size np.
c                          If points "a" and "b" are on the same side of the
c                          vertex of the cone, the signs of ar and br must be
c                          the same.  Otherwise, the signs must be different.
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 aw, ar, bw, br, qc, qu, qv, qw,
c                          quv, qvw, qwu, quu, qvv, qww.
c
c     qc        Output   Constant term in quadric cone equation.  Size np.
c                          Zero only if the vertex of the cone is at the origin.
c                          qc = -(aw * br - ar * bw)**2.
c
c     qu,qv,qw  Output   Coefficients of u, v, w in quadric cone equation.
c                          Size np.
c                          qu = qv = 0.
c                          qw = 2.0 * (aw * br - ar * bw) * (br - ar).
c
c     q..       Output   Coefficients of second-order terms in quadric
c                          cone equation:  quv, qvw, qwu, quu, qvv, qww.
c                          Coefficients of u*v, v*w, w*u, u**2, v**2, w**2,
c                          respectively.  Size np.
c                          quv = qvw = qwu = 0.
c                          quu = qvv = (bw - aw)**2.
c                          qww = -(br - ar)**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---- Coordinate r of point "a".
      dimension ar      (1)
c---- Coordinate w of point "a".
      dimension aw      (1)
c---- Coordinate r of point "b".
      dimension br      (1)
c---- Coordinate w of point "b".
      dimension bw      (1)
c---- Coefficient of term 1.0.
      dimension qc      (1)
c---- Coefficient of term u.
      dimension qu      (1)
c---- Coefficient of term u**2.
      dimension quu     (1)
c---- Coefficient of term u * v.
      dimension quv     (1)
c---- Coefficient of term v.
      dimension qv      (1)
c---- Coefficient of term v**2.
      dimension qvv     (1)
c---- Coefficient of term v * w.
      dimension qvw     (1)
c---- Coefficient of term w.
      dimension qw      (1)
c---- Coefficient of term w * u.
      dimension qwu     (1)
c---- Coefficient of term w**2.
      dimension qww     (1)

c.... Local variables.

c---- Distance from point "a" to point "b".
      common /laptcnis/ abd     (64)
c---- Cross product of vectors "a" and "b".
      common /laptcnis/ abm     (64)
c---- Difference br - ar.
      common /laptcnis/ abr     (64)
c---- Difference bw - aw.
      common /laptcnis/ abw     (64)
c---- Index in local array.
      common /laptcnis/ n
c---- First index of subset of data.
      common /laptcnis/ n1
c---- Last index of subset of data.
      common /laptcnis/ n2
c---- Index in external array.
      common /laptcnis/ nn
c---- Size of current subset of data.
      common /laptcnis/ ns

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcnis finding the quadric coefficients for a cone.',
cbug     &  '  tol=',1pe13.5 /
cbug     &  '  Cone passes through the points:')
cbug 9902 format (
cbug     &  i3,'  aw,ar=',1p2e22.14 /
cbug     &  '     bw,br=',1p2e22.14)
cbug      write ( 3, 9901) tol
cbug      write ( 3, 9902) (n, aw(n), ar(n), bw(n), br(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 (/ "aptcnis fatal.  Non-positive np =",i5)
cbug        write ( 3, 9903) np
cbugc***DEBUG ends.
        go to 210
      endif

c.... Set up the indices of the first subset of data.

      n1 = 1
      n2 = min (np, 64)

  110 ns = n2 - n1 + 1

c.... Find the 2-D vector from point "a" to point "b".

      call aptvdic (aw(n1), ar(n1), bw(n1), br(n1), ns, tol,
     &              abw, abr, abd, nerr)

c.... Find the cross product of the 2-D vectors "a" and "b".

      call aptvaxc (aw(n1), ar(n1), bw(n1), br(n1), ns, tol, abm, nerr)

c.... Find the coefficients of the equation of the cone.

c---- Loop over subset of data.
      do 120 n = 1, ns

        nn = n + n1 - 1

        qc(nn)  = -abm(n)**2
        qu(nn)  = 0.0
        qv(nn)  = 0.0
        qw(nn)  = 2.0 * abm(n) * abr(n)
        quv(nn) = 0.0
        qvw(nn) = 0.0
        qwu(nn) = 0.0
        quu(nn) = abw(n)**2
        qvv(nn) = quu(nn)
        qww(nn) = -abr(n)**2

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

c.... See if all data subsets are done.

c---- Do another subset of data.
      if (n2 .lt. np) then
        n1 = n2 + 1
        n2 = min (np, n1 + 63)
        go to 110
      endif
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptcnis results:' /
cbug     &  (i3,'  qc=      ',1pe22.14 /
cbug     &  '  qu,qv,qw=   ',1p3e22.14 /
cbug     &  '  quv,qvw,qwu=',1p3e22.14 /
cbug     &  '  quu,qvv,qww=',1p3e22.14))
cbug      write ( 3, 9904) (n, qc(n), qu(n), qv(n), qw(n),
cbug     &  quv(n), qvw(n), qwu(n), quu(n), qvv(n), qww(n),
cbug     &  n = 1, np)
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832