subroutine apttloc (au, av, bu, bv, cu, cv, np, kloc, tol,
& pu, pv, wa, wb, wc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTTLOC
c
c call apttloc (au, av, bu, bv, cu, cv, np, kloc, tol,
c & pu, pv, wa, wb, wc, nerr)
c
c Version: apttloc Updated 1990 December 3 16:20.
c apttloc 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 triangle in a major plane with
c vertices a = (au, av), b = (bu, bv), c = (cu, cv),
c in any order around the triangle, and optionally (kloc = 1), to
c return the local coordinates wa, wb and wc for each point "p":
c
c vn = (b - a) X (c - b) = (c - b) X (a - c) = (a - c) X (b - a)
c
c wa = (vn . (c - b) X (p - c)) / vn**2
c wb = (vn . (a - c) X (p - a)) / vn**2
c wb = (vn . (b - a) X (p - b)) / vn**2
c
c p = wa * a + wb * b + wc * c
c
c where "X" indicates the vector product, and "." indicates the
c dot product.
c
c The relationship of the triangle local coordinates wa, wb and wc
c to the quadrangle local coordinates fdk and fdl is as follows:
c
c wa = 1 - fdk fdk = wb + wc
c wb = fdk * (1 - fdl) fdl = wc / (wb + wc)
c wc = fdk * fdl
c
c where fdk is the fractional distance from edge "aa" to edge
c "bc", and fdl is the fractional distance from edge "ab" to
c edge "ca", and any point "p" is given by:
c
c p = a + fdk * (b - a) + fdk * fdl * (c - b)
c
c Flag nerr indicates any input error.
c
c History: 1990 May 10 17:00. Replaced arguments fdk, fdl with wa, wb, wc,
c and made calculation and return of wa, wb and wc optional,
c based on new input argument kloc.
c
c Input: au, av, bu, bv, cu, cv, np, kloc, tol.
c
c Output: pu, pv, wa, wb, wc, nerr.
c
c Calls: aptslid
c
c Glossary:
c
c au, av Input The u and v coordinates of vertex "a" of triangle.
c
c bu, bv Input The u and v coordinates of vertex "b" of triangle.
c
c cu, cv Input The u and v coordinates of vertex "c" of triangle.
c
c kloc Input Indicates local coordinates wa, wb and wc are to be
c returned, if 1. Otherwise, 0.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kloc is not 0 or 1.
c
c np Input Size of arrays fdk, fdl, pu, pv.
c Number of points "p" to sample.
c
c pu, pv Output Sampled point "p". Size np. Coordinates may be
c truncated to zero, if less than the estimated error
c in their calculation, based on tol.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c wa Output Fractional distance of point "p" to vertex "a" from
c edge "bc". Size np.
c
c wb Output Fractional distance of point "p" to vertex "b" from
c edge "ca". Size np.
c
c wc Output Fractional distance of point "p" to vertex "c" from
c edge "ab". Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cooodinate u of a sampled point "p".
dimension pu (1)
c---- Coordinate v of a sampled point "p".
dimension pv (1)
c---- Fractional distance to "a" from "bc".
dimension wa (1)
c---- Fractional distance to "b" from "ca".
dimension wb (1)
c---- Fractional distance to "c" from "ab".
dimension wc (1)
c.... Local variables.
c---- Fractional distance from "a" to "bc".
common /lapttloc/ fdk (64)
c---- Fractional distance from "ab" to "ca".
common /lapttloc/ fdl (64)
c---- Index in local array.
common /lapttloc/ n
c---- First index of subset of data.
common /lapttloc/ n1
c---- Last index of subset of data.
common /lapttloc/ n2
c---- Index in external array.
common /lapttloc/ nn
c---- Size of current subset of data.
common /lapttloc/ ns
c---- Estimated error in pu.
common /lapttloc/ puerr
c---- Estimated error in pv.
common /lapttloc/ pverr
cbugc***DEBUG begins.
cbug common /lapttloc/ avgu, avgv, devu, devv
cbug common /lapttloc/ avgu2, avgv2
cbug common /lapttloc/ avgwa, avgwb, avgwc
cbug common /lapttloc/ devwa, devwb, devwc
cbug common /lapttloc/ avgwa2, avgwb2, avgwc2
cbug 9901 format (/ 'apttloc sampling from a triangle. np=',i3 /
cbug & ' au,av=',1p2e22.14 /
cbug & ' bu,bv=',1p2e22.14 /
cbug & ' cu,cv=',1p2e22.14)
cbug write (3, 9901) np, au, av, bu, bv, cu, cv
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((kloc .lt. 0) .or. (kloc .gt. 1)) then
nerr = 2
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
c.... Loop over the data subsets.
110 ns = n2 - n1 + 1
c.... Sample ns values of fdk from a linear distribution, with
c.... probability 0.0 at vertex "a", probability 2.0 at edge "bc".
call aptslid (0.0, 0.0, 1.0, 2.0, ns, fdk, nerr)
c.... Sample ns values of fdl from a uniform distribution, with
c.... equal probabilities at sides "ab" and "ca".
c---- Loop over subset of data.
do 120 n = 1, ns
fdl(n) = ranf( )
c---- End of loop over subset of data.
120 continue
c.... Find the coordinates of the points.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
pu(nn) = au + fdk(n) * (bu - au) + fdk(n) * fdl(n) * (cu - bu)
pv(nn) = av + fdk(n) * (bv - av) + fdk(n) * fdl(n) * (cv - bv)
c---- End of loop over subset of data.
130 continue
c.... See if truncation allowed.
c---- Truncation to zero allowed.
if (tol .gt. 0.0) then
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
puerr = tol * (fdk(n) * abs (bu - au) +
& fdk(n) * fdl(n) * abs (cu - bu))
pverr = tol * (fdk(n) * abs (bv - av) +
& fdk(n) * fdl(n) * abs (cv - bv))
if (abs (pu(nn)) .lt. puerr) then
pu(nn) = 0.0
endif
if (abs (pv(nn)) .lt. pverr) then
pv(nn) = 0.0
else
endif
c---- End of loop over subset of data.
140 continue
c---- Tested tol.
endif
c.... See if the local triangular coordinates are to be returned.
c---- Find local triangular coordinates.
if (kloc .eq. 1) then
c---- Loop over subset of data.
do 150 n = 1, ns
nn = n + n1 - 1
wa(nn) = 1.0 - fdk(n)
wb(nn) = fdk(n) * (1.0 - fdl(n))
wc(nn) = fdk(n) * fdl(n)
c---- End of loop over subset of data.
150 continue
c---- Tested kloc.
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
cbug avgu = 0.0
cbug avgv = 0.0
cbug avgu2 = 0.0
cbug avgv2 = 0.0
cbug
cbug do 160 n = 1, np
cbug avgu = avgu + pu(n)
cbug avgv = avgv + pv(n)
cbug avgu2 = avgu2 + pu(n)**2
cbug avgv2 = avgv2 + pv(n)**2
cbug 160 continue
cbug
cbug avgu = avgu / np
cbug avgv = avgv / np
cbug avgu2 = avgu2 / np
cbug avgv2 = avgv2 / np
cbug devu = sqrt (avgu2 - avgu**2)
cbug devv = sqrt (avgv2 - avgv**2)
cbug
cbug 9902 format (/ 'apttloc results:')
cbug 9903 format (i3,' pu,pv=',1p2e22.14)
cbug 9904 format (/ 'apttloc mean and deviation:' /
cbug & ' avg(u,v)=',1p2e22.14 /
cbug & ' dev(u,v)=',1p2e22.14)
cbug
cbug write ( 3, 9902)
cbug write ( 3, 9903) (n, pu(n), pv(n), n = 1, np)
cbug write ( 3, 9904) avgu, avgv, devu, devv
cbug
cbugc---- Local coordinates returned.
cbug if (kloc .eq. 1) then
cbug
cbug avgwa = 0.0
cbug avgwb = 0.0
cbug avgwc = 0.0
cbug avgwa2 = 0.0
cbug avgwb2 = 0.0
cbug avgwc2 = 0.0
cbug
cbug do 170 n = 1, np
cbug avgwa = avgwa + wa(n)
cbug avgwb = avgwb + wb(n)
cbug avgwc = avgwc + wc(n)
cbug avgwa2 = avgwa2 + wa(n)**2
cbug avgwb2 = avgwb2 + wb(n)**2
cbug avgwc2 = avgwc2 + wc(n)**2
cbug 170 continue
cbug
cbug avgwa = avgwa / np
cbug avgwb = avgwb / np
cbug avgwc = avgwc / np
cbug avgwa2 = avgwa2 / np
cbug avgwb2 = avgwb2 / np
cbug avgwc2 = avgwc2 / np
cbug devwa = sqrt (avgwa2 - avgwa**2)
cbug devwb = sqrt (avgwb2 - avgwb**2)
cbug devwc = sqrt (avgwc2 - avgwc**2)
cbug
cbug 9905 format (i3,' wa,wb,wc=',1p3e22.14)
cbug 9906 format (/ 'apttloc mean and deviation:' /
cbug & ' avgw(a,b,c)=',1p3e22.14 /
cbug & ' devw(a,b,c)=',1p3e22.14)
cbug
cbug write ( 3, 9905) (n, wa(n), wb(n), wc(n), n = 1, np)
cbug write ( 3, 9906) avgwa, avgwb, avgwc,
cbug & devwa, devwb, devwc
cbug
cbugc---- Tested kloc.
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine apttloc. (+1 line.)
end
UCRL-WEB-209832