subroutine aptparx (px, py, pz, vx, vy, vz, ax, ay, az,
     &                    qc, qx, qy, qz, qxy, qyz, qzx,
     &                    qxx, qyy, qzz, tol,
     &                    nprox, atype, ptx, pty, ptz,
     &                    vtx, vty, vtz, vtlen, dprox, tprox,
     &                    dside, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTPARX
c
c     call aptparx (px, py, pz, vx, vy, vz, ax, ay, az,
c    &              qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol,
c    &              nprox, atype, ptx, pty, ptz,
c    &              vtx, vty, vtz, vtlen, dprox, tprox,
c    &              dside, nerr)
c
c     Version:  aptparx  Updated    1997 May 6 15:25.
c               aptparx  Originated 1997 May 6 15:25.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  For the particle with initial coordinates p = (px, py, pz),
c               initial velocity v = (vx, vy, vz) and constant acceleration
c               a = (ax, ay, az), and for the quadric surface with the implicit
c               equation F(x,y,z) = qc + qx*x + qy*y + qz*z + qxy*x*y +
c               qyz*y*z + qzx*z*x + qxx*x**2 + qyy*y**2 + qzz*z**2 = 0,
c               to find the number nprox (0, 1, 2 or 3) of extrema of F on
c               the particle path,
c               the extremum type ("minimum", "min/max" or "maximum"),
c               the coordinates pt(n) = (ptx(n), pty(n), ptz(n)),
c               the velocities vt(n) = (vtx(n), vty(n), vtz(n)),
c               the velocity magnitudes vtlen(n),
c               the path lengths to extremum dprox(n),
c               the times to extremum tprox(n), and
c               the value dside(n) of F(x,y,z), all for n = 1, nprox.
c               Flag nerr indicates any input error.
c
c               The x, y and z components of the particle velocity and
c               coordinates are given by:
c                 vtx = vx + ax*t
c                 vty = vy + ay*t
c                 vtz = vz + az*t
c                 ptx = px + vx*t + 0.5*ax*t
c                 pty = py + vy*t + 0.5*ay*t
c                 ptz = pz + vz*t + 0.5*az*t
c
c               The times of any extrema of the particle path with the
c               quadric surface, tprox, are obtained by substituting ptx, pty
c               and ptz for x, y and z in the equation of the quadric surface,
c               differentiating with respect to time, and solving the resulting
c               cubic, quadratic or linear equation for the times.
c               There may be 0, 1, 2 or 3 such extrema.
c               The path distances dprox are found in aptparb.
c
c               Flag nerr indicates any input error.
c
c     History:
c
c     Input:    px, py, pz, vx, vy, vz, qc, qx, qy, qz, qxy, qyz, qzx,
c               qxx, qyy, qzz, tol. 
c
c     Output:   nprox, atype, ptx, pty, ptz, vtx, vty, vtz, vtlen,
c               dprox, tprox, nerr.
c
c     Calls: aptcubs, aptparb, aptqrts, aptquar 
c               
c
c     Glossary:
c
c     atype     Output   The extrema types, 'minimum' or 'maximum'.
c                          Size = 3.
c
c     ax,ay,az  Input    The x, y and z components of the constant acceleration
c                          vector "a".
c
c     dprox      Output   The distance along the particle path to the quadric
c                          surface, if nprox is positive.  Size = 3.
c
c     nerr      Output   Indicates an input error, if not zero.
c                          1 if equation for extrema times is null, indicating
c                          particle path is a point.  If so, nprox = 1, and
c                          extremum is the particle point.
c                          2 if equation for extrema times is impossible.
c
c     nprox     Output   The number of extrema of the particle path and
c                          the quadric surface.  Range 0 to 3.
c
c     ptx,y,z   Output   The x, y and z coordinates of the extrema on the
c                          particle path, if nprox is positive.  Size = 3.
c
c     pz,py,pz  Input    The initial x, y and z coordinates of the particle.
c
c     qc-qzz    Input    The coefficients of the implicit equation for the
c                        quartic surface.
c
c     tprox      Output  The times along the particle path to the extrema,
c                          if nprox is positive.  Size = 3.
c
c     tol       Input    Numerical tolerance limit.  With 64-bit floating point
c                          numbers, recommend 1.e-5 to 1.e-11.
c
c     vtx,y,z   Output   The x, y and z components of the particle velocity
c                          at points pt, if nprox is positive.  Size = 3.
c
c     vx,vy,vz  Input    The initial x, y and z components of the particle
c                          velocity.
ccend.

c.... Dimensioned arguments.

      dimension atype(3)              ! Extremum type 'minimum' or 'maximum'.
      character*8 atype               ! Extremum type 'minimum' or 'maximum'.
      dimension dprox(3)              ! Path distance to extremum.
      dimension dside(3)              ! Value of quadric equation.
      dimension ptx(3)                ! Coordinate x of int point.
      dimension pty(3)                ! Coordinate y of int point.
      dimension ptz(3)                ! Coordinate z of int point.
      dimension tprox(3)              ! Time to extremum.
      dimension timag(3)              ! Imaginary component of time.
      dimension vtlen(3)              ! Magnitude of particle velocity.
      dimension vtx(3)                ! Component x of particle velocity.
      dimension vty(3)                ! Component y of particle velocity.
      dimension vtz(3)                ! Component z of particle velocity.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptparx finding extrema of quadric surface' /
cbug     &  '  equation on particle path.',
cbug     &  '  tol=      ',1pe22.14)
cbug 9902 format ('  Particle at p with velocity v, acceleration a:' /
cbug     &  '  px,py,pz=   ',1p3e22.14 /
cbug     &  '  vx,vy,vz=   ',1p3e22.14 /
cbug     &  '  ax,ay,az=   ',1p3e22.14 )
cbug 9903 format ('  Quadric surface with implicit equation',
cbug     &  ' coefficients:' /
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) px, py, pz, vx, vy, vz, ax, ay, az
cbug      write ( 3, 9903) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.

c.... Initialize.

      nprox    = 0
      do 100 n = 1, 3
        atype(n) = 'unknown'
        dprox(n) = -1.e99
        ptx(n)   = -1.e99
        pty(n)   = -1.e99
        ptz(n)   = -1.e99
        tprox(n) = -1.e99
        timag(n) = -1.e99
        vtlen(n) = -1.e99
        vtx(n)   = -1.e99
        vty(n)   = -1.e99
        vtz(n)   = -1.e99
  100 continue

c.... Find the coefficients of the times in the implicit quadric equation.

      c0  = qc +
     &      qx *px    + qy *py    + qz *pz    +
     &      qxy*px*py + qyz*py*pz + qzx*pz*px +
     &      qxx*px*px + qyy*py*py + qzz*pz*pz

      e0  = abs (qc) +
     &      abs (qx *px)    + abs (qy *py)    + abs (qz *pz)    +
     &      abs (qxy*px*py) + abs (qyz*py*pz) + abs (qzx*pz*px) +
     &      abs (qxx*px*px) + abs (qyy*py*py) + abs (qzz*pz*pz)

      if (abs (c0) .le. tol * e0) c0 = 0.0

      c1  = qx*vx + qy*vy + qz*vz +
     &      qxy*(px*vy + py*vx) +
     &      qyz*(py*vz + pz*vy) +
     &      qzx*(pz*vx + px*vz) +
     &      2.0*(qxx*px*vx + qyy*py*vy + qzz*pz*vz)

      e1  = abs (qx*vx) + abs (qy*vy) + abs (qz*vz) +
     &      abs (qxy)*(abs (px*vy) + abs (py*vx)) +
     &      abs (qyz)*(abs (py*vz) + abs (pz*vy)) +
     &      abs (qzx)*(abs (pz*vx) + abs (px*vz)) +
     &      2.0*(abs (qxx*px*vx) + abs (qyy*py*vy) + abs (qzz*pz*vz))

      if (abs (c1) .le. tol * e1) c1 = 0.0

      c2  = 0.5*    (qx*ax +     qy*ay + qz*az) +
     &      0.5*qxy*(px*ay + 2.0*vx*vy + py*ax) +
     &      0.5*qyz*(py*az + 2.0*vy*vz + pz*ay) +
     &      0.5*qzx*(pz*ax + 2.0*vz*vx + px*az) +
     &          qxx*(px*ax +     vx*vx)         +
     &          qyy*(py*ay +     vy*vy)         +
     &          qzz*(pz*az +     vz*vz)

      e2  = 0.5*     ((abs (qx*ax) +     abs (qy*ay) + abs (qz*az))  +
     &      abs (qxy)*(abs (px*ay) + 2.0*abs (vx*vy) + abs (py*ax))  +
     &      abs (qyz)*(abs (py*az) + 2.0*abs (vy*vz) + abs (pz*ay))  +
     &      abs (qzx)*(abs (pz*ax) + 2.0*abs (vz*vx) + abs (px*az))) +
     &      abs (qxx)*(abs (px*ax) +     abs (vx*vx))         +
     &      abs (qyy)*(abs (py*ay) +     abs (vy*vy))         +
     &      abs (qzz)*(abs (pz*az) +     abs (vz*vz))

      if (abs (c2) .le. tol * e2) c2 = 0.0

      c3  = 0.5*qxy*(vx*ay + vy*ax) +
     &      0.5*qyz*(vy*az + vz*ay) +
     &      0.5*qzx*(vz*ax + vx*az) +
     &          qxx*vx*ax + qyy*vy*ay + qzz*vz*az

      e3  = 0.5*abs (qxy)*(abs (vx*ay) + abs (vy*ax)) +
     &      0.5*abs (qyz)*(abs (vy*az) + abs (vz*ay)) +
     &      0.5*abs (qzx)*(abs (vz*ax) + abs (vx*az)) +
     &          abs (qxx*vx*ax) + abs (qyy*vy*ay) + abs (qzz*vz*az)

      if (abs (c3) .le. tol * e3) c3 = 0.0

      c4  = 0.25*(qxy*ax*ay + qyz*ay*az + qzx*az*ax +
     &            qxx*ax*ax + qyy*ay*ay + qzz*az*az)

      e4  = 0.25*(abs (qxy*ax*ay) + abs (qyz*ay*az) + abs (qzx*az*ax) +
     &            abs (qxx*ax*ax) + abs (qyy*ay*ay) + abs (qzz*az*az))

      if (abs (c4) .le. tol * e4) c4 = 0.0

c.... Find any extrema of the particle path and the quadric.

      d0 = c1
      d1 = 2.0 * c2
      d2 = 3.0 * c3
      d3 = 4.0 * c4

cbugcbugc***DEBUG begins.
cbugcbug 9904   format ('Coefficients of the equation for the',
cbugcbug     &    ' extremum times:' /
cbugcbug     &          '  set d0         =',1pe20.12 /
cbugcbug     &          '  set d1         =',1pe20.12 /
cbugcbug     &          '  set d2         =',1pe20.12 /
cbugcbug     &          '  set d3         =',1pe20.12 )
cbugcbug
cbugcbug      write ( 3, 9904) d0, d1, d2, d3
cbugcbugc***DEBUG ends.

      if (d3 .ne. 0.0) then       ! Solve cubic equation.
cbugcbugc***DEBUG begins.
cbugcbug 9906 format ('aptparx DEBUG.  Solving cubic.')
cbugcbug      write ( 3, 9906)
cbugcbugc***DEBUG ends
        d0 = d0 / d3
        d1 = d1 / d3
        d2 = d2 / d3

        call aptcubs (d0, d1, d2, tol,
     &                nprox, tprox, timag, atype, itrun)

      elseif (d2 .ne. 0.0) then       ! Solve quadratic equation.
cbugcbugc***DEBUG begins.
cbugcbug 9907 format ('aptparx DEBUG.  Solving quadratic.')
cbugcbug      write ( 3, 9907)
cbugcbugc***DEBUG ends

        call aptqrts (0, d2, d1, d0, qq, tol,
     &                nprox, tprox(1), tprox(2), itrun)

      elseif (d1 .ne. 0.0) then       ! Solve linear equation.
cbugcbugc***DEBUG begins.
cbugcbug 9908 format ('aptparx DEBUG.  Solving linear.')
cbugcbug      write ( 3, 9908)
cbugcbugc***DEBUG ends
        nprox   = 1
        tprox(1) = -d0 / d1
      elseif (d0 .ne. 0.0) then       ! Impossible equation.
        nerr = 2
cbugcbugc***DEBUG begins.
cbugcbug 9909 format ('aptparx DEBUG.  Impossible equation.')
cbugcbug      write ( 3, 9909)
cbugcbugc***DEBUG ends
      else                            ! Null equation.  Path is a point.
        nerr     = 1
        nprox    = 1
        tprox(1) = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 9910 format ('aptparx DEBUG.  Null equation.')
cbugcbug      write ( 3, 9910)
cbugcbugc***DEBUG ends
      endif                           ! Tested order of equation for time.

c.... Find path, speed, velocity and position at extrema.

      if (nprox .gt. 0) then

        do 850 n = 1, nprox           ! Loop over times of extrema.

          curve   = d1 + 2.0 * d2 * tprox(n) + 3.0 * d3 * tprox(n)**2
          curverr = tol * (abs (d1) + 2.0 * abs (d2 * tprox(n)) +
     &                     3.0 * abs (d3) * tprox(n)**2)
          if (abs (curve) .le. curverr) curve = 0.0

          if (curve .lt. 0.0) then
            atype(n) ='maximum'
          elseif (curve .eq. 0.0) then
            atype(n) = 'min/max'
          else
            atype(n) ='minimum'
          endif

          call aptparb (px, py, pz, vx, vy, vz, ax, ay, az,
     &                  tprox(n), tol,
     &                  ptx(n), pty(n), ptz(n),
     &                  vtx(n), vty(n), vtz(n),
     &                  vtlen(n), dprox(n), ntype)

          ta       = tprox(n)
          dside(n) = c0 + ta * (c1 + ta * (c2 + ta * (c3 + ta * c4)))
          ta       = abs (tprox(n))
          eside    = tol * (abs (c0) + ta * (abs (c1) + ta *
     &               (abs (c2) + ta * (abs (c3) + ta * abs (c4)))))
          if (abs (dside(n)) .lt. eside) dside(n) = 0.0

  850   continue                      ! End of loop over times of extrema.

      endif

cbugc***DEBUG begins.
cbug 9912 format (/ 'aptparx results:  nerr=',i2,' nprox=',i2 )
cbug      write ( 3, 9912) nerr, nprox
cbug      if (nprox .gt. 0) then          ! One or more extrema.
cbug 9913   format (/ 'Extremum type   ',a8 /
cbug     &            'Time, path, speed ',1p3e20.12 /
cbug     &            'Velocity (xyz)    ',1p3e20.12 /
cbug     &            'Position (xyz)    ',1p3e20.12 /
cbug     &            'dside             ',1pe20.12 )
cbug          write ( 3, 9913) (atype(n), tprox(n), dprox(n), vtlen(n),
cbug     &                        vtx(n), vty(n), vtz(n),
cbug     &                        ptx(n), pty(n), ptz(n), dside(n),
cbug     &                        n = 1, nprox)
cbug      endif                           ! Tested number of extrema.
cbugc***DEBUG ends.

      return

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

UCRL-WEB-209832