subroutine aptinvp (ax, ay, az, px, py, pz, np, tol, rinv)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTINVP
c
c call aptinvp (ax, ay, az, px, py, pz, np, tol, rinv)
c
c Version: aptinvp Updated 1990 November 28 10:00.
c aptinvp Originated 1989 November 2 14:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the matrix operator rinv for inversion through the
c origin, and to do an equivalent inversion through the point
c a = (ax, ay, az), of the np points or vectors
c p = (px, py, pz). If "p" are unbound vectors, point "a" must
c be at the origin. The new components of "p" will be truncated
c to zero if less than the estimated error in their calculation,
c based on tol.
c
c Input: itype, ax, ay, az, px, py, pz, np, tol.
c
c Output: px, py, pz, rinv.
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of the inversion point "a".
c
c np Input Number of points or vectors (px, py, pz). May be 0.
c
c px,py,pz In/Out The x, y, z coordinates of a point, or
c the x, y, z components of a vector. Size np.
c Truncated to zero if less than the estimated error
c in their calculation. See tol.
c
c rinv Output Array rinv(3,3). Inversion operator. Diagonal
c elements are -1. Off-diagonal elements are 0.
c
c tol Input Numerical tolerance limit. Used to test and adjust
c point components.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate or component px. Size np.
dimension px (1)
c---- Coordinate or component py. Size np.
dimension py (1)
c---- Coordinate or component pz. Size np.
dimension pz (1)
c---- Inversion matrix operator.
dimension rinv (3,3)
c.... Local variables.
c---- Index of point or vector.
common /laptinvp/ n
c---- Estimated error in px.
common /laptinvp/ pxerr
c---- Estimated error in py.
common /laptinvp/ pyerr
c---- Estimated error in pz.
common /laptinvp/ pzerr
cbugc***DEBUG begins.
cbugc---- Row index of inversion matrix.
cbug common /laptinvp/ i
cbugc---- Column index of inversion matrix.
cbug common /laptinvp/ j
cbug 9901 format (/ 'aptinvp. Inverting through point' /
cbug & ' ax,ay,az=',1p3e22.14)
cbug write ( 3, 9901) ax, ay, az
cbugc***DEBUG ends.
c.... Form the components of the inversion matrix.
c.... Row and column vectors are orthogonal unit vectors.
rinv(1,1) = -1.0
rinv(1,2) = 0.0
rinv(1,3) = 0.0
rinv(2,1) = 0.0
rinv(2,2) = -1.0
rinv(2,3) = 0.0
rinv(3,1) = 0.0
rinv(3,2) = 0.0
rinv(3,3) = -1.0
cbugc***DEBUG begins.
cbug 9902 format (/ ' rinv=',3(/ 1p3e22.14))
cbug write ( 3, 9902) ((rinv(i,j), j = 1, 3), i = 1, 3)
cbugc***DEBUG ends.
c.... Do the inversion operation on the np points (px, py, pz).
c.... Truncate very small coordinates or components to zero.
c---- Points or vectors exist.
if (np .gt. 0) then
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptinvp inverting points:')
cbug 9904 format (i3,' px,py,pz=',1p3e22.14)
cbug write ( 3, 9903)
cbug write ( 3, 9904) (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
c---- No truncation of (px, py, pz).
if (tol .le. 0.0) then
c---- Loop over points or vectors.
do 110 n = 1, np
px(n) = 2.0 * ax - px(n)
py(n) = 2.0 * ay - py(n)
pz(n) = 2.0 * az - pz(n)
c---- End of loop over points or vectors.
110 continue
c---- If tol is positive.
else
c---- Loop over points or vectors.
do 120 n = 1, np
pxerr = tol * (2.0 * abs (ax) + abs (px(n)))
pyerr = tol * (2.0 * abs (ay) + abs (py(n)))
pzerr = tol * (2.0 * abs (az) + abs (pz(n)))
px(n) = 2.0 * ax - px(n)
py(n) = 2.0 * ay - py(n)
pz(n) = 2.0 * az - pz(n)
if (abs (px(n)) .lt. pxerr) then
px(n) = 0.0
endif
if (abs (py(n)) .lt. pyerr) then
py(n) = 0.0
endif
if (abs (pz(n)) .lt. pzerr) then
pz(n) = 0.0
endif
c---- End of loop over points or vectors.
120 continue
c---- Tested tol.
endif
cbugc***DEBUG begins.
cbug 9905 format (/ 'aptinvp inverted points:')
cbug write ( 3, 9905)
cbug write ( 3, 9904) (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
c---- Tested np.
endif
return
c.... End of subroutine aptinvp. (+1 line.)
end
UCRL-WEB-209832