subroutine aptmvcy (kth, psz, psr, cths, sths, usz, usr, ust,
     &                    dpmove, np, tol,
     &                    pz, pr, cth, sth, uz, ur, ut, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTMVCY
c
c     call aptmvcy (kth, psz, psr, cths, sths, usz, usr, ust,
c    &              dpmove, np, tol,
c    &              pz, pr, cth, sth, uz, ur, ut, nerr)
c
c     Version:  aptmvcy  Updated    1990 November 28 14:50.
c               aptmvcy  Originated 1989 December 4 11:00.
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, the new
c               point p = (pz, pr, cth, sth) and unit direction vector
c               u = (uz, ur, ut), resulting from moving from the point
c               ps = (psz, psr, cths, sths) in the direction of the unit
c               vector us = (usz, usr, ust) for a distance dpmove, in
c               cylindrical coordinates.
c               If kth = 0, all cths = 1.0, all sths = 0.0, and cth and sth
c               will not be calculated, and none need be dimensioned.
c               Any component of point "p" or vector "u" less than the
c               estimated error in its calculation, based on tol, will be
c               truncated to zero.  If tol = 0, no truncation tests are done.
c               Flag nerr indicates any input error.
c
c     History:  1990 March 14.  Changed tol to 0.0 in call to unit vector
c               subroutine.  Allows small magnitudes.
c               1990 March 30.  Fixed array index error which affected problems
c               with np .gt. 64.
c
c     Input:    kth, psz, psr, cths, sths, usz, usr, ust, dpmove, np, tol.
c
c     Output:   pz, pr, cth, sth, uz, ur, ut, nerr.
c
c     Calls: aptvuna 
c
c     Glossary:
c                          Size np.
c
c     cth       Output   The cosine of the final value of theta (angle in the
c                          xy plane counterclockwise from x axis).
c                          May be truncated to zero, if less than the estimated
c                          numerical error in their calculation based on tol.
c                          Size np, if kth = 1.  Otherwise, not calculated.
c
c     cths      Input    The cosine of the initial value of theta.  See cth.
c                          Size np, if kth = 1.  Otherwise, scalar 1.0.
c
c     dpmove    Input    The distance from point "ps" to point "p".  Size np.
c                          (Assuming vector "us" is a unit vector.)
c
c     kth       Input    Indicates size of arrays cth, cths, sth, sths:
c                          0 if array size is 1, with cths = 1.0, sths = 0.0,
c                            and cth and sth are not to be calculated.
c                          1 if array size is np, input values of cths and sths
c                            will be used, and cth and sth will be calculated.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if kth is not 0 or 1.
c
c     np        Input    Size of arrays.
c
c     pr, pz    Output   The r and z components of final point "p".  Size np.
c                          May be truncated to zero, if less than the estimated
c                          numerical error in their calculation based on tol.
c
c     psr, psz  Input    The r and z coordinates of initial point "ps".
c                          Size np.
c
c     sth       Output   The sine of the final value of theta.  See cth.
c                          May be truncated to zero, if less than the estimated
c                          numerical error in their calculation based on tol.
c                          Size np, if kth = 1.  Otherwise, not calculated.
c
c     sths      Input    The sine of the initial value of theta.  See cth.
c                          Size np, if kth = 1.  Otherwise, scalar 0.0.
c
c     tol       Input    Numerical tolerance limit.  If zero, no tests done.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
c     ur        Output   The r component of final unit direction vector "u".
c                          May be truncated to zero, if less than the estimated
c                          numerical error in their calculation based on tol.
c                          Size np.
c
c     usr       Input    The r component of initial unit direction vector "us".
c                          Size np.
c
c     ust       Input    The theta component of initial unit direction vector
c                          "us".  See cth.  Size np.
c
c     usz       Input    The z component of initial unit direction vector "us".
c                          Size np.
c
c     ut        Output   The theta component of final unit direction vector "u".
c                          May be truncated to zero, if less than the estimated
c                          numerical error in their calculation based on tol.
c                          See cth.  Size np.
c
c     uz        Output   The z component of final unit direction vector "u".
c                          Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Cosine of final theta.
      dimension cth     (1)
c---- Cosine of initial theta.
      dimension cths    (1)
c---- The distance from "ps" to "p".
      dimension dpmove  (1)
c---- Coordinate r of point "p".
      dimension pr      (1)
c---- Coordinate r of point "ps".
      dimension psr     (1)
c---- Coordinate z of point "ps".
      dimension psz     (1)
c---- Coordinate z of point "p".
      dimension pz      (1)
c---- Sine of final theta.
      dimension sth     (1)
c---- Sine of initial theta.
      dimension sths    (1)
c---- Component r of unit vector "u".
      dimension ur      (1)
c---- Component r of unit vector "us".
      dimension usr     (1)
c---- Theta component of unit vector "us".
      dimension ust     (1)
c---- Component z of unit vector "us".
      dimension usz     (1)
c---- Theta component of unit vector "u".
      dimension ut      (1)
c---- Component z of unit vector "u".
      dimension uz      (1)

c.... Local variables.

c---- Cosine of angle theta.
      common /laptmvcy/ costh

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

c---- Index in arrays.
      common /laptmvcy/ n
c---- First index of subset of data.
      common /laptmvcy/ n1
c---- Last index of subset of data.
      common /laptmvcy/ n2
c---- Size of current subset of data.
      common /laptmvcy/ ns
c---- Coordinate x at point "ps".
      common /laptmvcy/ px
c---- Estimated error in px.
      common /laptmvcy/ pxerr
c---- Coordinate y at point "ps".
      common /laptmvcy/ py
c---- Estimated error in py.
      common /laptmvcy/ pyerr
c---- Estimated error in pz.
      common /laptmvcy/ pzerr
c---- Sine of angle theta.
      common /laptmvcy/ sinth
c---- Estimated error in ur.
      common /laptmvcy/ urerr
c---- Estimated error in ut.
      common /laptmvcy/ uterr
c---- Component x of unit vector "us".
      common /laptmvcy/ ux
c---- Estimated error in ux.
      common /laptmvcy/ uxerr
c---- Component y of unit vector "us".
      common /laptmvcy/ uy
c---- Estimated error in uy.
      common /laptmvcy/ uyerr
c---- Magnitude of new direction vector.
      common /laptmvcy/ vlen    (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptmvcy initial point, direction and distance',
cbug     &  ' (kth=',i2,'):' /
cbug     &  (i3,' pz,pr=   ',1p2e22.14 /
cbug     &  '    uz,ur,ut=',1p3e22.14 /
cbug     &  '    dpmove=  ',1pe22.14))
cbug 9902 format (i3,' cth,sth= ',1p2e22.14)
cbug      write ( 3, 9901) kth, (n, psz(n), psr(n),
cbug     &  usz(n), usr(n), ust(n), dpmove(n), n = 1, np)
cbug      if (kth .eq. 1) then
cbug        write ( 3, 9902) (n, cths(n), sths(n), n = 1, np)
cbug      endif
cbugc***DEBUG ends.

c.... Initialize.

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

      nerr = 0

c.... Test for input errors.

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

      if ((kth .lt. 0) .or. (kth .gt. 1)) then
        nerr = 2
        go to 210
      endif

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

      n1 = 1
      n2 = min (np, 64)

c.... Loop over subsets of data.

  110 ns = n2 - n1 + 1

c.... Find new point "p" and new unit direction vector "u".

c---- No truncation.
      if (tol .le. 0.0) then

c---- No calculation of cth, sth.
        if (kth .eq. 0) then

c---- Loop over subset of data.
          do 120 n = n1, n2

            pz(n) = psz(n) + dpmove(n) * usz(n)

            px    = psr(n) + dpmove(n) * usr(n)
            py    = dpmove(n) * ust(n)
            pr(n) = sqrt (px**2 + py**2)

            uz(n) = usz(n)
            ur(n) = ( usr(n) * (px + fuz) + ust(n) * py) /
     &              (pr(n) + fuz)
            ut(n) = (-usr(n) * py + ust(n) * (px + fuz)) /
     &              (pr(n) + fuz)

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

c---- Calculate cth, sth.
        else

c---- Loop over subset of data.
          do 130 n = n1, n2

            pz(n)  = psz(n) + dpmove(n) * usz(n)

            ux     = usr(n) * cths(n) - ust(n) * sths(n)
            uy     = usr(n) * sths(n) + ust(n) * cths(n)

            px     = psr(n) * cths(n) + dpmove(n) * ux
            py     = psr(n) * sths(n) + dpmove(n) * uy
            pr(n)  = sqrt (px**2 + py**2)

            cth(n) = (px + fuz) / (pr(n) + fuz)
            sth(n) = py / (pr(n) + fuz)

            uz(n)  = usz(n)
            ur(n)  =  ux * cth(n) + uy * sth(n)
            ut(n)  = -ux * sth(n) + uy * cth(n)

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

c---- Tested kth.
        endif

c---- Truncate small values to zero.
      else

c---- No calculation of cth, sth.
        if (kth .eq. 0) then

c---- Loop over subset of data.
          do 140 n = n1, n2

            pz(n)  = psz(n) + dpmove(n) * usz(n)
            pzerr  = tol * (abs (psz(n)) + abs (dpmove(n) * usz(n)))
            if (abs (pz(n)) .lt. pzerr) then
              pz(n) = 0.0
            endif

            px     = psr(n) + dpmove(n) * usr(n)
            pxerr  = tol * (abs (psr(n)) + abs (dpmove(n) * usr(n)))
            if (abs (px) .lt. pxerr) then
              px = 0.0
            endif

            py     = dpmove(n) * ust(n)

            pr(n)  = sqrt (px**2 + py**2)

            costh  = (px + fuz) / (pr(n) + fuz)
            if (abs (costh) .lt. tol) then
              costh = 0.0
            endif
            sinth  = py / (pr(n) + fuz)
            if (abs (sinth) .lt. tol) then
              sinth = 0.0
            endif

            uz(n)  = usz(n)

            ur(n)  = usr(n) * costh + ust(n) * sinth
            urerr  = tol * (abs (usr(n) * costh) +
     &                      abs (ust(n) * sinth))
            if (abs (ur(n)) .lt. urerr) then
              ur(n) = 0.0
            endif

            ut(n)  = -usr(n) * sinth + ust(n) * costh
            uterr  = tol * (abs (usr(n) * sinth) +
     &                      abs (ust(n) * costh))
            if (abs (ut(n)) .lt. uterr) then
              ut(n) = 0.0
            endif

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

c---- Calculate cth, sth.
        else

c---- Loop over subset of data.
          do 150 n = n1, n2

            pz(n)  = psz(n) + dpmove(n) * usz(n)
            pzerr  = tol * (abs (psz(n)) + abs (dpmove(n) * usz(n)))
            if (abs (pz(n)) .lt. pzerr) then
              pz(n) = 0.0
            endif

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

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

            ux     = usr(n) * costh - ust(n) * sinth
            uxerr  = tol * (abs (usr(n) * costh) +
     &                      abs (ust(n) * sinth))
            if (abs (ux) .lt. uxerr) then
              ux = 0.0
            endif

            uy     = usr(n) * sinth + ust(n) * costh
            uyerr  = tol * (abs (usr(n) * sinth) +
     &                      abs (ust(n) * costh))
            if (abs (uy) .lt. uyerr) then
              uy = 0.0
            endif

            px     = psr(n) * costh + dpmove(n) * ux
            pxerr  = tol * (abs (psr(n) * costh) +
     &                      abs (dpmove(n) * ux))
            if (abs (px) .lt. pxerr) then
              px = 0.0
            endif

            py     = psr(n) * sinth + dpmove(n) * uy
            pyerr  = tol * (abs (psr(n) * sinth) +
     &                      abs (dpmove(n) * uy))
            if (abs (py) .lt. pyerr) then
              py = 0.0
            endif

            pr(n)  = sqrt (px**2 + py**2)

            uz(n)  = usz(n)

            costh  = (px + fuz) / (pr(n) + fuz)
            if (abs (costh) .lt. tol) then
              costh = 0.0
            endif

            sinth  = py / (pr(n) + fuz)
            if (abs (sinth) .lt. tol) then
              sinth = 0.0
            endif

            cth(n) = costh
            sth(n) = sinth

            ur(n)  =  ux * costh + uy * sinth
            urerr  = tol * (abs (ux * costh) + abs (uy * sinth))
            if (abs (ur(n)) .lt. urerr) then
              ur(n) = 0.0
            endif

            ut(n)  = -ux * sinth + uy * costh
            uterr  = tol * (abs (ux * sinth) + abs (uy * costh))
            if (abs (ut(n)) .lt. uterr) then
              ut(n) = 0.0
            endif

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

c---- Tested kth.
        endif

c---- Tested tol.
      endif

c.... Renormalize unit direction vector.

      call aptvuna (ur(n1), ut(n1), uz(n1), ns, 0., vlen, nerr)

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 9903 format (/ 'aptmvcy results:  new point and direction:' /
cbug     &  (i3,' pz,pr=   ',1p2e22.14 /
cbug     &  '    uz,ur,ut=',1p3e22.14))
cbug      write ( 3, 9903) (n, pz(n), pr(n),
cbug     &  uz(n), ur(n), ut(n), n = 1, np)
cbug      if (kth .eq. 1) then
cbug        write ( 3, 9902) (n, cths(n), sths(n), n = 1, np)
cbug      endif
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832