subroutine aptqprr (ax, ay, az, bx, by, bz, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, nran, tol,
     &                    nprox, cx, cy, cz, dist, dx, dy, dz,
     &                    snx, sny, snz, costh, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTQPRR
c
c     call aptqprr (ax, ay, az, bx, by, bz, qc, qx, qy, qz,
c                   qxy, qyz, qzx, qxx, qyy, qzz, nran, tol,
c                   nprox, cx, cy, cz, dist, dx, dy, dz,
c                   snx, sny, snz, costh, nerr)
c
c     Version:  aptqprr  Updated    1994 December 20 12:50.
c               aptqprr  Originated 1994 October 31 13:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To try to find the proximal point "c" = (cx, cy, cz) on a
c               quadric surface, nearest point "a" = (ax, ay, az), with an
c               initial trial direction "b" = (bx, by, bz).
c
c               Note:  use APT subroutine aptwhis to find the proximal point.
c               Use this subroutine only if aptwhis fails.
c
c               The quadric surface is represented by:
c
c                 f(x,y,z) = qc + qx  * x     + qy  * y     + qz  * z     +
c                                 qxy * x * y + qyz * y * z + qzx * z * x +
c                                 qxx * x**2  + qyy * y**2  + qzz * z**2 = 0,
c
c               The method is to try various directions from point "a", starting
c               with direction "b", and the directions to any extrema and axial
c               intercept points on the quadric surface, and the x, y and z
c               directions, then trying nran directions randomly sampled from
c               an isotropic distribution, and then trying up to 10 * nran
c               (or until the result converges, based on tol) directions
c               randomly sampled from a cosine-power distribution directed
c               toward the last best direction, with the power geometrically
c               increasing from 1.0 to 1.e10.
c
c               The final vector from point "a" to point "c" is vector
c               "d" = (dx, dy, dz) = (cx - ax, cy - ay, cz - az).
c
c               The vector normal to the quadric surface at point "c" is
c               "sn" = (snx, sny, snz) = grad f:
c                 snx  = qx + 2.0 * qxx * cx + qxy * cy + qzx * cz,
c                 sny  = qy + qxy * cx + 2.0 * qyy * cy + qyz * cz,
c                 snz  = qz + qzx * cx + qyz * cy + 2.0 * qzz * cz,
c               and the cosine of the angle between vector "d" and the normal
c               vector "sn" is costh, which is -1 or 1 if the method converges
c               (nprox = 1).
c
c               Any result less than the estimated error in its calculation,
c               based on tol, will be truncated to zero.
c
c               Flag nerr indicates any input error.
c
c     Input:    ax, ay, az, bx, by, bz, qc, qx, qy, qz, qxy, qyz, qzx,
c               qxx, qyy, qzz, nran, tol.
c
c     Output:   nprox, cx, cy, cz, dist, dx, dy, dz, snx, sny, snz, costh, nerr.
c
c     Calls: aptrkis, aptscat, aptscav, aptvusz 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z coordinates of point "a".
c
c     bx,by,bz  Input    The x, y, z components of the first trial vector "b".
c
c     costh     Output   The cosine of the angle between vector "d" (from point
c                          "a" to point "c") and the vector normal to the
c                          quadric surface at point "c".  If nprox = 1, costh
c                          should be either -1.0 or 1.0.
c
c     cx,cy,cz  Output   The x, y, z coordinates of point "c" on the quadric
c                          surface, found nearest to point "a".  If nprox = 1,
c                          the proximal point, to within tolerance limit tol.
c
c     dx,dy,dz  Output   The x, y, z components of the vector from point "a" to
c                          point "c".
c
c     dist      Output   The distance from point "a" to point "c".
c                          If nprox = 1, the proximal distance, to within
c                          tolerance limit tol.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if vector "b" is too small.
c                          2 if nran is not positive.
c
c     nprox     Output   The number of exact proximal points "c" and proximal
c                          distances dist found.
c                         -1 if none, and none of the trial directions from
c                            point "a" intersect the quadric surface.
c                            The coefficients qc, ..., qzz may be bad.
c                          0 if none, but the surface point "c" in direction "d"
c                            from point "a" produced the smallest distance dist
c                            of any directions tried.
c                          1 if an exact proximal point "c" was found, within
c                            tolerance limit tol.
c
c     q..       Input    Coefficients of the implicit equation of a quadric
c                          surface in xyz space (qc, qx, qy, qz, qxy, qyz, qzx,
c                          qxx, qyy, qzz).
c
c     snx,...   Output   The x, y, z components of the vector normal to the
c                          quadric surface at point "c".
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Local variables.

      common /laptqprr/ ccx(4)        ! X coordinate of surface point.
      common /laptqprr/ ccy(4)        ! Y coordinate of surface point.
      common /laptqprr/ ccz(4)        ! Z coordinate of surface point.
      common /laptqprr/ ddist(4)      ! Distance to quadric surface.
      common /laptqprr/ aqp(25)       ! Type of surface point.
      character*8       aqp
      common /laptqprr/ qpx(25)       ! Coordinate x of surface point.
      common /laptqprr/ qpy(25)       ! Coordinate y of surface point.
      common /laptqprr/ qpz(25)       ! Coordinate z of surface point.

cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqprr finding (nran=',i8',) the distance from',
cbug     &  ' point "a"' /
cbug     &  ' to the quadric surface with the coefficients qc-qzz,' /
cbug     &  ' and initial trial direction "b":' /
cbug     &  ' ax, ay, az =',1p3e22.14 /
cbug     &  ' bx, by, bz =',1p3e22.14 /
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)  nran, ax, ay, az, bx, by, bz,
cbug     &  qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.

c.... Initialize.

      nerr    = 0
      nprox   = 0
      ndecr   = 0                     ! Total number of improved guesses.
      distmin =  1.e99                ! Minimum distance found.
      dist    = -1.e99
      cx      = -1.e99
      cy      = -1.e99
      cz      = -1.e99
      dx      = -1.e99
      dy      = -1.e99
      dz      = -1.e99
      snx     = -1.e99
      sny     = -1.e99
      snz     = -1.e99
      costh   = -1.e99

c.... Test for errors.

      call aptvusz (bx, by, bz, tol, ubx, uby, ubz, blen)

      if (blen .le. tol) then
        nerr  = 1
        nprox = -1
        go to 710
      endif

      if (nran .lt. 1) then
        nerr  = 2
        nprox = -1
        go to 710
      endif

c=======================================================================********

c.... Find any extrema and axis intercepts on the quadric surface.

      call aptqupt (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
     &              tol, aqp, qpx, qpy, qpz, nqp, nerr)

c...  Add specified direction "b", and x, y and z directions.

      nqp      = nqp + 1
      qpx(nqp) = ax + bx
      qpy(nqp) = ay + by
      qpz(nqp) = az + bz
      nqp      = nqp + 1
      qpx(nqp) = ax + 1.0
      qpy(nqp) = ay + 0.0
      qpz(nqp) = az + 0.0
      nqp      = nqp + 1
      qpx(nqp) = ax + 0.0
      qpy(nqp) = ay + 1.0
      qpz(nqp) = az + 0.0
      nqp      = nqp + 1
      qpx(nqp) = ax + 0.0
      qpy(nqp) = ay + 0.0
      qpz(nqp) = az + 1.0

c=======================================================================********

c.... Loop over input vector and directions from point "a" to surface points,
c....   then loop over nran random directions.

      nranp   = nran + nqp

      do 210 nd = 1, nranp            ! Loop over random directions.

        if (nd .le. nqp) then       ! Use directions to surface points, and "b".

          call aptvdis (ax, ay, az, qpx(nd), qpy(nd), qpz(nd), 1, tol,
     &                  vx, vy, vz, vlen, nerra)

          if (vlen .eq. 0.0) go to 210

          vx = vx / vlen
          vy = vy / vlen
          vz = vz / vlen

        else                          ! Use a random direction.

          call aptscat (1, vx, vy, vz, nerra)

        endif

c....   Find the intersection of the track from point "a" in direction "v".

        dintmn = -distmin
        dintmx =  distmin

        call aptrkis (ax, ay, az, vx, vy, vz,
     &                qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
     &                dintmn, dintmx, 1, tol,
     &                nint, sx, sy, sz, distest,
     &                ssnx, ssny, ssnz, cosn, nerra)

        if (nint .eq. 0) go to 210

        if (abs (distest) .le. distmin) then

          if (distest .lt. 0.0) then
            distest = -distest
            cosn    = -cosn
            vx      = -vx
            vy      = -vy
            vz      = -vz
          endif

          ndecr   = ndecr + 1
          cx      = sx
          cy      = sy
          cz      = sz
          distmin = abs (distest)
          dist    = distest
          dx      = vx * distest
          dy      = vy * distest
          dz      = vz * distest
          snx     = ssnx
          sny     = ssny
          snz     = ssnz
          costh   = cosn
cbugc***DEBUG begins.
cbug 9902   format (/ 'Proximal x,y,z   ?',1p3e20.12 /
cbug     &            'Proximal dx,dy,dz?',1p3e20.12 /
cbug     &            '  Approx distance?',1pe20.12,'  Test',i8 )
cbug 9903   format (  '  Cosine of angle ',1pe20.12,'  (from normal)',
cbug     &          5x,1pe20.12 )
cbug        power = 0.0
cbug        write ( 3, 9902) cx, cy, cz, dx, dy, dz, dist, ndecr
cbug        write ( 3, 9903) costh, power
cbugc***DEBUG ends.

          if (distmin .eq. 0.0) then
            nprox = 1
            costh = 1.0
            go to 710
          endif

        endif

  210 continue                        ! End of loop over random directions.

      if (ndecr .eq. 0) then          ! No intersections.
        nprox = -1
        go to 710
      endif

      delcos  = 1.0 - abs (costh)

      if (delcos  .le. tol) then      ! Found an exact proximal point.
        nprox = 1
        go to 710
      endif

c=======================================================================********

c.... Search in cosine-power distribution around the best direction.

      numtot  = 10 * nran
      power   = 1.0
      powmult = 10.0**(1.0 / nran)

      do 250 nd = 1, numtot           ! Loop over cosine-power distribution.

        call aptvusz (dx, dy, dz, tol,
     &                ugoodx, ugoody, ugoodz, dlen)

        call aptscav (1, power, ugoodx, ugoody, ugoodz,
     &                vx, vy, vz, nerra)

        power = power * powmult

c....   Find the intersection of the track from point P in direction V.

        dintmn = -distmin
        dintmx =  distmin

        call aptrkis (ax, ay, az, vx, vy, vz,
     &                qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
     &                dintmn, dintmx, 1, tol,
     &                nint, sx, sy, sz, distest,
     &                ssnx, ssny, ssnz, cosn, nerra)

        if (nint .eq. 0) go to 250    ! No intersection.

        if (abs (distest) .le. distmin) then

          if (distest .lt. 0.0) then
            distest = -distest
            cosn    = -cosn
            vx      = -vx
            vy      = -vy
            vz      = -vz
          endif

          ndecr   = ndecr + 1
          cx      = sx
          cy      = sy
          cz      = sz
          distmin = abs (distest)
          dist    = distest
          dx      = vx * distest
          dy      = vy * distest
          dz      = vz * distest
          snx     = ssnx
          sny     = ssny
          snz     = ssnz
          costh   = cosn
cbugc***DEBUG begins.
cbug        write ( 3, 9902) cx, cy, cz, dx, dy, dz, dist, ndecr
cbug        write ( 3, 9903) costh, power
cbugc***DEBUG ends.

          deldist = abs (abs (distest) - distmin)
          delcos  = 1.0 - abs (costh)

c....     Test for final convergence.

          if ((deldist .le. tol * distmin) .and.
     &        (delcos  .le. tol))          then

            nprox = 1
            go to 710

          endif

        endif                         ! Improved guess.

  250 continue                        ! End of loop over cosine-power distr.

c=======================================================================********

  710 continue

cbugc***DEBUG begins.
cbug 9904 format (/ 'aptqprr results:  nprox=',i2,' ndecr=',i8,' nerr=',i3 /
cbug     &  ' dist       =',1pe22.14 /
cbug     &  ' cx, cy, cz =',1p3e22.14 /
cbug     &  ' dx, dy, dz =',1p3e22.14 /
cbug     &  ' snx,sny,snz=',1p3e22.14 /
cbug     &  ' costh      =',1pe22.14 )
cbug      write ( 3, 9904)  nprox, ndecr, nerr, dist, cx, cy, cz,
cbug     &  dx, dy, dz, snx, sny, snz, costh
cbugc***DEBUG ends.

      return

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

UCRL-WEB-209832