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