subroutine aptquwt (fdk, fdl, ax, ay, az, bx, by, bz,
& cx, cy, cz, dx, dy, dz, np, tol,
& px, py, pz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQUWT
c
c call aptquwt (fdk, fdl, ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, px, py, pz, nerr)
c
c Version: aptquwt Updated 1990 November 29 15:30.
c aptquwt Originated 1990 May 17 11:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np sets of local quadrangle coordinates
c fdk and fdl, the point p = (px, py, pz), interpolated between
c vertices a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz)
c and d = (dx, dy, dz) of a quadrangle:
c
c p = a + fdk*(b - a) + fdl*(d - a) + fdk*fdl*(a - b + c - d)
c
c p = (1-fdk)*(1-fdl)*a + fdk*(1-fdl)*b + fdk*fdl*c +
c (1-fdk)*fdl*d
c
c p = wa * a + wb * b + wc * c + wd * d
c
c fdk = wb + wc = 1 - wa - wd
c fdl = wc + wd = 1 - wa - wb
c
c wa + wb + wc + wd = 1
c
c Flag nerr indicates any input error.
c
c Input: fdk, fdl, ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz,
c np, tol.
c
c Output: px, py, pz, nerr.
c
c Calls: None
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a" of quadrangle.
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of quadrangle.
c Size np.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of quadrangle.
c Size np.
c
c dx,dy,dz Input The x, y, z coordinates of vertex "d" of quadrangle.
c Size np.
c
c fdk Input Fractional distance of point "p" between sides "da"
c and edge "bc" of quadrangle. Size np.
c
c fdl Input Fractional distance of point "p" between sides "ab"
c and "cd" of quadrangle. 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 fdk, fdl, ax, ay, az, bx, by, bz,
c cx, cy, cz, dx, dy, dz, px, py, pz.
c
c px,py,pz Output Interpolated 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. With 64-bit floating
c point numbers, 1.e-5 to 1.e-11 is recommended.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cooodinate x of point "a".
dimension ax (1)
c---- Coordinate y of point "a".
dimension ay (1)
c---- Coordinate z of point "a".
dimension az (1)
c---- Cooodinate x of point "b".
dimension bx (1)
c---- Coordinate y of point "b".
dimension by (1)
c---- Coordinate z of point "b".
dimension bz (1)
c---- Cooodinate x of point "c".
dimension cx (1)
c---- Coordinate y of point "c".
dimension cy (1)
c---- Coordinate z of point "c".
dimension cz (1)
c---- Cooodinate x of point "d".
dimension dx (1)
c---- Coordinate y of point "d".
dimension dy (1)
c---- Coordinate z of point "d".
dimension dz (1)
c---- Fract. dist. between "da" and "bc".
dimension fdk (1)
c---- Fract. dist. between "ab" and "cd".
dimension fdl (1)
c---- Cooodinate x of point "p".
dimension px (1)
c---- Coordinate y of point "p".
dimension py (1)
c---- Coordinate z of point "p".
dimension pz (1)
c.... Local variables.
c---- Index in local array.
common /laptquwt/ n
c---- First index of subset of data.
common /laptquwt/ n1
c---- Last index of subset of data.
common /laptquwt/ n2
c---- Index in external array.
common /laptquwt/ nn
c---- Size of current subset of data.
common /laptquwt/ ns
c---- Estimated error in px.
common /laptquwt/ pxerr
c---- Estimated error in py.
common /laptquwt/ pyerr
c---- Estimated error in pz.
common /laptquwt/ pzerr
c---- Weight of vertex "a".
common /laptquwt/ wa
c---- Weight of vertex "b".
common /laptquwt/ wb
c---- Weight of vertex "c".
common /laptquwt/ wc
c---- Weight of vertex "d".
common /laptquwt/ wd
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptquwt interpolating in a quadrangle.' /
cbug & (i3,' fdk,fdl= ',1p2e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' dx,dy,dz=',1p3e22.14))
cbug write (3, 9901) (n, fdk(n), fdl(n),
cbug & ax(n), ay(n), az(n), bx(n), by(n), bz(n),
cbug & cx(n), cy(n), cz(n), dx(n), dy(n), dz(n), n = 1, np)
cbugc***DEBUG ends.
c.... initialize.
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)
c.... Loop over the data subsets.
110 ns = n2 - n1 + 1
c.... Find the coordinates of the points.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
wa = (1.0 - fdk(nn)) * (1.0 - fdl(nn))
wb = fdk(nn) * (1.0 - fdl(nn))
wc = fdk(nn) * fdl(nn)
wd = (1.0 - fdk(nn)) * fdl(nn)
px(nn) = wa * ax(nn) + wb * bx(nn) + wc * cx(nn) + wd * dx(nn)
py(nn) = wa * ay(nn) + wb * by(nn) + wc * cy(nn) + wd * dy(nn)
pz(nn) = wa * az(nn) + wb * bz(nn) + wc * cz(nn) + wd * dz(nn)
c---- End of loop over subset of data.
120 continue
c.... See if truncation allowed.
c---- Truncation to zero allowed.
if (tol .gt. 0.0) then
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
wa = (1.0 - fdk(nn)) * (1.0 - fdl(nn))
wb = fdk(nn) * (1.0 - fdl(nn))
wc = fdk(nn) * fdl(nn)
wd = (1.0 - fdk(nn)) * fdl(nn)
pxerr = tol * (abs (wa * ax(nn)) + abs (wb * bx(nn)) +
& abs (wc * cx(nn)) + abs (wd * dx(nn)))
pyerr = tol * (abs (wa * ay(nn)) + abs (wb * by(nn)) +
& abs (wc * cy(nn)) + abs (wd * dy(nn)))
pzerr = tol * (abs (wa * az(nn)) + abs (wb * bz(nn)) +
& abs (wc * cz(nn)) + abs (wd * dz(nn)))
if (abs (px(nn)) .lt. pxerr) then
px(nn) = 0.0
endif
if (abs (py(nn)) .lt. pyerr) then
py(nn) = 0.0
endif
if (abs (pz(nn)) .lt. pzerr) then
pz(nn) = 0.0
endif
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
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 9902 format (/ 'aptquwt results:' /
cbug & (i3,' px,py,pz=',1p3e22.14))
cbug write ( 3, 9902) (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptquwt. (+1 line.)
end
UCRL-WEB-209832