subroutine aptrkcy (pz, pr, uz, ur, ut, az, ar, bz, br,
     &                    dintmn, dintmx, np, noptd, tolf, tols,
     &                    nint, pinz, pinr, dint, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTRKCY
c
c     call aptrkcy (pz, pr, uz, ur, ut, az, ar, bz, br,
c    &              dintmn, dintmx, np, noptd, tolf, tols,
c    &              nint, pinz, pinr, dint, nerr)
c
c     Version:  aptrkcy  Updated    1993 December 6 15:20.
c               aptrkcy  Originated 1989 December 7 16:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find, for each of the np sets of input data, any acceptable
c               exit intersection point pin = (pinz, pinr) of the linear track
c               through point p = (pz, pr) with the unit direction vector
c               u = (uz, ur, ut), with the conical surface symmetric around the
c               z axis, and through the points a = (az, ar) and b = (bz, br),
c               for which:  (1) the distance from point "p" to point "pin" is
c               between the limits dintmn and dintmx, (2) the intersection is
c               between points "a" and "b", and (3) the crossing at the
c               intersection is from left to right in the z-r plane (right to
c               left in the r-z plane).  This is a zone exit if the points "a"
c               and "b" are the vertices of a zone edge, counterclockwise around
c               the zone in the the z-r plane (clockwise in the r-z plane).
c               Flag nerr indicates any input error, if not zero.
c
c               The coordinate system is 3-D cylindrical, with the axial
c               coordinate z, the radial coordinate r = sqrt (x**2 + y**2),
c               and the angular coordinate theta (counterclockwise in the x-y
c               plane from the x axis.
c
c     History:  1990 January 23 10:10.  Implemented use of aptfdav to adjust
c               values of fint1 and fint2 near 0.0 or 1.0.  Vectorized the
c               intersection test loop.
c
c     Input:    pz, pr, uz, ur, ut, az, ar, bz, br, dintmn, dintmx, np,
c               noptd, tolf, tols.
c
c     Output:   nint, pinz, pinr, dint, nerr.
c
c     Calls: aptqrtv, aptvlim, aptfdav 
c               
c
c     Glossary:
c
c     az, ar    Input    Coordinates z and r of the beginning of the surface.
c                          Size np.
c
c     bz, br    Input    Coordinates z and r of the end of the surface.
c                          Size np.
c
c     dint      Output   Distance from point "p" to the intersection at point
c                          "pin" along the track (not distance in z-r plane).
c                          Size np.
c
c     dintmn    Input    Minimum allowed value of the distance to intersection
c                          dint.  Size np.
c
c     dintmx    Input    Maximum allowed value of the distance to intersection
c                          dint.  Size np.
c
c     nerr      Output   If not 0, indicates an input error.
c                          1 if np is not positive.
c
c     nint      Output   O if no acceptable exit intersection was found.
c                          1 if an acceptable exit intersection was found on the
c                          surface segment between points "a" and "b", with a
c                          distance to the intersection dint between dintmn and
c                          dintmx.  Size np.
c
c     noptd     Input    Option to limit the minimum magnitude of the components
c                          of the direction vector u = (uz, ur, ut).  The method
c                          used here will fail if uz or ut is zero.
c                          0 for no magnitude test.  Must be done elsewhere.
c                          1 to limit the minimum magnitudes of ut and uz to no
c                            less than tols.
c
c     np        Input    Size of arrays pz, pr, uz, ur, ut, az, ar, bz, br,
c                          dintmn, dintmx, nint, pinz, pinr.
c
c     pinr      Output   Coordinate r of the intersection point "pin".  Size np.
c
c     pinz      Output   Coordinate z of the intersection point "pin".  Size np.
c
c     pz, pr    Input    Coordinates z and r of the initial point on the linear
c                          track.  Size np.
c
c     uz,ur,ut  Input    Initial z, r and theta components of the unit direction
c                          vector along the linear track.  Size np.
c
c     tolf      Input    Truncation error limit to be imposed on the fractional
c                          distance of the intersection point "pin" along
c                          the line segment from point "a" to point "b".  Values
c                          less than -tolf or greater than 1.0 + tolf will not
c                          be accepted.  Values from -tolf to tolf will be
c                          changed to tolf.  Values from 1.0 - tolf to
c                          1.0 + tolf will be changed to 1.0 - tolf.
c                          Also used to test for the intersection being nearly
c                          tangent, and for the accuracy of the intersection.
c                          Must not be zero.
c                          On Cray computers, recommend 1.e-11.
c
c     tols      Input    Truncation error limit to be imposed on uz, ut and
c                          pinr.  Magnitudes of ut and uz less than tols
c                          will be increased to tols.  Values of pinr less
c                          than tols * pr will be increased to tols * pr.
c                          On Cray computers, recommend 1.e-5.
c                          A value of zero may produce unpredictable results.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Coordinate r at start of edge "ab".
      dimension ar      (1)
c---- Coordinate z at start of edge "ab".
      dimension az      (1)
c---- Coordinate r at end of edge "ab".
      dimension br      (1)
c---- Coordinate z at end of edge "ab".
      dimension bz      (1)
c---- Distance to intersection "pin".
      dimension dint    (1)
c---- Minimum distance to intersection.
      dimension dintmn  (1)
c---- Maximum distance to intersection.
      dimension dintmx  (1)
c---- Number of intersections (0 or 1).
      dimension nint    (1)
c---- Coordinate r of intersection "pin".
      dimension pinr    (1)
c---- Coordinate z of intersection "pin".
      dimension pinz    (1)
c---- Coordinate r of beginning of track.
      dimension pr      (1)
c---- Coordinate z of beginning of track.
      dimension pz      (1)
c---- Component r of direction vector "u".
      dimension ur      (1)
c---- Component theta of "u".
      dimension ut      (1)
c---- Component z of direction vector "u".
      dimension uz      (1)

c.... Local variables.

c---- Coefficient of fint**2 in quadratic.
      common /laptrkcy/ aa      (64)
c---- Coefficient of fint in quadratic.
      common /laptrkcy/ bb      (64)

c---- A very big number.
      common /laptrkcy/ big

c---- Factor 2.0 * pr * ur / uz.
      common /laptrkcy/ c1      (64)
c---- Factor (ur**2 + ut**2) / uz**2.
      common /laptrkcy/ c2      (64)
c---- Constant term in quadratic.
      common /laptrkcy/ cc      (64)
c---- Distance to intersection 1.
      common /laptrkcy/ dint1
c---- Distance to intersection 2.
      common /laptrkcy/ dint2
c---- Difference br - ar.
      common /laptrkcy/ dr      (64)
c---- Difference bz - az.
      common /laptrkcy/ dz      (64)
c---- First root of quadratic.
      common /laptrkcy/ fint1   (64)
c---- Second root of quadratic.
      common /laptrkcy/ fint2   (64)

c---- A very small number.
      common /laptrkcy/ fuz

c---- 1 = tangent, 2 = bad qq.
      common /laptrkcy/ itrun   (64)
c---- Index in local array.
      common /laptrkcy/ n
c---- First index of subset of data.
      common /laptrkcy/ n1
c---- Last index of subset of data.
      common /laptrkcy/ n2
c---- 1 if intersection 1 good exit.
      common /laptrkcy/ nexit1
c---- 1 if intersection 2 good exit.
      common /laptrkcy/ nexit2
c---- 1 if fint1 or fint 2 adjusted.
      common /laptrkcy/ nlim    (64)
c---- Index in external array.
      common /laptrkcy/ nn
c---- Number of real roots of quadratic.
      common /laptrkcy/ nroots  (64)
c---- Size of current subset of data.
      common /laptrkcy/ ns
c---- Coordinate r at intersection 1.
      common /laptrkcy/ pinr1
c---- Coordinate r at intersection 2.
      common /laptrkcy/ pinr2
c---- Coordinate z at intersection 1.
      common /laptrkcy/ pinz1
c---- Coordinate z at intersection 2.
      common /laptrkcy/ pinz2
c---- Factor bb**2 - 4.0 * aa * cc.
      common /laptrkcy/ qq      (64)
c---- New ur at intersection 1.
      common /laptrkcy/ ur1
c---- New ur at intersection 2.
      common /laptrkcy/ ur2
c---- Magnitude of direction vector.
      common /laptrkcy/ vlen    (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrkcy finding intersections.  np=',i3,
cbug     &  ' noptd=',i2 / '  tolf=',1pe22.14,' tols=',1pe22.14)
cbug 9902 format (i3,' pz,pr=    ',1p2e22.14 /
cbug     &  '    uz,ur,ut= ',1p3e22.14 /
cbug     &  '    az,ar=    ',1p2e22.14 /
cbug     &  '    bz,br=    ',1p2e22.14 /
cbug     &  '    dintmn=   ',1pe22.14,' dintmx=',1pe22.14)
cbug      write ( 3, 9901) np, noptd, tolf, tols
cbug      write ( 3, 9902) (nn, pz(nn), pr(nn), uz(nn), ur(nn), ut(nn),
cbug     &  az(nn), ar(nn), bz(nn), br(nn), dintmn(nn), dintmx(nn),
cbug     &  nn = 1, np)
cbugc***DEBUG ends.

c.... Initialize.

c---- A very big number.
      big = 1.e+99

c---- A very small number.
      fuz = 1.e-99

c.... Test for input errors.

      nerr = 0
      if (np .le. 0) then
        nerr = 1
        go to 210
      endif

c.... Do data sets in groups of 64.

      n1 = 1
      n2 = min (np, 64)
  110 ns = n2 - n1 + 1

c.... Do not allow any tracks that are linear in the z-r plane.

c---- Limit minimum size of uz, ut.
      if (noptd .ne. 0) then

        call aptvlim (uz(n1), ur(n1), ut(n1), ns, tols, 0., tols,
     &                vlen, nerr)

c---- Tested noptd.
      endif

c.... Find the coefficients of aa * fint**2 + bb * fint + cc = 0, where
c....   fint = ((pinz - az) * dz + (pinr - ar) * dr) / (dr**2 + dz**2),
c....   by substituting pinz = az + fint * dz, pinr = ar + fint * dr, in
c....   pinr = sqrt (pr**2 + c1 * (pinz - pz) + c2 * (pinz - pz)**2).

c---- Loop over local arrays.
      do 120 n = 1, ns
        nn    = n + n1 - 1
        c1(n) = 2.0 * pr(nn) * ur(nn) / uz(nn)
        c2(n) = (ur(nn)**2 + ut(nn)**2) / uz(nn)**2
        dr(n) = br(nn) - ar(nn)
        dz(n) = bz(nn) - az(nn)
        aa(n) = c2(n) * dz(n)**2 - dr(n)**2
        bb(n) = (c1(n) + 2.0 * c2(n) * (az(nn) - pz(nn))) * dz(n) -
     &          2.0 * ar(nn) * dr(n)
        cc(n) = pr(nn)**2 - ar(nn)**2 + (az(nn) - pz(nn)) *
     &          (c1(n) + c2(n) * (az(nn) - pz(nn)))

c.... The following calculation of qq is more accurate than
c....   qq(n) = bb(n)**2 - 4.0 * aa(n) * cc(n), due to terms cancelling.

        qq(n) = 4.0 * dr(n) * (dr(n) *
     &          (pr(nn)**2 + (az(nn) - pz(nn)) *
     &          (c1(n) + c2(n) * (az(nn) - pz(nn)))) -
     &          dz(n) * ar(nn) * (c1(n) +
     &          2.0 * c2(n) * (az(nn) - pz(nn)))) + dz(n)**2 *
     &          (c1(n)**2 - 4.0 * c2(n) * (pr(nn)**2 - ar(nn)**2))
c---- End of loop over local arrays.
  120 continue

c.... Find any intersections.

      call aptqrtv (1, aa, bb, cc, qq, ns, tolf,
     &              nroots, fint1, fint2, itrun, nerr)

c.... Adjust any values of fint1 and fint2 near 0.0 or 1.0.

      call aptfdav (fint1, ns, 1, tolf, nlim, nerr)
      call aptfdav (fint2, ns, 1, tolf, nlim, nerr)

c.... See if either intersection is on zone edge, on track,
c....   and an exit from zone.

c---- Loop over local arrays.
      do 130 n = 1, ns

        nn    = n + n1 - 1
        pinz1 = az(nn) + fint1(n) * dz(n)
        pinz2 = az(nn) + fint2(n) * dz(n)
        pinr1 = amax1 (pr(nn) * tols, ar(nn) + fint1(n) * dr(n))
        pinr2 = amax1 (pr(nn) * tols, ar(nn) + fint2(n) * dr(n))
        dint1 = (pinz1 - pz(nn)) / uz(nn)
        dint2 = (pinz2 - pz(nn)) / uz(nn)

        ur1   = (pr(nn) * ur(nn) + (ur(nn)**2 + ut(nn)**2) * dint1) /
     &          (pinr1 + fuz)
        ur2   = (pr(nn) * ur(nn) + (ur(nn)**2 + ut(nn)**2) * dint2) /
     &          (pinr2 + fuz)

        if ((nroots(n) .ge. 1)      .and.
     &      (fint1(n) .ge. 0.0)     .and.
     &      (fint1(n) .le. 1.0)     .and.
     &      (dint1 .ge. dintmn(nn)) .and.
     &      (dint1 .le. dintmx(nn)) .and.
     &      ((uz(nn) * dr(n)) .gt. (ur1 * dz(n)))) then
          nexit1 = 1
        else
          nexit1 = 0
        endif

        if (nexit1 .eq. 1) then
          nroots(n) = 1
        endif

        if (nexit1 .eq. 1) then
          nint(nn) = 1
          pinr(nn) = pinr1
          pinz(nn) = pinz1
          dint(nn) = dint1
        else
          nint(nn) = 0
          pinr(nn) = -big
          pinz(nn) = -big
          dint(nn) = -big
        endif

        if ((nroots(n) .eq. 2)      .and.
     &      (fint2(n) .ge. 0.0)     .and.
     &      (fint2(n) .le. 1.0)     .and.
     &      (dint2 .ge. dintmn(nn)) .and.
     &      (dint2 .le. dintmx(nn)) .and.
     &      ((uz(nn) * dr(n)) .gt. (ur2 * dz(n)))) then
          nexit2 = 1
        else
          nexit2 = 0
        endif

        if (nexit2 .eq. 1) then
          nint(nn) = 1
          pinr(nn) = pinr2
          pinz(nn) = pinz2
          dint(nn) = dint2
        endif

c---- End of loop over local arrays.
  130 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptrkcy results:')
cbug 9904 format (i3,' nint=',i2 /
cbug     &  '    z,r,dint= ',1p3e22.14 /
cbug     &  '    nn=',i3 /
cbug     &  '    fint(1,2)=',1p2e22.14 /
cbug     &  '    c1,c2,qq= ',1p3e22.14 /
cbug     &  '    aa,bb,cc= ',1p3e22.14)
cbug      write ( 3, 9903)
cbug      do 140 n = 1, ns
cbug        nn = n + n1 - 1
cbug        write ( 3, 9904) nn, nint(nn), pinz(nn), pinr(nn), dint(nn),
cbug     &    n, fint1(n), fint2(n),
cbug     &    c1(n), c2(n), qq(n), aa(n), bb(n), cc(n)
cbug  140 continue
cbugc***DEBUG ends.

c.... See if all sets are done.

c---- Do another data set.
      if (n2 .lt. np) then
        n1 = n2 + 1
        n2 = min (np, n1 + 63)
        go to 110
      endif

  210 return

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

UCRL-WEB-209832