subroutine aptetrp (ax, ay, az, bx, by, bz, cx, cy, cz,
     &                    dx, dy, dz, np, kloc, tol,
     &                    px, py, pz, wa, wb, wc, wd, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTETRP
c
c     call aptetrp (ax, ay, az, bx, by, bz, cx, cy, cz,
c    &              dx, dy, dz, np, kloc, tol,
c    &              px, py, pz, wa, wb, wc, wd, nerr)
c
c     Version:  aptetrp  Updated    1990 November 27 14:00.
c               aptetrp  Originated 1990 May 9 9:30.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To randomly sample np points p = (px, py, pz) from a uniform
c               spatial distribution inside the tetrahedron with vertices
c               a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz) and
c               d = (dx, dy, dz) in 3-D space, and optionally (kloc = 1),
c               to find the local coordinates or vertex weights wa, wb,
c               wc, wd, corresponding to the point "p" in the tetrahedron:
c
c                 p = wa * a + wb * b + wc * c + wd * d
c
c               Coordinates of "p" 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     Input:    ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np,
c               kloc, tol.
c
c     Output:   px, py, pz, wa, wb, wc, wd, nerr.
c
c     Calls: aptspod, apttlod 
c               
c
c     Glossary:
c
c     ax,ay,az  Input    The x, y, z coordinates of tetrahedron vertex "a".
c
c     bx,by,bz  Input    The x, y, z coordinates of tetrahedron vertex "b".
c
c     cx,cy,cz  Input    The x, y, z coordinates of tetrahedron vertex "c".
c
c     dx,dy,dz  Input    The x, y, z coordinates of tetrahedron vertex "d".
c
c     kloc      Input    Option flag to return local tetrahedron coordinates
c                          for each point "p", if 1.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if kloc is not 0 or 1.
c
c     np        Input    Size of arrays px, py, pz.
c                          Number of points "p" to sample.
c
c     px,py,pz  Output   The coordinates of spatial points "p" randomly sampled
c                          from a uniform distribution over the volume of the
c                          tetrahedron.  Size np.  Values less than the
c                          estimated error in their calculation, based on tol,
c                          will be truncated to zero.
c
c     tol       Input    Numerical tolerance limit.
c                          On Cray computers, recommend 1.e-5 to 1.e-11.
c
c     wa, ...   Output   Fractional distances of point "p" from the opposite
c                          triangular faces to the tetrahedron vertices "a", ...
c                          Size np only if kloc = 1.  Otherwise, not returned.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Coordinate x of point "p".
      dimension px      (1)
c---- Coordinate y of point "p".
      dimension py      (1)
c---- Coordinate z of point "p".
      dimension pz      (1)
c---- Fractional distance of "p" to "a".
      dimension wa      (1)
c---- Fractional distance of "p" to "b".
      dimension wb      (1)
c---- Fractional distance of "p" to "c".
      dimension wc      (1)
c---- Fractional distance of "p" to "d".
      dimension wd      (1)

c.... Local variables.

c---- Component x of vector "e".
      common /laptetrp/ ex      (64)
c---- Component y of vector "e".
      common /laptetrp/ ey      (64)
c---- Component z of vector "e".
      common /laptetrp/ ez      (64)
c---- Fractional distance from "d" to "abc".
      common /laptetrp/ fdis    (64)
c---- Index in local array.
      common /laptetrp/ n
c---- First index of subset of data.
      common /laptetrp/ n1
c---- Last index of subset of data.
      common /laptetrp/ n2
c---- Index in external array.
      common /laptetrp/ nn
c---- Size of current subset of data.
      common /laptetrp/ ns
c---- Estimated error in px.
      common /laptetrp/ pxerr
c---- Estimated error in py.
      common /laptetrp/ pyerr
c---- Estimated error in pz.
      common /laptetrp/ pzerr
c---- Weight of triangle vertex "a".
      common /laptetrp/ twa     (64)
c---- Weight of triangle vertex "b".
      common /laptetrp/ twb     (64)
c---- Weight of triangle vertex "c".
      common /laptetrp/ twc     (64)
cbugc***DEBUG begins.
cbug      common /laptetrp/ avgx, avgy, avgz, devx, devy, devz
cbug      common /laptetrp/ avgx2, avgy2, avgz2
cbug      common /laptetrp/ avgwa, avgwb, avgwc, avgwd
cbug      common /laptetrp/ devwa, devwb, devwc, devwd
cbug      common /laptetrp/ avgwa2, avgwb2, avgwc2, avgwd2
cbug 9901 format (/ 'aptetrp sampling inside a tetrahedron.  np=',i3 /
cbug     &  '  ax,ay,az=',1p3e22.14 /
cbug     &  '  bx,by,bz=',1p3e22.14 /
cbug     &  '  cx,cy,cz=',1p3e22.14 /
cbug     &  '  dx,dy,dz=',1p3e22.14)
cbug      write (3, 9901) np, ax, ay, az, bx, by, bz,
cbug     &  cx, cy, cz, dx, dy, dz
cbugc***DEBUG ends.

c.... initialize.

      nerr = 0

c.... Test for input errors.

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

      if ((kloc .lt. 0) .or. (kloc .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 the data subsets.

  110 ns = n2 - n1 + 1

c.... Randomly sample points in the triangle "abc".

      call apttlod (ax, ay, az, bx, by, bz, cx, cy, cz, ns, 1, tol,
     &              ex, ey, ez, twa, twb, twc, nerr)

c.... Randomly sample fractional distances from a quadratic distribution.

      call aptspod (2., 0., 1., ns, fdis, nerr)

c.... Find points "p" randomly sampled uniformly over the tetrahedron.

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

        nn     = n + n1 - 1
        px(nn) = dx + fdis(n) * (ex(n) - dx)
        py(nn) = dy + fdis(n) * (ey(n) - dy)
        pz(nn) = dz + fdis(n) * (ez(n) - dz)

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

          pxerr  = tol * (abs (dx) + fdis(n) * abs (ex(n) - dx))
          pyerr  = tol * (abs (dy) + fdis(n) * abs (ey(n) - dy))
          pzerr  = tol * (abs (dz) + fdis(n) * abs (ez(n) - dz))

          if (abs(px(nn)) .lt. pxerr) then
            px(nn) = 0.0
          endif
          if (abs(py(nn)) .lt. pyerr) then
            py(nn) = 0.0
          endif
          if (abs(pz(nn)) .lt. pzerr) then
            pz(nn) = 0.0
          endif

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

c---- Tested tol.
      endif

c.... See if local coordinates are to be returned.

c---- Local coordinates needed.
      if (kloc .eq. 1) then

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

          nn     = n + n1 - 1

          wa(nn) = fdis(n) * twa(n)
          wb(nn) = fdis(n) * twb(n)
          wc(nn) = fdis(n) * twc(n)
          wd(nn) = 1.0 - fdis(n)

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

c---- Tested kloc.
      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 (/ 'aptetrp results:')
cbug 9903 format (i3,' px,py,pz=',1p3e22.14)
cbug 9904 format (/ 'aptetrp mean and deviation:' /
cbug     &  '  avg(x,y,z)=',1p3e22.14 /
cbug     &  '  dev(x,y,z)=',1p3e22.14)
cbug 9905 format (i3,' wa,wb,wc,wd=',4f15.12)
cbug 9906 format (/ 'aptetrp mean and deviation:' /
cbug     &  '  avg(wa,wb,wc,wd)=',4f15.12 /
cbug     &  '  dev(wa,wb,wc,wd)=',4f15.12)
cbug
cbug      write ( 3, 9902)
cbug
cbug      avgx  = 0.0
cbug      avgy  = 0.0
cbug      avgz  = 0.0
cbug      avgx2 = 0.0
cbug      avgy2 = 0.0
cbug      avgz2 = 0.0
cbug
cbug      do 150 n = 1, np
cbug        avgx   = avgx  + px(n)
cbug        avgy   = avgy  + py(n)
cbug        avgz   = avgz  + pz(n)
cbug        avgx2  = avgx2 + px(n)**2
cbug        avgy2  = avgy2 + py(n)**2
cbug        avgz2  = avgz2 + pz(n)**2
cbug        write ( 3, 9903) n, px(n), py(n), pz(n)
cbug  150 continue
cbug
cbug      avgx  = avgx  / np
cbug      avgy  = avgy  / np
cbug      avgz  = avgz  / np
cbug      avgx2 = avgx2 / np
cbug      avgy2 = avgy2 / np
cbug      avgz2 = avgz2 / np
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
cbug
cbugc---- Local coordinates returned.
cbug      if (kloc .eq. 1) then
cbug
cbug        avgwa  = 0.0
cbug        avgwb  = 0.0
cbug        avgwc  = 0.0
cbug        avgwd  = 0.0
cbug        avgwa2 = 0.0
cbug        avgwb2 = 0.0
cbug        avgwc2 = 0.0
cbug        avgwd2 = 0.0
cbug
cbug        do 160 n = 1, np
cbug          avgwa  = avgwa  + wa(n)
cbug          avgwb  = avgwb  + wb(n)
cbug          avgwc  = avgwc  + wc(n)
cbug          avgwd  = avgwd  + wd(n)
cbug          avgwa2 = avgwa2 + wa(n)**2
cbug          avgwb2 = avgwb2 + wb(n)**2
cbug          avgwc2 = avgwc2 + wc(n)**2
cbug          avgwd2 = avgwd2 + wd(n)**2
cbug          write ( 3, 9905) n, wa(n), wb(n), wc(n), wd(n)
cbug  160   continue
cbug
cbug        avgwa  = avgwa  / np
cbug        avgwb  = avgwb  / np
cbug        avgwc  = avgwc  / np
cbug        avgwd  = avgwd  / np
cbug        avgwa2 = avgwa2 / np
cbug        avgwb2 = avgwb2 / np
cbug        avgwc2 = avgwc2 / np
cbug        avgwd2 = avgwd2 / np
cbug        devwa  = sqrt (avgwa2 - avgwa**2)
cbug        devwb  = sqrt (avgwb2 - avgwb**2)
cbug        devwc  = sqrt (avgwc2 - avgwc**2)
cbug        devwd  = sqrt (avgwd2 - avgwd**2)
cbug
cbug        write ( 3, 9906) avgwa, avgwb, avgwc, avgwd,
cbug     &    devwa, devwb, devwc, devwd
cbug
cbugc---- Tested kloc.
cbug      endif
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832