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

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTFDQC
c
c     call aptfdqc (au, av, bu, bv, cu, cv, du, dv, pu, pv, np, tol,
c    &              fdk, fdl, ngood, nerr)
c
c     Version:  aptfdqc  Updated    1993 December 6 15:20..
c               aptfdqc  Originated 1990 January 26 16:20.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find, for each of the np sets of input data, the fractional
c               distances fdk and fdl of the point p = (pu, pv), between the
c               opposite edges "da" and "bc", and "ab" and "cd", respectively,
c               of the 2-D quadrangle in the uv plane with vertices
c               a = (au, av), b = (bu, bv), c = (cu, cv), and d = (du, dv).
c
c                 p = a + fdk * (b - a) + fdl * (d - a)
c                     + fdk * fdl * (a - b + c - d)
c
c                 p = (1 - fdk) * (1 - fdl) * a + fdk * (1 - fdl) * b
c                     + fdk * fdl * c + (1 - fdk) * fdl * d
c
c               Flag ngood indicates if the point is in the quadrangle.
c               Flag nerr indicates any input error.
c
c     Input:    au, av, bu, bv, cu, cv, du, dv, pu, pv, np, tol.
c
c     Output:   fdk, fdl, ngood, nerr.
c
c     Calls: aptqrtv 
c
c     Glossary:
c
c     au, av    Input    The u and v coordinates of point "a" in the uv plane.
c                          Size np.
c
c     bu, bv    Input    The u and v coordinates of point "b" in the uv plane.
c                          Size np.
c
c     cu, cv    Input    The u and v coordinates of point "c" in the uv plane.
c                          Size np.
c
c     du, dv    Input    The u and v coordinates of point "d" in the uv plane.
c                          Size np.
c
c     fdk       Output   Fractional distance of point "p" from the line segment
c                          "da" to the line segment "bc".  Size np.
c                          Values between -tol and tol will be adjusted to tol.
c                          Values between 1.0 - tol and 1.0 + tol will be
c                          adjusted to 1.0 - tol.  See ngood.
c
c     fdl       Output   Fractional distance of point "p" from the line segment
c                          "ab" to the line segment "cd".  Size np.
c                          Values between -tol and tol will be adjusted to tol.
c                          Values between 1.0 - tol and 1.0 + tol will be
c                          adjusted to 1.0 - tol.  See ngood.
c
c     ngood     Output   Indicates range of fractional distances fdk and fdl:
c                          0 if either fdk or fdl is outside the range from
c                            -tol to 1.0 + tol.  Also if both are.
c                          1 if fdk and fdl are both between -tol and 1.0 + tol.
c                            This can be true even when point "p" is outside
c                            a boomeranged or bowtied quadrangle.
c                          2 if fdk and fdl are both between -tol and 1.0 + tol,
c                            and two possible solutions exist.  This can happen
c                            when the quadrangle is a boomerang or bowtie.
c                            Only one of the solutions is returned.
c                          Size np.
c
c     nerr      Output   Indicates an input error, if not 0.
c                          1 if np is not positive.
c
c     np        Input    Size of arrays pu, pv, au, av, bu, bv, cu, cv, du, dv,
c                          fdk, fdl, ngood.
c
c     pu, pv    Input    The u and v coordinates of point "p" in the uv plane.
c                          Size np.
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 u of vertex "a" of quad.
      dimension au      (1)
c---- Coordinate v of vertex "a" of quad.
      dimension av      (1)
c---- Coordinate u of vertex "b" of quad.
      dimension bu      (1)
c---- Coordinate v of vertex "b" of quad.
      dimension bv      (1)
c---- Coordinate u of vertex "c" of quad.
      dimension cu      (1)
c---- Coordinate v of vertex "c" of quad.
      dimension cv      (1)
c---- Coordinate u of vertex "d" of quad.
      dimension du      (1)
c---- Coordinate v of vertex "d" of quad.
      dimension dv      (1)
c---- Frac dist of "p" from "da" to "bc".
      dimension fdk     (1)
c---- Frac dist of "p" from "ab" to "cd".
      dimension fdl     (1)
c---- 1 if good fdk, fdl found.
      dimension ngood   (1)
c---- Coordinate u of point "p".
      dimension pu      (1)
c---- Coordinate v of point "p".
      dimension pv      (1)

c.... Local variables.

c---- Coefficient of fdk**2 in quadratic.
      common /laptfdqc/ aa      (64)
c---- Component u of vector "ab".
      common /laptfdqc/ abu     (64)
c---- Component v of vector "ab".
      common /laptfdqc/ abv     (64)
c---- Coefficient of fdk in quadratic.
      common /laptfdqc/ bb      (64)
c---- Constant term in quadratic.
      common /laptfdqc/ cc      (64)
c---- Component u of vector "cd".
      common /laptfdqc/ cdu     (64)
c---- Component v of vector "cd".
      common /laptfdqc/ cdv     (64)
c---- Component u of vector "da".
      common /laptfdqc/ dau     (64)
c---- Component v of vector "da".
      common /laptfdqc/ dav     (64)
c---- Intermediate value.
      common /laptfdqc/ denu1
c---- Intermediate value.
      common /laptfdqc/ denu2
c---- Intermediate value.
      common /laptfdqc/ denv1
c---- Intermediate value.
      common /laptfdqc/ denv2
c---- First real root (-1.e+99 if none).
      common /laptfdqc/ fdk1    (64)
c---- Second real root (-1.e+99 if none).
      common /laptfdqc/ fdk2    (64)
c---- Value of fdl from fdk1.
      common /laptfdqc/ fdl1    (64)
c---- Value of fdl from fdk2.
      common /laptfdqc/ fdl2    (64)

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

c---- Truncation error indicator.
      common /laptfdqc/ itrun   (64)
c---- Index in arrays.
      common /laptfdqc/ n
c---- First index of subset of data.
      common /laptfdqc/ n1
c---- Last index of subset of data.
      common /laptfdqc/ n2
c---- 1 if fdk1, fdl1 acceptable.
      common /laptfdqc/ ngood1  (64)
c---- 1 if fdk2, fdl2 acceptable.
      common /laptfdqc/ ngood2  (64)
c---- 1 if fdk or fdl range adjusted.
      common /laptfdqc/ nlim    (64)
c---- Index in external array.
      common /laptfdqc/ nn
c---- Size of current subset of data.
      common /laptfdqc/ ns
c---- Number of real roots of quadratic.
      common /laptfdqc/ nroots  (64)
c---- Component u of vector "pa".
      common /laptfdqc/ pau     (64)
c---- Component v of vector "pa".
      common /laptfdqc/ pav     (64)
c---- Component u of vector "pd".
      common /laptfdqc/ pdu     (64)
c---- Component v of vector "pd".
      common /laptfdqc/ pdv     (64)
c---- Factor bb**2 - 4.0 * aa * cc.
      common /laptfdqc/ qq      (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptfdqc finding fdk and fdl in quadrangle:' /
cbug     &  (i3,' pu,pv=  ',1p2e22.14,' (point)' /
cbug     &  '    au,av=  ',1p2e22.14,' (vertices)' /
cbug     &  '    bu,bv=  ',1p2e22.14 /
cbug     &  '    cu,cv=  ',1p2e22.14 /
cbug     &  '    du,dv=  ',1p2e22.14))
cbug      write ( 3, 9901) (n, pu(n), pv, au(n), av(n),
cbug     &  bu(n), bv(n), cu(n), cv(n), du(n), dv(n),
cbug     &  n = 1, np)
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.... Set up the indices of the first subset of data.

      n1 = 1
      n2 = min (np, 64)

  110 ns = n2 - n1 + 1

c.... Find the vectors of the point position, quadrangle edges.

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

        nn = n + n1 - 1

        pau(n) = au(nn) - pu(nn)
        pav(n) = av(nn) - pv(nn)
        pdu(n) = du(nn) - pu(nn)
        pdv(n) = dv(nn) - pv(nn)

        abu(n) = bu(nn) - au(nn)
        abv(n) = bv(nn) - av(nn)
        cdu(n) = du(nn) - cu(nn)
        cdv(n) = dv(nn) - cv(nn)
        dau(n) = au(nn) - du(nn)
        dav(n) = av(nn) - dv(nn)

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

c.... Find the coefficients of the quadratic equation for fdk.
c....   aa * fdk**2 + bb * fdk + cc = 0.

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

        aa(n) = cdv(n) * abu(n) - abv(n) * cdu(n)
        bb(n) = abv(n) * pdu(n) - abu(n) * pdv(n) +
     &          cdv(n) * pau(n) - cdu(n) * pav(n)
        cc(n) = dav(n) * pau(n) - dau(n) * pav(n)

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

c.... Find any possible values of fdk.

      call aptqrtv (0, aa, bb, cc, qq, ns, tol,
     &              nroots, fdk1, fdk2, itrun, nerr)

c.... Find any acceptable value of fdl corresponding to fdk1.

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

        nn = n + n1 - 1

        denu1   = dau(n) + fdk1(n) * (abu(n) + cdu(n))
        denv1   = dav(n) + fdk1(n) * (abv(n) + cdv(n))

        fdl1(n) = ((pau(n) + fdk1(n) * abu(n)) * denu1 +
     &             (pav(n) + fdk1(n) * abv(n)) * denv1) /
     &            (denu1**2 + denv1**2 + fuz)

        if ((fdk1(n) .ge. -tol)        .and.
     &      (fdk1(n) .le. (1.0 + tol)) .and.
     &      (fdl1(n) .ge. -tol)        .and.
     &      (fdl1(n) .le. (1.0 + tol))       ) then
          ngood1(n) = 1
        else
          ngood1(n) = 0
        endif


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

c.... Find any acceptable value of fdl corresponding to fdk2.

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

        denu2   = dau(n) + fdk2(n) * (abu(n) + cdu(n))
        denv2   = dav(n) + fdk2(n) * (abv(n) + cdv(n))

        fdl2(n) = ((pau(n) + fdk2(n) * abu(n)) * denu2 +
     &             (pav(n) + fdk2(n) * abv(n)) * denv2) /
     &            (denu2**2 + denv2**2 + fuz)

        if ((fdk2(n) .ge. -tol)        .and.
     &      (fdk2(n) .le. (1.0 + tol)) .and.
     &      (fdl2(n) .ge. -tol)        .and.
     &      (fdl2(n) .le. (1.0 + tol))       ) then
          ngood2(n) = 1
        else
          ngood2(n) = 0
        endif

        nn = n + n1 - 1

        if ((ngood2(n) .eq. 1) .and.
     &      (ngood1(n) .ne. 1)       ) then
          fdk(nn) = fdk2(n)
        else
          fdk(nn) = fdk1(n)
        endif

        if ((ngood2(n) .eq. 1) .and.
     &      (ngood1(n) .ne. 1)       ) then
          fdl(nn) = fdl2(n)
        else
          fdl(nn) = fdl1(n)
        endif

        ngood(nn) = ngood1(n) + ngood2(n)

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

c.... See if the ranges of fdk, fdl should be adjusted.

c---- Adjust the ranges of fdk, fdl.
      if (tol .gt. 0.0) then

        call aptfdav (fdk(n1), ns, 1, tol, nlim, nerr)
        call aptfdav (fdl(n1), ns, 1, tol, nlim, nerr)

c---- Tested tol.
      endif
cbugc***DEBUG begins.
cbug 9801 format (/ 'aptfdqc intermediate results:' /
cbug     &  (i3,' fdk1,fdl1=',1p2e22.14 /
cbug     &  '    fdk2,fdl2=',1p2e22.14))
cbug      write ( 3, 9801) (n, fdk1(n), fdl1(n), fdk2(n), fdl2(n),
cbug     &  n = 1, np)
cbugc***DEBUG ends.

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 9902 format (/ 'aptfdqc results:' /
cbug     &  (i3,' fdk,fdl=',1p2e22.14,' ngood=',i2))
cbug      write ( 3, 9902) (n, fdk(n), fdl(n), ngood(n), n = 1, np)
cbugc***DEBUG ends.

  210 return

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

UCRL-WEB-209832