subroutine aptqfdc (au, av, bu, bv, cu, cv, du, dv, pu, pv, np, & tol, fdk, fdl, ngood, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTQFDC c c call aptqfdc (au, av, bu, bv, cu, cv, du, dv, pu, pv, np, tol, c & fdk, fdl, ngood, nerr) c c Version: aptqfdc Updated 1993 December 6 15:20. c aptqfdc Originated 1990 February 14 16:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for each of the np points p = (pu, pv), the fractional c distances fdk and fdl of point "p" between the opposite edges c "da" and "bc", and "ab" and "cd", respectively, of the 2-D c quadrangle in the uv plane with vertices a = (au, av), c 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 c bu, bv Input The u and v coordinates of point "b" in the uv plane. c c cu, cv Input The u and v coordinates of point "c" in the uv plane. c c du, dv Input The u and v coordinates of point "d" in the uv plane. 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 nerr Output Indicates an input error, if not 0. c 1 if np is not positive. 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 np Input Size of arrays pu, pv, 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---- 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 /laptqfdc/ aa (64) c---- Component u of vector "ab". common /laptqfdc/ abu c---- Component v of vector "ab". common /laptqfdc/ abv c---- Coefficient of fdk in quadratic. common /laptqfdc/ bb (64) c---- Constant term in quadratic. common /laptqfdc/ cc (64) c---- Component u of vector "cd". common /laptqfdc/ cdu c---- Component v of vector "cd". common /laptqfdc/ cdv c---- Component u of vector "da". common /laptqfdc/ dau c---- Component v of vector "da". common /laptqfdc/ dav c---- Intermediate value. common /laptqfdc/ denu1 c---- Intermediate value. common /laptqfdc/ denu2 c---- Intermediate value. common /laptqfdc/ denv1 c---- Intermediate value. common /laptqfdc/ denv2 c---- First real root (-1.e+99 if none). common /laptqfdc/ fdk1 (64) c---- Second real root (-1.e+99 if none). common /laptqfdc/ fdk2 (64) c---- Value of fdl from fdk1. common /laptqfdc/ fdl1 (64) c---- Value of fdl from fdk2. common /laptqfdc/ fdl2 (64) c---- A very small number. common /laptqfdc/ fuz c---- Truncation error indicator. common /laptqfdc/ itrun (64) c---- Index in arrays. common /laptqfdc/ n c---- First index of subset of data. common /laptqfdc/ n1 c---- Last index of subset of data. common /laptqfdc/ n2 c---- 1 if fdk1, fdl1 acceptable. common /laptqfdc/ ngood1 (64) c---- 1 if fdk2, fdl2 acceptable. common /laptqfdc/ ngood2 (64) c---- 1 if fdk or fdl range adjusted. common /laptqfdc/ nlim (64) c---- Index in external array. common /laptqfdc/ nn c---- Size of current subset of data. common /laptqfdc/ ns c---- Number of real roots of quadratic. common /laptqfdc/ nroots (64) c---- Component u of vector "pa". common /laptqfdc/ pau (64) c---- Component v of vector "pa". common /laptqfdc/ pav (64) c---- Component u of vector "pd". common /laptqfdc/ pdu (64) c---- Component v of vector "pd". common /laptqfdc/ pdv (64) c---- Factor bb**2 - 4.0 * aa * cc. common /laptqfdc/ qq (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptqfdc finding fdk and fdl in quadrangle:' / cbug & ' au,av= ',1p2e22.14,' (vertices)' / cbug & ' bu,bv= ',1p2e22.14 / cbug & ' cu,cv= ',1p2e22.14 / cbug & ' du,dv= ',1p2e22.14 / cbug & (i3,' pu,pv= ',1p2e22.14,' (point)')) cbug write ( 3, 9901) au, av, bu, bv, cu, cv, du, dv, cbug & (n, pu(n), pv(n), 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 quadrangle edge vectors. abu = bu - au abv = bv - av cdu = du - cu cdv = dv - cv dau = au - du dav = av - dv c.... Find the point position vectors. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 pau(n) = au - pu(nn) pav(n) = av - pv(nn) pdu(n) = du - pu(nn) pdv(n) = dv - pv(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 * abu - abv * cdu bb(n) = abv * pdu(n) - abu * pdv(n) + & cdv * pau(n) - cdu * pav(n) cc(n) = dav * pau(n) - dau * 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 + fdk1(n) * (abu + cdu) denv1 = dav + fdk1(n) * (abv + cdv) fdl1(n) = ((pau(n) + fdk1(n) * abu) * denu1 + & (pav(n) + fdk1(n) * abv) * 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 nn = n + n1 - 1 denu2 = dau + fdk2(n) * (abu + cdu) denv2 = dav + fdk2(n) * (abv + cdv) fdl2(n) = ((pau(n) + fdk2(n) * abu) * denu2 + & (pav(n) + fdk2(n) * abv) * 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 c---- End of loop over subset of data. 150 continue c.... Select any good solution. c---- Loop over subset of data. do 160 n = 1, ns nn = n + n1 - 1 if ((ngood2(n) .eq. 1) .and. & (ngood1(n) .ne. 1) ) then fdk(nn) = fdk2(n) fdl(nn) = fdl2(n) else fdk(nn) = fdk1(n) fdl(nn) = fdl1(n) endif ngood(nn) = ngood1(n) + ngood2(n) c---- End of loop over subset of data. 160 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 (/ 'aptqfdc 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 (/ 'aptqfdc 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 aptqfdc. (+1 line.) end UCRL-WEB-209832