subroutine aptqloc (au, av, bu, bv, cu, cv, du, dv, np, tol,
     &                    fdk, fdl, pu, pv, nerr)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTQLOC
c
c     call aptqloc (au, av, bu, bv, cu, cv, du, dv, np, tol,
c    &              fdk, fdl, pu, pv, nerr)
c
c     Version:  aptqloc  Updated    1990 November 29 15:30.
c               aptqloc  Originated 1990 February 8 14:30.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To randomly sample np spatial points p = (pu, pv) from a
c               uniform distribution over the planar quadrangle in the uv plane,
c               with vertices a = (au, av), b = (bu, bv), c = (cu, cv),
c               d = (du, dv), in counterclockwise order around the quadrangle.
c               For a triangle, make du = au, dv = av, or call apttloc.
c               If the quadrangle has a non-positive area, no points will be
c               sampled.
c
c               Variables fdk and fdl are the local coordinates in the
c               quadrangle.  For any point p = (pu, pv) in the quadrangle:
c
c                 pu = au + fdk * (bu - au) + fdl * (du - au) +
c                      fdk * fdl * (au - bu + cu - du)
c                 pv = av + fdk * (bv - av) + fdl * (dv - av) +
c                      fdk * fdl * (av - bv + cv - dv)
c
c               Given p = (pu, pv), fdk and fdl may be found by calling
c               subroutine aptqfdc.
c
c               Flag nerr indicates any input error.
c
c     History:  1990 February 14 17:30.  Added capability to sample from
c               boomerangs and bowties.
c
c     Input:    au, av, bu, bv, cu, cv, du, dv, np, tol.
c
c     Output:   fdk, fdl, pu, pv, nerr.
c
c     Calls: aptqfdc, aptqvac, aptslid, aptsliv, apttloc 
c               
c
c     Glossary:
c
c     au, av    Input    The u and v coordinates of vertex "a" of quadrangle.
c
c     bu, bv    Input    The u and v coordinates of vertex "b" of quadrangle.
c
c     cu, cv    Input    The u and v coordinates of vertex "c" of quadrangle.
c
c     du, dv    Input    The u and v coordinates of vertex "d" of quadrangle.
c
c     fdk       Output   Fractional distance of point "p" between the quadrangle
c                          sides "da" and "bc".  Range 0.0 to 1.0.  Size np.
c                          Range may be less in a boomerang or bowtie.
c
c     fdl       Output   Fractional distance of point "p" between the quadrangle
c                          sides "ab" and "cd".  Range 0.0 to 1.0.  Size np.
c                          Range may be less in a boomerang or bowtie.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c                          2 if quadrangle area is not positive (no sampling).
c
c     np        Input    Size of arrays fdk, fdl, pu, pv.
c                          Number of points "p" to sample.
c
c     pu, pv    Output   The u and v coordinates of a point sampled randomly
c                          from a quadrangle.  Size np.
c
c     tol       Input    Numerical tolerance limit.  With 64-bit floating
c                          point numbers, 1.e-5 to 1.e-11 is recommended.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

c---- Fractional distance from "da" to "bc".
      dimension fdk     (1)
c---- Fractional distance from "ab" to "cd".
      dimension fdl     (1)
c---- Cooodinate u of a sampled point "p".
      dimension pu      (1)
c---- Coordinate v of a sampled point "p".
      dimension pv      (1)

c.... Local variables.

c---- Coordinate u of sampled pt in "abc".
      common /laptqloc/ abcu    (64)
c---- Coordinate v of sampled pt in "abc".
      common /laptqloc/ abcv    (64)
c---- Parallelogram area at vertex "a".
      common /laptqloc/ aw
c---- Coordinate u of sampled pt in "bcd".
      common /laptqloc/ bcdu    (64)
c---- Coordinate v of sampled pt in "bcd".
      common /laptqloc/ bcdv    (64)
c---- Parallelogram area at vertex "b".
      common /laptqloc/ bw
c---- Coordinate u of sampled pt in "cda".
      common /laptqloc/ cdau    (64)
c---- Coordinate v of sampled pt in "cda".
      common /laptqloc/ cdav    (64)
c---- Parallelogram area at vertex "c".
      common /laptqloc/ cw
c---- Coordinate u of sampled pt in "dab".
      common /laptqloc/ dabu    (64)
c---- Coordinate v of sampled pt in "dab".
      common /laptqloc/ dabv    (64)
c---- Parallelogram area at vertex "d".
      common /laptqloc/ dw
c---- Constant 0.0.
      common /laptqloc/ fl0     (64)
c---- Constant 1.0.
      common /laptqloc/ fl1     (64)

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

c---- Index in local array.
      common /laptqloc/ n
c---- First index of subset of data.
      common /laptqloc/ n1
c---- Last index of subset of data.
      common /laptqloc/ n2
c---- 1 if good fdk, fdl found.
      common /laptqloc/ ngood   (64)
c---- Index in external array.
      common /laptqloc/ nn
c---- Size of current subset of data.
      common /laptqloc/ ns
c---- Shape type of quadrangle.
      common /laptqloc/ ntype
c---- Probability of fdk at 0.0.
      common /laptqloc/ pk0
c---- Probability of fdk at 1.0.
      common /laptqloc/ pk1
c---- Probability of fdl at 0.0.
      common /laptqloc/ pl0     (64)
c---- Probability of fdl at 1.0.
      common /laptqloc/ pl1     (64)
c---- Coordinate u of bowtie point.
      common /laptqloc/ qu
c---- Coordinate v of bowtie point.
      common /laptqloc/ qv
c---- Random number (0.0 to 1.0).
      common /laptqloc/ ranfp   (64)
c---- Dummy argument in apttloc call.
      common /laptqloc/ wa
c---- Dummy argument in apttloc call.
      common /laptqloc/ wb
c---- Dummy argument in apttloc call.
      common /laptqloc/ wc
cbugc***DEBUG begins.
cbugc---- Mean value.
cbug      common /laptqloc/ argmean
cbugc---- Standard deviation from mean.
cbug      common /laptqloc/ argdev
cbug 9901 format (/ 'aptqloc sampling from a quadrangle.  np=',i3 /
cbug     &  '  au=',1pe22.14,' av=',1pe22.14 /
cbug     &  '  bu=',1pe22.14,' bv=',1pe22.14 /
cbug     &  '  cu=',1pe22.14,' cv=',1pe22.14 /
cbug     &  '  du=',1pe22.14,' dv=',1pe22.14)
cbug      write (3, 9901) np, au, av, bu, bv, cu, cv, du, dv
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

c.... Find the areas of the equivalent parallelograms at each vertex,
c....   the shape type, and any bowtie intersection point.
c....   Quadrangle area = 0.5 * (aw + cw) = 0.5 * (bw + dw).
c....   Triangle (du = au, dv = av) area = 0.5 * bw = 0.5 * cw.

      call aptqvac (au, av, bu, bv, cu, cv, du, dv, 1, tol, 2,
     &              aw, bw, cw, dw, ntype, qu, qv, nerr)

c.... Reject quadrangles with non-positive area.

      if ((aw + cw) .le. 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)
  110 ns = n2 - n1 + 1

c---- Quadrangle is regular.
      if (ntype .eq. 0) then

c....   Sample from the quadrangle.

c....   Sample ns values of fdk from the (unnormalized) probability distribution
c....     obtained by integrating prob(fdk,fdl) over fdl from 0.0 to 1.0:
c....     prob(fdk) = 0.5 * (aw + dw) + (bw - aw) * fdk.
c....     For a triangle (du = au, dv = av), prob(fdk) = bw * fdk.

        pk0 = 0.5 * (aw + dw)
        pk1 = 0.5 * (bw + cw)

        call aptslid (0.0, pk0, 1.0, pk1, ns, fdk(n1), nerr)

c....   For each value of fdk, sample a value of fdl from the unnormalized
c....     probability distribution:
c....     prob(fdl) = aw + (bw - aw) * fdk + (dw - aw) * fdl
c....     For a triangle (du = au, dv = av), prob(fdl) = bw * fdk.

c---- Loop over local arrays.
        do 120 n = 1, ns

          nn     = n + n1 - 1
          fl0(n) = 0.0
          fl1(n) = 1.0
          pl0(n) = aw + (bw - aw) * fdk(nn)
          pl1(n) = dw + (cw - dw) * fdk(nn)

c---- End of loop over local arrays.
  120   continue

        call aptsliv (fl0, pl0, fl1, pl1, ns, fdl(n1), nerr)

c....   Find the coordinates of the points.

c---- Loop over points.
        do 130 nn = n1, n2

          pu(nn) = au + fdk(nn) * (bu - au) + fdl(nn) * (du - au) +
     &             fdk(nn) * fdl(nn) * (au - bu + cu - du)
          pv(nn) = av + fdk(nn) * (bv - av) + fdl(nn) * (dv - av) +
     &             fdk(nn) * fdl(nn) * (av - bv + cv - dv)

c---- End of loop over points.
  130   continue

      elseif ((ntype .eq. 12) .or.
     &        (ntype .eq. 14)) then

c....   Boomerang at "b" or "d".

c....   If ntype = 12, no points with fdk between aw / (aw - bw) and 1.0,
c....                             and fdl between 0.0 and -bw / (cw - bw).

c....   If ntype = 14, no points with fdk between 0.0 and -dw / (cw - dw),
c....                             and fdl between aw / (aw - dw) and 1.0.

c....   Sample from triangles "dab" and "bcd" separately.

        call apttloc (au, av, bu, bv, cu, cv, ns, 0, tol,
     &                dabu, dabv, wa, wb, wc, nerr)
        call apttloc (cu, cv, du, dv, au, av, ns, 0, tol,
     &                bcdu, bcdv, wa, wb, wc, nerr)

c....   Find a set of random numbers.

c---- Loop over local arrays.
        do 140 n = 1, ns
          ranfp(n) = ranf( )
c---- End of loop over local arrays.
  140   continue

c....   Select randomly from "dab" (weight aw) or "bcd" (weight cw).

c---- Loop over local arrays.
        do 150 n = 1, ns

          nn     = n + n1 - 1
          if (aw .ge. ranfp(n) * (aw + cw)) then
            pu(nn) = dabu(n)
            pv(nn) = dabv(n)
          else
            pu(nn) = bcdu(n)
            pv(nn) = bcdv(n)
          endif

c---- End of loop over local arrays.
  150   continue

      elseif ((ntype .eq. 11) .or.
     &        (ntype .eq. 13)) then

c....   Boomerang at "a" or "c".

c....   If ntype = 11, no points with fdk between 0.0 and -aw / (bw - aw),
c....                             and fdl between 0.0 and -aw / (dw - aw).

c....   If ntype = 13, no points with fdk between dw / (dw - cw) and 1.0,
c....                             and fdl between bw / (bw - cw) and 1.0.

c....   Sample from triangles "abc" and "cda" separately.

        call apttloc (au, av, bu, bv, cu, cv, ns, 0, tol,
     &                abcu, abcv, wa, wb, wc, nerr)
        call apttloc (cu, cv, du, dv, au, av, ns, 0, tol,
     &                cdau, cdav, wa, wb, wc, nerr)

c....   Find a set of random numbers.

c---- Loop over local arrays.
        do 160 n = 1, ns
          ranfp(n) = ranf( )
c---- End of loop over local arrays.
  160   continue

c....   Select randomly from "abc" (weight bw) or "cda" (weight dw).

c---- Loop over local arrays.
        do 170 n = 1, ns

          nn     = n + n1 - 1
          if (bw .ge. ranfp(n) * (bw + dw)) then
            pu(nn) = abcu(n)
            pv(nn) = abcv(n)
          else
            pu(nn) = cdau(n)
            pv(nn) = cdav(n)
          endif

c---- End of loop over local arrays.
  170   continue

c---- Bowtie with good edge "ab".
      elseif (ntype .eq. 21) then

c....   No points with fdl greater than min (aw / (bw - cw), bw / (bw - cw)).

c....   Sample from triangle "abq".

        call apttloc (au, av, bu, bv, qu, qv, ns, 0, tol,
     &                pu(n1), pv(n1), wa, wb, wc, nerr)

c---- Bowtie with good edge "bc".
      elseif (ntype .eq. 22) then

c....   No points with fdk less than max (-aw / (bw - aw), -dw / (bw - aw)).

c....   Sample from triangle "bcq".

        call apttloc (bu, bv, cu, cv, qu, qv, ns, 0, tol,
     &                pu(n1), pv(n1), wa, wb, wc, nerr)

c---- Bowtie with good edge "cd".
      elseif (ntype .eq. 23) then

c....   No points with fdl less than max (-aw / (dw - aw), -bw / (dw - aw)).

c....   Sample from triangle "cdq".

        call apttloc (cu, cv, du, dv, qu, qv, ns, 0, tol,
     &                pu(n1), pv(n1), wa, wb, wc, nerr)

c---- Bowtie with good edge "da".
      elseif (ntype .eq. 24) then

c....   No points with fdk greater than min (aw / (aw - bw), dw / (aw - bw)).

c....   Sample from triangle "daq".

        call apttloc (du, dv, au, av, qu, qv, ns, 0, tol,
     &                pu(n1), pv(n1), wa, wb, wc, nerr)

c---- Tested ntype.
      endif

c.... Find fdk, fdl for boomerang or bowtie.

c---- Quadrangle a boomerang or bowtie.
      if (ntype .ne. 0) then

        call aptqfdc (au, av, bu, bv, cu, cv, du, dv, pu(n1), pv(n1),
     &                ns, tol, fdk(n1), fdl(n1), ngood, nerr)

c---- Tested ntype.
      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 9903 format (/ 'aptqloc results:  fdk, fdl, pu, pv')
cbug      write ( 3, 9903)
cbug      call aptmean (fdk, np, 1.e-11, argmean, argdev, nerr)
cbug      call aptmean (fdl, np, 1.e-11, argmean, argdev, nerr)
cbug      call aptmean (pu,  np, 1.e-11, argmean, argdev, nerr)
cbug      call aptmean (pv,  np, 1.e-11, argmean, argdev, nerr)
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832