subroutine aptmopv (smat, inv, ax, ay, az, px, py, pz, np, tol, & nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTMOPV c c call aptmopv (smat, inv, ax, ay, az, px, py, pz, np, tol, nerr) c c Version: aptmopv Updated 1990 November 29 11:40. c aptmopv Originated 1989 November 2 14:10. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To transform the np points or vectors p = (px, py, pz) c by subtracting a = (ax, ay, az), then operating on the result c with the 3 x 3 matrix operator smat (inv = 0) or its transpose c (inv = 1), then adding "a". c Components of "p" smaller than the estimated error in their c calculation, based on tol, will be truncated to zero. c If np is not positive, nerr = 1 will be returned. c If inv is not 0 or 1, nerr = 2 will be returned. c c (1) if smat is a rotation operator (unitary, orthogonal, c the transpose = the inverse, and the the cosine of the c rotation angle is (trace - 1) / 2), the rotation is around c an axis through the point "a". c Option inv = 0 rotates the 3 row vectors of smat to be c parallel to the major axes, or rotates the major axes c to be parallel to the 3 column vectors of smat. c Option inv = 1 rotates the 3 column vectors of smat to c be parallel to the major axes, or rotates the major axes c to be parallel to the 3 row vectors of smat. c c (2) if smat is a reflection operator (unitary, symmetric, c its own inverse, and the trace = 1) the reflection is in a c plane through the point "a"; inv has no effect. If the values c "p" are unbound vectors, then "a" must be (0., 0., 0.). c The sequential application of two reflections is a rotation c and a possible translation. c c (3) if smat is an inversion operator (diagonals = -1, other c elements = 0), the inversion is through the point "a". c c History: 1990 March 13. Deleted truncation of components based on c total magnitude of vector. c c Input: smat, inv, ax, ay, az, px, py, pz, np, tol. c c Output: px, py, pz, nerr. c c Calls: apttran c c Glossary: c c ax,ay,az Input The x, y, z coordinates of invariant point "a". c Should be (0., 0., 0.) if "p" is a vector. c c inv Input 0 to operate with the matrix smat, c 1 to operate with its transpose. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c 2 if inv is not 0 or 1. c c np Input Size of arrays px, py, pz. Must be positive. c c px,py,pz In/Out The x, y, z coordinates of a point, or the c x, y, z components of a vector. Size np. c Components smaller than the estimated error in their c calculation, based on tol, will be truncated to zero. c c smat Input Array smat(3,3). A matrix operator. 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 or component x. Size np. dimension px (1) c---- Coordinate or component y. Size np. dimension py (1) c---- Coordinate or component z. Size np. dimension pz (1) c---- Matrix operator. dimension smat (3,3) c.... Local variables. c---- Index of point or vector. common /laptmopv/ n c---- Initial value of px. common /laptmopv/ ppx c---- Initial value of py. common /laptmopv/ ppy c---- Initial value of pz. common /laptmopv/ ppz c---- Estimated error in px. common /laptmopv/ pxerr c---- Estimated error in py. common /laptmopv/ pyerr c---- Estimated error in pz. common /laptmopv/ pzerr cbugc***DEBUG begins. cbugc---- Array index. cbug common /laptmopv/ i cbugc---- Array index. cbug common /laptmopv/ j cbug 9901 format (/ 'aptmopv operating with a matrix operator.' / cbug & ' np=',i3,' tol=',1pe13.5 / cbug & ' smat=',3(/ 2x,1p3e22.14)) cbug 9902 format (' inv=',i2, cbug & ' (0 for direct, 1 for transpose operation.)' / cbug & 'invariant pt ',1p3e22.14) cbug 9903 format (i3,' px,py,pz=',1p3e22.14) cbug write ( 3, 9901) np, tol, ((smat(i,j), j = 1, 3), i = 1, 3) cbug write ( 3, 9902) inv, ax, ay, az cbug write ( 3, 9903) (n, px(n), py(n), pz(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 cbugc***DEBUG begins. cbug write ( 3, '(/ "aptmopv fatal. bad np=",i3)') np cbugc***DEBUG ends. go to 210 endif if ((inv .lt. 0) .or. (inv .gt. 1)) then nerr = 2 cbugc***DEBUG begins. cbug write ( 3, '(/ "aptmopv fatal. bad inv=",i3)') inv cbugc***DEBUG ends. go to 210 endif c.... Translate the origin to point "a". call apttran (ax, ay, az, px, py, pz, np, tol, nerr) c.... Operate with smat or its transpose. c---- Operate with smat. if (inv .eq. 0) then c---- No truncation tests. if (tol .le. 0.0) then c---- Loop over points or vectors. do 110 n = 1, np ppx = px(n) ppy = py(n) ppz = pz(n) px(n) = smat(1,1) * ppx + smat(1,2) * ppy + smat(1,3) * ppz py(n) = smat(2,1) * ppx + smat(2,2) * ppy + smat(2,3) * ppz pz(n) = smat(3,1) * ppx + smat(3,2) * ppy + smat(3,3) * ppz c---- End of loop over points or vectors. 110 continue c---- Truncate small components to zero. else c---- Loop over points or vectors. do 120 n = 1, np ppx = px(n) ppy = py(n) ppz = pz(n) px(n) = smat(1,1) * ppx + smat(1,2) * ppy + smat(1,3) * ppz py(n) = smat(2,1) * ppx + smat(2,2) * ppy + smat(2,3) * ppz pz(n) = smat(3,1) * ppx + smat(3,2) * ppy + smat(3,3) * ppz pxerr = tol * (abs (smat(1,1) * ppx) + & abs (smat(1,2) * ppy) + & abs (smat(1,3) * ppz)) pyerr = tol * (abs (smat(2,1) * ppx) + & abs (smat(2,2) * ppy) + & abs (smat(2,3) * ppz)) pzerr = tol * (abs (smat(3,1) * ppx) + & abs (smat(3,2) * ppy) + & abs (smat(3,3) * ppz)) 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 c---- Operate with transpose of smat. else c---- No truncation tests. if (tol .le. 0.0) then c---- Loop over points or vectors. do 130 n = 1, np ppx = px(n) ppy = py(n) ppz = pz(n) px(n) = smat(1,1) * ppx + smat(2,1) * ppy + smat(3,1) * ppz py(n) = smat(1,2) * ppx + smat(2,2) * ppy + smat(3,2) * ppz pz(n) = smat(1,3) * ppx + smat(2,3) * ppy + smat(3,3) * ppz c---- End of loop over points or vectors. 130 continue c---- Truncate small components to zero. else c---- Loop over points or vectors. do 140 n = 1, np ppx = px(n) ppy = py(n) ppz = pz(n) px(n) = smat(1,1) * ppx + smat(2,1) * ppy + smat(3,1) * ppz py(n) = smat(1,2) * ppx + smat(2,2) * ppy + smat(3,2) * ppz pz(n) = smat(1,3) * ppx + smat(2,3) * ppy + smat(3,3) * ppz pxerr = tol * (abs (smat(1,1) * ppx) + & abs (smat(2,1) * ppy) + & abs (smat(3,1) * ppz)) pyerr = tol * (abs (smat(1,2) * ppx) + & abs (smat(2,2) * ppy) + & abs (smat(3,2) * ppz)) pzerr = tol * (abs (smat(1,3) * ppx) + & abs (smat(2,3) * ppy) + & abs (smat(3,3) * ppz)) 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. 140 continue c---- Tested tol. endif c---- Tested inv. endif c.... Translate the origin from point "a" back to its original position. call apttran (-ax, -ay, -az, px, py, pz, np, tol, nerr) cbugc***DEBUG begins. cbug 9904 format (/ 'aptmopv results:' / cbug & (i3,' px,py,pz=',1p3e22.14)) cbug write ( 3, 9904) (n, px(n), py(n), pz(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptmopv. (+1 line.) end UCRL-WEB-209832