subroutine aptcois (kang, theta, aw, np, tol, qc, qu, qv, qw,
     &                    quv, qvw, qwu, quu, qvv, qww, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTCOIS
c
c     call aptcois (kang, theta, aw, np, tol, qc, qu, qv, qw,
c    &              quv, qvw, qwu, quu, qvv, qww, nerr)
c
c     Version:  aptcois  Updated    1990 December 12 14:30.
c               aptcois  Originated 1990 December 7 11:40.
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, with the vertex at w = aw,
c               and the angle theta between the w axis and the conical surface.
c               The coordinates (u,v,w) are along the three major orthogonal
c               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               Special cases:
c
c               The result will be a line if theta is zero plus or minus
c               any integer multiple of 180 degrees (pi radians).
c
c               The result will be a plane if theta is 90 degrees
c               (pi/2 radians) plus or minus any integer multiple of
c               180 degrees (pi radians).
c
c     Note:     Use apttris and aptrois to translate and rotate the cone
c               into any desired position and orientation.
c               Use aptcnis to find the equation of a cone on a given axis,
c               and passing through two given points.
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:    kang, theta, aw, np, tol.
c
c     Output:   qc, qu, qv, qw, quv, qvw, qwu, quu, qvv, qww, nerr.
c
c     Glossary:
c
c     aw        Input    Coordinate w of the vertex of the cone on the w axis.
c                          Size np.
c
c     kang      Input    Indicates theta units are degrees (0) or radians (1).
c
c     nerr      Output   Indicates an input error, if not zero:
c                          1 if np is not positive.
c                          2 if kang is not 0 or 1.
c
c     np        Input    Size of arrays theta, aw, 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 cone passes through the origin.
c                          qc = -aw**2 * (sin (theta))**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 * (sin (theta))**2.
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 = (cos (theta))**2.
c                          qww = -(sin (theta))**2.
c
c     theta     Input    Angle from w axis to conical surface.  Units are
c                          degrees (kang = 0) or radians (kang = 1).  Size np.
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.  Size np.

c---- Coordinate w of vertex of cone.
      dimension aw      (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---- Angle from w axis to cone.
      dimension theta   (1)

c.... Local variables.

c---- Vertex angle in radians.
      common /laptcois/ angle
c---- Square of cos (theta).
      common /laptcois/ cthsq   (64)
c---- Index in local array.
      common /laptcois/ n
c---- First index of subset of data.
      common /laptcois/ n1
c---- Last index of subset of data.
      common /laptcois/ n2
c---- Index in external array.
      common /laptcois/ nn
c---- Size of current subset of data.
      common /laptcois/ ns

c---- Mathematical constant, pi.
      common /laptcois/ pi

c---- Square of sin (theta).
      common /laptcois/ sthsq   (64)

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcois finding the quadric coefficients for a cone.',
cbug     &  '  tol=',1pe13.5,'  kang=',i2)
cbug 9902 format (i3,'  theta=',1pe22.14,'  aw=',1pe22.14)
cbug      write ( 3, 9901) tol, kang
cbug      write ( 3, 9902) (n, theta(n), aw(n), n = 1, np)
cbugc***DEBUG ends.

c.... Initialize.

c---- Mathematical constant, pi.
      pi = 3.14159265358979323

      nerr = 0

c.... Test for input errors.

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

      if ((kang .lt. 0) .or. (kang .gt. 1)) then
        nerr = 2
cbugc***DEBUG begins.
cbug        write ( 3, '(/ "aptcois fatal.  Bad kang=",i3)') kang
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 vertex angle functions.

c---- Angles are in degrees.
      if (kang .eq. 0) then

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

          nn       = n + n1 - 1
          angle    = theta(nn) * pi / 180.0
          cthsq(n) = (cos (angle))**2
          sthsq(n) = (sin (angle))**2

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

c---- Angles are in radians.
      else

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

          nn       = n + n1 - 1
          angle    = theta(nn)
          cthsq(n) = (cos (angle))**2
          sthsq(n) = (sin (angle))**2

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

      endif

c.... See if truncation tests should be made.

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

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

          if (abs (sthsq(n)) .lt. tol) then
            sthsq(n) = 0.0
          endif

          if (abs (cthsq(n)) .lt. tol) then
            cthsq(n) = 0.0
          endif

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

c---- Tested tol.
      endif

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

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

        nn = n + n1 - 1

        qc(nn)  = -aw(nn)**2 * sthsq(n)
        qu(nn)  = 0.0
        qv(nn)  = 0.0
        qw(nn)  = 2.0 * aw(nn) * sthsq(n)
        quv(nn) = 0.0
        qvw(nn) = 0.0
        qwu(nn) = 0.0
        quu(nn) = cthsq(n)
        qvv(nn) = quu(nn)
        qww(nn) = -sthsq(n)

c---- End of loop over subset of data.
  150 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 9905 format (/ 'aptcois 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, 9905) (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 aptcois.      (+1 line.)
      end

UCRL-WEB-209832