subroutine aptlocs (rin, rout, ax, ay, az, np, tol,
     &                    bx, by, bz, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTLOCS
c
c     call aptlocs (rin, rout, ax, ay, az, np, tol,
c    &              bx, by, bz, nerr)
c
c     Version:  aptlocs  Updated    1990 November 29 11:40.
c               aptlocs  Originated 1990 May 4 10:00.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To randomly sample np points b = (bx, by, bz) from a uniform
c               spatial distribution inside the spherical shell with inner
c               radius rin, outer radius rout, centered at point
c               a = (ax, ay, az).
c               Coordinates of "b" less than the estimated error in their
c               calculation, based on tol, will be truncated to zero.
c               Flag nerr indicates any input error.
c
c     History:  1990 May 7 16:00.  Deleted arguments tolx, toly, tolz from call
c               to aptscat.
c
c     Input:    rin, rout, ax, ay, az, np, tol.
c
c     Output:   bx, by, bz, nerr.
c
c     Calls: aptscat, aptspod 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z coordinates of point "a" at the center of
c                          the spherical shell of inner radius rin, outer
c                          radius rout.
c
c     bx,by,bz  Output   The coordinates of spatial points "b" sampled randomly
c                          from a uniform distribution over the volume of the
c                          spherical shell.  Size np.  Values less than the
c                          estimated error in their calculation, based on tol,
c                          will be truncated to zero.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if rin is negative, or rout is less than rin.
c
c     np        Input    Size of arrays bx, by, bz.
c                          Number of points "b" to sample.
c
c     rin       Input    The inner radius of the spherical shell.
c
c     rout      Input    The outer radius of the spherical shell.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Coordinate x of point "b".
      dimension bx      (1)
c---- Coordinate y of point "b".
      dimension by      (1)
c---- Coordinate z of point "b".
      dimension bz      (1)

c.... Local variables.

c---- Estimated error in bx.
      common /laptlocs/ bxerr
c---- Estimated error in by.
      common /laptlocs/ byerr
c---- Estimated error in bz.
      common /laptlocs/ bzerr
c---- Component x of vector "c".
      common /laptlocs/ cx      (64)
c---- Component y of vector "c".
      common /laptlocs/ cy      (64)
c---- Component z of vector "c".
      common /laptlocs/ cz      (64)
c---- Index in local array.
      common /laptlocs/ n
c---- First index of subset of data.
      common /laptlocs/ n1
c---- Last index of subset of data.
      common /laptlocs/ n2
c---- Index in external array.
      common /laptlocs/ nn
c---- Size of current subset of data.
      common /laptlocs/ ns
c---- Randomly sampled radius.
      common /laptlocs/ rsam    (64)
cbugc***DEBUG begins.
cbug      common /laptlocs/ avgx, avgy, avgz, devx, devy, devz
cbug      common /laptlocs/ avgx2, avgy2, avgz2
cbug      common /laptlocs/ avgr, avgr2, devr, rad, rad2
cbug 9901 format (/ 'aptlocs sampling inside a spherical shell.',
cbug     &  '  np=',i3 /
cbug     &  '  rin,rout=',1p2e22.14 /
cbug     &  '  ax,ay,az=',1p3e22.14)
cbug      write (3, 9901) np, rin, rout, ax, ay, az
cbugc***DEBUG ends.

c.... initialize.

      nerr = 0

c.... Test for input errors.

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

      if ((rout .lt. rin) .or. (rin .lt. 0.0)) 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 the data subsets.

  110 ns = n2 - n1 + 1

c.... Randomly sample unit vectors "c" in the radial direction of the sphere.

      call aptscat (ns, cx, cy, cz, nerr)

c.... Randomly sample radii from a quadratic distribution from rin to rout.

      call aptspod (2.0, rin, rout, ns, rsam, nerr)

c.... Find points "b" sampled uniformly over the spherical shell.

c---- Loop over subset of data.
      do 120 n = 1, ns

        nn     = n + n1 - 1
        bx(nn) = ax + rsam(n) * cx(n)
        by(nn) = ay + rsam(n) * cy(n)
        bz(nn) = az + rsam(n) * cz(n)

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

c.... See if truncation to zero is allowed.

c---- Truncation to zero allowed.
      if (tol .gt. 0.0) then

c---- Loop over subset of data.
        do 130 n = 1, ns

          nn     = n + n1 - 1

          bxerr  = tol * (abs (ax) + rsam(n) * abs (cx(n)))
          byerr  = tol * (abs (ay) + rsam(n) * abs (cy(n)))
          bzerr  = tol * (abs (az) + rsam(n) * abs (cz(n)))

          if (abs(bx(nn)) .lt. bxerr) then
            bx(nn) = 0.0
          endif
          if (abs(by(nn)) .lt. byerr) then
            by(nn) = 0.0
          endif
          if (abs(bz(nn)) .lt. bzerr) then
            bz(nn) = 0.0
          endif

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

      endif

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
cbug 9902 format (/ 'aptlocs results:')
cbug 9903 format (i3,' bx,by,bz=',1p3e22.14 /
cbug     &  '    rad=     ',1pe22.14)
cbug 9904 format (/ 'aptlocs mean and deviation:' /
cbug     &  '  avg(x,y,z)=',1p3e22.14 /
cbug     &  '  dev(x,y,z)=',1p3e22.14 /
cbug     &  '  avgr,devr= ',1p2e22.14)
cbug
cbug      write ( 3, 9902)
cbug
cbug      avgr  = 0.0
cbug      avgx  = 0.0
cbug      avgy  = 0.0
cbug      avgz  = 0.0
cbug      avgr2 = 0.0
cbug      avgx2 = 0.0
cbug      avgy2 = 0.0
cbug      avgz2 = 0.0
cbug
cbug      do 180 n = 1, np
cbug        rad2   = (bx(n) - ax)**2 + (by(n) - ay)**2 +
cbug     &           (bz(n) - az)**2
cbug        rad    = sqrt (rad2)
cbug        avgr   = avgr + rad
cbug        avgx   = avgx + bx(n)
cbug        avgy   = avgy + by(n)
cbug        avgz   = avgz + bz(n)
cbug        avgr2  = avgr2 + rad2
cbug        avgx2  = avgx2 + bx(n)**2
cbug        avgy2  = avgy2 + by(n)**2
cbug        avgz2  = avgz2 + bz(n)**2
cbug        write ( 3, 9903) n, bx(n), by(n), bz(n), rad
cbug  180 continue
cbug
cbug      avgr  = avgr / np
cbug      avgx  = avgx / np
cbug      avgy  = avgy / np
cbug      avgz  = avgz / np
cbug      avgr2 = avgr2 / np
cbug      avgx2 = avgx2 / np
cbug      avgy2 = avgy2 / np
cbug      avgz2 = avgz2 / np
cbug      devr  = sqrt (avgr2 - avgr**2)
cbug      devx  = sqrt (avgx2 - avgx**2)
cbug      devy  = sqrt (avgy2 - avgy**2)
cbug      devz  = sqrt (avgz2 - avgz**2)
cbug
cbug      write ( 3, 9904) avgx, avgy, avgz, devx, devy, devz, avgr, devr
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832