subroutine aptrotv (ax, ay, az, bx, by, bz, tol, rotm, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTROTV c c call aptrotv (ax, ay, az, bx, by, bz, tol, rotm, nerr) c c Version: aptrotv Updated 1991 July 25 14:30. c aptrotv 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 rotation matrix operator "rotm", for rotating c the vector a = (ax, ay, az) to be parallel to the vector c b = (bx, by, bz), around an axis perpendicular to both. c Any components of "rotm" within tol of -1.0, 0.0, or 1.0, c will be truncated to those values. c Flag nerr indicates any input error. c c History: 1990 March 14. Changed tol to 0.0 in call to unit vector c module aptvunb. Allows small magnitudes. c Protected costh from absolute values exceeding 1.0. c 1991 July 25. Added coding to find an arbitrary rotation c axis perpendicular to vectors "a" and "b" when they are c parallel or antiparallel. Previously, the resulting operator c "rotm" was an inversion operator, not a rotation operator. c c Input: ax, ay, az, bx, by, bz, tol. c c Output: rotm, nerr. c c Calls: aptvxun, aptvdot, aptvunb c c c c Glossary: c c ax,ay,az Input The x, y, z components of vector "a". c c bx,by,bz Input The x, y, z components of vector "b". c c nerr Output Indicates an input error, if not 0. c 1 if the magnitude of vector "a" or "b" is too small. c c rotm Output Rotation operator (a unitary 3 x 3 matrix). c Must be sized rotm(3,3). c c tol Input Numerical tolerance limit. Used to test and adjust c unit vector components and point coordinates. c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Rotation matrix operator. dimension rotm (3,3) c.... Local variables. c---- Cosine of rotation angle. common /laptrotv/ costh c---- 1.0 - cosine of rotation angle. common /laptrotv/ costhm c---- Index of matrix row. common /laptrotv/ i c---- Index of matrix column. common /laptrotv/ j c---- Sine of rotation angle. common /laptrotv/ sinth c---- Component x of vector "ua". common /laptrotv/ uax c---- Component y of vector "ua". common /laptrotv/ uay c---- Component z of vector "ua". common /laptrotv/ uaz c---- Component x of vector "ub". common /laptrotv/ ubx c---- Component y of vector "ub". common /laptrotv/ uby c---- Component z of vector "ub". common /laptrotv/ ubz c---- Component x of unit axis vector. common /laptrotv/ ux c---- Component y of unit axis vector. common /laptrotv/ uy c---- Component z of unit axis vector. common /laptrotv/ uz c---- Magnitude of a vector. common /laptrotv/ vlen cbugc***DEBUG begins. cbug 9901 format (/ 'aptrotv. tol=',1pe13.5,' Rotating vector' / cbug & ' ax,ay,az=',1p3e22.14 / cbug & ' around mutually perpendicular axis onto vector' / cbug & ' bx,by,bz=',1p3e22.14) cbug write ( 3, 9901) tol, ax, ay, az, bx, by, bz cbugc***DEBUG ends. c.... Initialize. nerr = 0 c.... Find the unit vector parallel to "a". call aptvunb (ax, ay, az, 1, 0., uax, uay, uaz, vlen, nerr) c---- Vector magnitude too small. if (vlen .le. tol) then nerr = 1 go to 210 endif c.... Find the unit vector parallel to "b". call aptvunb (bx, by, bz, 1, 0., ubx, uby, ubz, vlen, nerr) c---- Vector magnitude too small. if (vlen .le. tol) then nerr = 1 go to 210 endif c.... Find the unit vector parallel to the axis of rotation. call aptvxun (uax, uay, uaz, ubx, uby, ubz, 1, tol, & ux, uy, uz, vlen, nerr) c.... Find an arbitrary rotation axis if "a" and "b" are parallel c.... or antiparallel. c---- Vectors "a" and "b" aligned. if (vlen .eq. 0.0) then if (uax .ne. 0.0) then vlen = sqrt (uax**2 + uay**2) ux = -uay / vlen uy = uax / vlen elseif (uay .ne. 0.0) then vlen = sqrt (uay**2 + uaz**2) uy = -uaz / vlen uz = uay / vlen elseif (uaz .ne. 0.0) then vlen = sqrt (uaz**2 + uax**2) uz = -uax / vlen ux = uaz / vlen endif c---- Tested rotation axis vector. endif c.... Find the needed functions of the angle of rotation. c.... costh = ua .dot. ub (scalar product). call aptvdot (uax, uay, uaz, ubx, uby, ubz, 1, tol, costh, nerr) costh = amax1 (-1.0, amin1 (costh, 1.0)) costhm = 1.0 - costh sinth = sqrt (costhm * (1.0 + costh)) c.... Form the components of the rotation matrix. rotm(1,1) = costhm * ux * ux + costh rotm(1,2) = costhm * ux * uy - sinth * uz rotm(1,3) = costhm * ux * uz + sinth * uy rotm(2,1) = costhm * uy * ux + sinth * uz rotm(2,2) = costhm * uy * uy + costh rotm(2,3) = costhm * uy * uz - sinth * ux rotm(3,1) = costhm * uz * ux - sinth * uy rotm(3,2) = costhm * uz * uy + sinth * ux rotm(3,3) = costhm * uz * uz + costh c---- Adjust components near 0 or 1. if (tol .gt. 0.0) then do 120 i = 1, 3 do 110 j = 1, 3 if (abs (rotm(i,j)) .le. tol) then rotm(i,j) = 0.0 elseif ((abs (abs (rotm(i,j)) - 1.0)) .le. tol) then rotm(i,j) = sign (1.0, rotm(i,j)) endif 110 continue 120 continue c---- Tested tol. endif 210 continue cbugc***DEBUG begins. cbug 9902 format (/ 'aptrotv results: rotm=',3(/ 11x,1p3e22.14)) cbug write ( 3, 9902) ((rotm(i,j), j = 1, 3), i = 1, 3) cbugc***DEBUG ends. return c.... End of subroutine aptrotv. (+1 line.) end UCRL-WEB-209832