subroutine aptparp (px, py, pz, vx, vy, vz, ax, ay, az,
     &                    bx, by, bz, cx, cy, cz, tol,
     &                    nints, ptx, pty, ptz, vtx, vty, vtz, vtlen,
     &                    dint, tint, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTPARP
c
c     call aptparp (px, py, pz, vx, vy, vz, ax, ay, az,
c    &              bx, by, bz, cx, cy, cz, tol,
c    &              nints, ptx, pty, ptz, vtx, vty, vtz, vtlen,
c    &              dint, tint, nerr)
c
c     Version:  aptparp  Updated    1997 April 16 16:20.
c               aptparp  Originated 1996 December 19 18:00.
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 plane through point
c               b = (bx, by, bz) with normal vector c = (cx, cy, cz), to find
c               the number nints (0, 1 or 2) of intersections of the particle
c               path with the plane, and for each such intersection point,
c               the coordinates pt(n) = (ptx(n), pty(n), ptz(n)), n = 1, nints,
c               the velocities vt(n) = (vtx(n), vty(n), vtz(n)), n = 1, nints,
c               the velocity magnitudes vtlen(n), n = 1, nints,
c               the path lengths to intersection dint(n), n = 1, nints, and
c               the times to intersection tint(n), n = 1, nints.
c               If no such intersections occur (nints = 0), to find the proximal
c               points pt(1) = (ptx(1), pty(1), ptz(1)) on the particle path and
c               pt(2) = (ptx(2), pty(2), ptz(2)) on the plane, the velocity
c               vt(1) = (vtx(1), vty(1), vtz(1)) and its magnitude vtlen(1) at
c               the proximal point, the path length to the proximal point
c               dint(1), the distance between the proximal points dint(2),
c               and the time tint(1) to reach the proximal point.
c               Flag nerr indicates any input error.
c
c               The particle velocity and coordinates are given by:
c                 vt = v + a*t,
c                 pt = p + v*t + 0.5*a*t**2.
c
c               The distance from the particle to the plane is given by:
c                 dp = (pt - b) dot c / |c|, or
c                 dp = [(p - b) dot c + (v dot c)*t + 0.5*(a dot c)*t**2] / |c|,
c               where "dot" indicates the vector dot product, and |c| is the
c               absolute magnitude of vector c.
c
c               If any intersections between the particle path and the plane
c               occur, they may be found by setting dp = 0 and solving for
c               t = tint in the equation above, then substituting t = tint
c               in the equation for pt to find d.
c               There may be zero, one or two such solutions.
c               If no solutions for t exist, then either the path is parallel
c               to the plane, with (v dot c) = 0 and (a dot c) = 0, or the
c               path is at its nearest point to the plane when is it parallel
c               to the plane:
c                 vt dot c = 0, or
c                 v dot c + (a dot c) tint = 0, or
c                 tint = -(v dot c) / (a dot c),
c               which may be used to find d = pt at the proximal point.
c
c               See subroutine aptparb for finding dint.
c
c     History:
c
c     Input:    px, py, pz, vx, vy, vz, ax, ay, az, bx, by, bz, cx, cy, cz, tol.
c
c     Output:   nints, ptx, pty, ptz, vtx, vty, vtz, vtlen, dint, tint, nerr.
c
c     Calls: aptparb, aptptpl, aptqrts, aptvdis, aptvdot, aptvsum, aptvunb 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y and z components of the constant acceleration
c                          vector "a".
c
c     bx,by,bz  Input    The x, y and z coordinates of point "b" in the plane.
c
c     cx,cy,cz  Input    The x, y and z components of the vector "c" normal to
c                          the plane.
c
c     dint      Output   The distance along the particle path to the proximal
c                          point to the plane, and from the plane to the
c                          particle path (nints = 0), or to the intersection(s)
c                          with the plane (nints = 1 or 2).  Size = 2.
c
c     nerr      Output   Indicates an input error, if not zero.
c                          1 if normal vector c has zero magnitude.
c
c     nints     Output   The number of intersections of the particle path and
c                          the plane.  0 if no intersections occur.
c                          -1 if the intersection is a tangent point.
c
c     ptx,y,z   Output   The x, y and z coordinates of the proximal points on
c                          the particle path and in the plane (nints = 0),
c                          or the point(s) of intersection of the particle path
c                          and the plane (nints = 1 or 2).  Size 2.
c
c     pz,py,pz  Input    The initial x, y and z coordinates of the particle.
c
c     tint      Output   The time along the particle path to the proximal
c                          point to the plane (nints = 0), or to the
c                          intersection(s) with the plane (nints = 1 or 2).
c                          Size = 2.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
c     vtlen     Output   The magnitudes of velocities vt.  Size 2.
c
c     vtx,y,z   Output   The x, y and z components of the particle velocity
c                          at points pt.  Size = 2.
c
c     vx,vy,vz  Input    The initial x, y and z components of the particle
c                          velocity.
ccend.

c.... Dimensioned arguments.

      dimension ptx(2)                ! Coordinate x of prox or int point.
      dimension pty(2)                ! Coordinate y of prox or int point.
      dimension ptz(2)                ! Coordinate z of prox or int point.
      dimension vtlen(2)              ! Magnitude of particle velocity.
      dimension vtx(2)                ! Component x of particle velocity.
      dimension vty(2)                ! Component y of particle velocity.
      dimension vtz(2)                ! Component z of particle velocity.
      dimension dint(2)               ! Path distance to prox or intersection.
      dimension tint(2)               ! Time to proximity or intersection.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptparp finding intersection or proximal point' /
cbug     &  '  of particle path and plane.  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 ('  Plane through point b with normal vector c:' /
cbug     &  '  bx,by,bz=   ',1p3e22.14 /
cbug     &  '  cx,cy,cz=   ',1p3e22.14 )
cbug      write ( 3, 9901) tol
cbug      write ( 3, 9902) px, py, pz, vx, vy, vz, ax, ay, az
cbug      write ( 3, 9903) bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.

c.... Initialize.

      nints    = 0
      do 100 n = 1, 2
        dint(n)  = -1.e99
        ptx(n)   = -1.e99
        pty(n)   = -1.e99
        ptz(n)   = -1.e99
        tint(n)  = -1.e99
        vtlen(n) = -1.e99
        vtx(n)   = -1.e99
        vty(n)   = -1.e99
        vtz(n)   = -1.e99
  100 continue

c.... Find the magnitude of the normal vector "c".  Test for error.

      call aptvunb (cx, cy, cz, 1, tol, unx, uny, unz, clen, nerr)

      if (clen .eq. 0.0) then
        nerr = 1
        go to 210
      endif

c.... Find the magnitude of the initial velocity vector "v".

      call aptvunb (vx, vy, vz, 1, tol, uvx, uvy, uvz, vlen, nerr)

c.... Find the magnitude of the acceleration vector "a".

      call aptvunb (ax, ay, az, 1, tol, uax, uay, uaz, alen, nerr)

c.... Find the vector pb from point "p" to point "b".

      call aptvdis (px, py, pz, bx, by, bz, 1, tol,
     &              pbx, pby, pbz, pblen, nerr)

c.... Find the dot product of vector "pb" and the unit normal vector "un".

      call aptvdot (pbx, pby, pbz, unx, uny, unz, 1, tol, pbdotu, nerr)

c.... Find the dot product of vector "v" and the unit normal vector "un".

      call aptvdot (vx, vy, vz, unx, uny, unz, 1, tol, vdotu, nerr)

c.... Find the dot product of vector "a" and the unit normal vector "un".

      call aptvdot (ax, ay, az, unx, uny, unz, 1, tol, adotu, nerr)

c.... Find any intersections, or if none, find the proximal points.

      if (adotu .eq. 0.0) then        ! Acceleration is parallel to plane.
        if (vdotu .eq. 0.0) then      ! Path always parallel to plane.
          nints    = 0
          tint(1)  = 0.0
          vtx(1)   = vx
          vty(1)   = vy
          vtz(1)   = vz
          vtlen(1) = vlen
          dint(1)  = 0.0
          ptx(1)   = px
          pty(1)   = py
          ptz(1)   = pz

          call aptptpl (px, py, pz, bx, by, bz, cx, cy, cz, 1, tol,
     &                  dint(2), ptx(2), pty(2), ptz(2), itrun, nerr)

        else                          ! Path intersects plane once.
          nints    = 1
          tint(1)  = pbdotu / vdotu

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

        endif                         ! Tested direction of initial velocity.
      else                            ! Acceleration is not parallel to plane.
        if (vdotu .eq. 0.0) then      ! Path initially parallel to plane.
          arg = 2.0 * pbdotu / adotu
          if (arg .lt. 0.0) then      ! No intersection.
            nints    = 0
            tint(1)  = 0.0
            vtx(1)   = vx
            vty(1)   = vy
            vtz(1)   = vz
            vtlen(1) = vlen
            dint(1)  = 0.0
            dint(2)  = -pbdotu
            ptx(1)   = px
            pty(1)   = py
            ptz(1)   = pz

            call aptptpl (px, py, pz, bx, by, bz, cx, cy, cz, 1, tol,
     &                    dint(2), ptx(2), pty(2), ptz(2), itrun, nerr)

          else                        ! Two intersections.
            nints    = 2
            tint(1)  = -sqrt (arg)
            tint(2)  =  sqrt (arg)

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

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

          endif                       ! Tested number of intersections.
        else                          ! Path is parabolic, not at vertex.
          aa = 0.5 * adotu
          bb = vdotu
          cc = -pbdotu

          call aptqrts (0, aa, bb, cc, qq, tol,
     &                  nroots, troot1, troot2, itrun)

          if (nroots .eq. 0) then     ! No intersections.
            nints   = 0
            tint(1) = -vdotu / adotu

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

            call aptptpl (ptx(1), pty(1), ptz(1),
     &                    bx, by, bz, cx, cy, cz, 1, tol,
     &                    dint(2), ptx(2), pty(2), ptz(2), itrun, nerr)


          elseif (nroots .eq. 1) then
            nints   = 1
            tint(1) = troot1

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

c....       See if the intersection is a tangent point.

            call aptvdot (vtx(1), vty(1), vtz(1), unx, uny, unz,
     &                    1, tol, vtdotu, nerr)

            if (vtdotu .eq. 0.0) nints = -1

          else
            nints   = 2
            tint(1) = troot1
            tint(2) = troot2

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

            call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(2),
     &                    tol, ptx(2), pty(2), ptz(2),
     &                    vtx(2), vty(2), vtz(2), vtlen(2),
     &                    dint(2), ntype)
          endif                       ! Tested number of roots.
        endif                         ! Tested direction of initial velocity.
      endif                           ! Tested direction of acceleration.

  210 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptparp results:  nerr=',i2,' nints=',i2 /
cbug     &  '  path1,tint1=',1p2e22.14 /
cbug     &  '  d1(x,y,z)=  ',1p3e22.14 /
cbug     &  '  path2,tint2=',1p2e22.14 /
cbug     &  '  d2(x,y,z)=  ',1p3e22.14 /
cbug     &  '  v1(x,y,z)=  ',1p3e22.14 /
cbug     &  '  vtlen1=     ',1p1e22.14 /
cbug     &  '  v2(x,y,z)=  ',1p3e22.14 /
cbug     &  '  vtlen2=     ',1p1e22.14 )
cbug      write ( 3, 9904) nerr, nints,
cbug     &                 dint(1), tint(1), ptx(1), pty(1), ptz(1),
cbug     &                 dint(2), tint(2), ptx(2), pty(2), ptz(2),
cbug     &                 vtx(1), vty(1), vtz(1), vtlen(1),
cbug     &                 vtx(2), vty(2), vtz(2), vtlen(2)
cbugc***DEBUG ends.

      return

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

UCRL-WEB-209832