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