subroutine aptscnv (frad, ra, ax, ay, az, rb, bx, by, bz,
     &                    np, tol, cx, cy, cz, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTSCNV
c
c     call aptscnv (frad, ra, ax, ay, az, rb, bx, by, bz, np, tol,
c    &              cx, cy, cz, nerr)
c
c     Version:  aptscnv  Updated    1991 June 6 13:40.
c               aptscnv  Originated 1991 June 6 13:40.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To randomly sample np points c = (cx, cy, cz) from a uniform
c               spatial distribution in the volume bounded by the surfaces of
c               two right circular cones with the same vertex and the same axis,
c               and by two parallel planes perpendicular to that axis.  The two
c               planar boundaries are the annular disk centered at point
c               a = (ax, ay, az) with inner and outer radii frad * ra and ra,
c               and the annular disk centered at point b = (bx, by, bz) with
c               inner and outer radii frad * rb and rb.  If frad = 1.0, the
c               points "c" will be sampled uniformly over the conical surface
c               between the two circles.
c               Coordinates of "c" 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     Notes:    The sampled volume is equal to
c                 (1.0 - frad**2) * pi * (ra**2 + rb**2 + ra * rb) * h / 3,
c               and the area of the outer conical surface between the disks is
c                 pi * (ra + rb) * sqrt ((rb - ra)**2 + h**2),
c               where pi = 3.14159265358979..., and h is the distance between
c               points "a" and "b", equal to
c                 sqrt ((bx - ax)**2 + (by - ay)**2 + (bz - az)**2).
c               The vertex of the cones is at point v = (vx, vy, vz), given by
c                 v = (rb * a - ra * b) / (rb - ra),
c               where v, a and b are vectors.
c               The vertex half-angles of the inner and outer conical surfaces
c               are equal to
c                 arctan (frad * (rb - ra) / h) and
c                 arctan ((rb - ra) / h).
c
c     Input:    frad, ra, ax, ay, az, rb, bx, by, bz, np, tol.
c
c     Output:   cx, cy, cz, nerr.
c
c     Calls: aptlocd, aptscad, aptslid, aptspod, aptvdis 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z coordinates of point "a" on the axis of
c                          the two conical surfaces.
c
c     bx,by,bz  Input    The x, y, z coordinates of point "b" on the axis of
c                          the two conical surfaces.
c
c     cx,cy,cz  Output   The coordinates of spatial points "c" sampled randomly
c                          from a uniform distribution over the volume bounded
c                          by two conical surfaces and two planes.  Size np.
c                          Values less than the estimated error in their
c                          calculation, based on tol, will be truncated to zero.
c
c     frad      Input    The ratio of the inner radius to the outer radius of
c                          each of the two annular disks centered at points
c                          "a" and "b".  0.0 .le. frad .le. 1.0.
c                          Use 0.0 to randomly sample the entire volume, and
c                          use 1.0 to randomly sample only on the outer surface.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if ra or rb is negative.
c                          3 if the cone is too short (point "a" is too
c                            close to point "b").
c                          4 if frad .lt. 0.0 or frad .gt. 1.0.
c
c     np        Input    Size of arrays cx, cy, cz.
c                          Number of points "c" to sample.
c
c     ra        Input    The outer radius of the annular disk centered at point
c                          "a", and perpendicular to the axis "ab".
c                          The inner radius is frad * ra.
c
c     rb        Input    The outer radius of the annular disk centered at point
c                          "b", and perpendicular to the axis "ab".
c                          The inner radius is frad * rb.
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 "c".
      dimension cx      (1)
c---- Coordinate y of point "c".
      dimension cy      (1)
c---- Coordinate z of point "c".
      dimension cz      (1)

c.... Local variables.

c---- Component x of vector "ab".
      common /laptscnv/ abx
c---- Component y of vector "ab".
      common /laptscnv/ aby
c---- Component z of vector "ab".
      common /laptscnv/ abz
c---- Estimated error in cx.
      common /laptscnv/ cxerr
c---- Estimated error in cy.
      common /laptscnv/ cyerr
c---- Estimated error in cz.
      common /laptscnv/ czerr
c---- Distance between points "a" and "b".
      common /laptscnv/ dab
c---- Difference 1.0 - frad.
      common /laptscnv/ df
c---- Estimated error in df.
      common /laptscnv/ dferr
c---- Difference rb - ra.
      common /laptscnv/ dr
c---- Estimated error in dr.
      common /laptscnv/ drerr
c---- Component x of vector "d".
      common /laptscnv/ dx      (64)
c---- Component y of vector "d".
      common /laptscnv/ dy      (64)
c---- Component z of vector "d".
      common /laptscnv/ dz      (64)
c---- Randomly sampled fractional distance.
      common /laptscnv/ fab     (64)

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

c---- Index in local array.
      common /laptscnv/ n
c---- First index of subset of data.
      common /laptscnv/ n1
c---- Last index of subset of data.
      common /laptscnv/ n2
c---- Index in external array.
      common /laptscnv/ nn
c---- Size of current subset of data.
      common /laptscnv/ ns
c---- Relative probability of ra.
      common /laptscnv/ pa
c---- Relative probability of rb.
      common /laptscnv/ pb
c---- Radius of cone at point c.
      common /laptscnv/ rc      (64)
cbugc***DEBUG begins.
cbug      common /laptscnv/ avgr, avgx, avgy, avgz, devr, devx, devy, devz
cbug      common /laptscnv/ avgr2, avgx2, avgy2, avgz2
cbug 9901 format (/ 'aptscnv sampling in a conical annulus.  np=',i3 /
cbug     &  '  frad    =',1pe22.14 /
cbug     &  '  ra,rb   =',1p2e22.14 /
cbug     &  '  ax,ay,az=',1p3e22.14 /
cbug     &  '  bx,by,bz=',1p3e22.14)
cbug      write (3, 9901) np, frad, ra, rb, ax, ay, az, bx, by, bz
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 ((ra .lt. 0.0) .or. (rb .lt. 0.0)) then
        nerr = 2
        go to 210
      endif

c.... Find vector "ab" between points "a" and "b".

      call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
     &              abx, aby, abz, dab, nerr)

      if (dab .le. tol) then
        nerr = 3
        go to 210
      endif

      if ((frad .lt. 0.0) .or. (frad .gt. 1.0)) then
        nerr = 4
        go to 210
      endif

c.... Find the nature of the sampling space.

      dr    = rb - ra
      drerr = tol * (ra + rb)

c---- Shape is a cylinder or a line.
      if (abs (dr) .le. drerr) then
        dr = 0.0
      endif

      df    = 1.0 - frad
      dferr = tol * frad

c---- Shape is a surface or a line.
      if (abs (df) .le. dferr) then
        df = 0.0
      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.... Find the fractional distances of points "c" along the axis "ab".

c---- Shape is a cylinder or a line.
      if (dr .eq. 0.0) then

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

c---- Uniform distribution.
          fab(n) = ranf( )

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

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

          rc(n)  = ra + fab(n) * dr

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

c---- Sample on a surface.
      elseif (df .eq. 0.0) then

        pa = ra + fuz
        pb = rb + fuz

c++++ Linear distribution.
        call aptslid (0.0, pa, 1.0, pb, ns, fab, nerr)

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

          rc(n)  = frad * (ra + fab(n) * dr)

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

c---- Sample in volume, ra > rb.
      elseif (dr .lt. 0.0) then

c++++ Quadratic distribution.
        call aptspod (2.0, rb, ra, ns, rc, nerr)

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

          fab(n) = (rc(n) - ra) / dr

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

c---- Sample in volume, rb > ra.
      else

c++++ Quadratic distribution.
        call aptspod (2.0, ra, rb, ns, rc, nerr)

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

          fab(n) = (rc(n) - ra) / dr

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

c---- Tested dr for shape of cone.
      endif

c.... Randomly sample vectors perpendicular to the axis of the cone.

c---- On outer surface.
      if (df .eq. 0.0) then

        call aptscad (ns, abx, aby, abz, dx, dy, dz, nerr)

c---- On annular disk.
      else

        call aptlocd (frad, 1.0, 0., 0., 0., abx, aby, abz, ns, tol,
     &                dx, dy, dz, nerr)

      endif

c....                   _   _         __        _
c.... Find points "c".  c = a + fab * ab + rc * d.

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

        nn     = n + n1 - 1
        cx(nn) = ax + fab(n) * abx + rc(n) * dx(n)
        cy(nn) = ay + fab(n) * aby + rc(n) * dy(n)
        cz(nn) = az + fab(n) * abz + rc(n) * dz(n)

c---- End of loop over subset of data.
  170 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 180 n = 1, ns

          nn     = n + n1 - 1

          cxerr = tol * (abs (ax) + fab(n) * (abs (ax) + abs (bx)) +
     &                   2.0 * (ra + fab(n) * (ra + rb)) * abs (dx(n)))

          cyerr = tol * (abs (ay) + fab(n) * (abs (ay) + abs (by)) +
     &                   2.0 * (ra + fab(n) * (ra + rb)) * abs (dy(n)))

          czerr = tol * (abs (az) + fab(n) * (abs (az) + abs (bz)) +
     &                   2.0 * (ra + fab(n) * (ra + rb)) * abs (dz(n)))


          if (abs(cx(nn)) .lt. cxerr) then
            cx(nn) = 0.0
          endif
          if (abs(cy(nn)) .lt. cyerr) then
            cy(nn) = 0.0
          endif
          if (abs(cz(nn)) .lt. czerr) then
            cz(nn) = 0.0
          endif

c---- End of loop over subset of data.
  180   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 (/ 'aptscnv results:')
cbug 9903 format (i4,' cx,cy,cz=',1p3e22.14)
cbug 9904 format (/ 'aptscnv 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 190 n = 1, np
cbug        avgr   = avgr + sqrt (cx(n)**2 + cy(n)**2 + cz(n)**2)
cbug        avgx   = avgx + cx(n)
cbug        avgy   = avgy + cy(n)
cbug        avgz   = avgz + cz(n)
cbug        avgr2  = avgr2 + cx(n)**2 + cy(n)**2 + cz(n)**2
cbug        avgx2  = avgx2 + cx(n)**2
cbug        avgy2  = avgy2 + cy(n)**2
cbug        avgz2  = avgz2 + cz(n)**2
cbug        write ( 3, 9903) n, cx(n), cy(n), cz(n)
cbug  190 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 (amax1 (0.0, avgr2 - avgr**2))
cbug      devx  = sqrt (amax1 (0.0, avgx2 - avgx**2))
cbug      devy  = sqrt (amax1 (0.0, avgy2 - avgy**2))
cbug      devz  = sqrt (amax1 (0.0, 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 aptscnv.      (+1 line.)
      end

UCRL-WEB-209832