subroutine aptrots (n1, th1, n2, th2, n3, th3, ku, tol, & rotm, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTROTS c c call aptrots (n1, th1, n2, th2, n3, th3, ku, tol, c & rotm, nerr) c c Version: aptrots Updated 1989 November 13 15:20. c aptrots 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, to do a sequential c rotation of angle th1 around axis n1, angle th2 around axis n2, c and angle th3 around axis n3. All angles are measured c counterclockwise, with the axis pointed at the observer. c The axes may be x (1), y (2), or z (3). c Angles may be in degrees (ku = 0) or radians (ku = 1). c Flag nerr indicates any input error. c c Input: n1, th1, n2, th2, n3, th3, ku, tol. c c Output: rotm, nerr. c c Glossary: c c ku Input Indicates theta units are degrees (0) or radians (1). c c n1 Input Indicates first axis is x (1), y (2), or z (3). c May not be 0, but th1 may be 0. c c n2 Input Indicates second axis is x (1), y (2), or z (3). c May not be 0, but th2 may be 0. c c n3 Input Indicates third axis is x (1), y (2), or z (3). c May not be 0, but th3 may be 0. c c nerr Output Indicates an input error, if not 0. c 1 if n1, n2, or n3 not in range 1-3, or not unique. c 2 if ku is not 0 or 1. c c rotm Output Rotation operator (a unitary 3 x 3 matrix). c Must be sized rotm(3,3). c c th1 Input Angle of rotation around axis n1, counterclockwise c when rotation axis is pointed at observer. Units c are degrees (ku = 0) or radians (ku = 1). c c th2 Input Angle of rotation around axis n2. See th1. c c th3 Input Angle of rotation around axis n3. See th1. 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---- Rotation angle around axis. common /laptrots/ ang (3) c---- Cosine of rotation angle ang. common /laptrots/ costh (3) c---- Index of matrix row. common /laptrots/ i c---- Index of matrix column. common /laptrots/ j c---- Parity of axis sequence. common /laptrots/ npar c---- Sum of axis numbers. common /laptrots/ nsum c---- Parity of axis sequence. common /laptrots/ par c---- Mathematical constant, pi. common /laptrots/ pi c---- Sine of rotation angle ang. common /laptrots/ sinth(3) cbugc***DEBUG begins. cbug 9901 format (/ 'aptrots. Rotate sequentially around axis by angle:' / cbug & 5x,3(i2,f13.6)) cbug write ( 3, 9901) n1, th1, n2, th2, n3, th3 cbugc***DEBUG ends. c.... Initialize. c---- Mathematical constant, pi. pi = 3.14159265358979323 nerr = 0 c.... Find parity of sequence of axes. nsum = n1 + n2 + n3 c---- Bad axis numbers. if (nsum .ne. 6) then nerr = 1 go to 210 endif npar = (n2 - n1) * (n3 - n2) * (n3 - n1) c---- 1,2,3 or 2,3,1 or 3,1,2. if (npar .eq. 2) then par = 1.0 c---- 3,2,1 or 2,1,3 or 1,3,2. elseif (npar .eq. -2) then par = -1.0 else nerr = 1 go to 210 endif c.... Get rotation angles in radians, and the sines and cosines. if (ku .eq. 0) then ang(n1) = th1 * pi / 180.0 ang(n2) = th2 * pi / 180.0 ang(n3) = th3 * pi / 180.0 elseif (ku .eq. 1) then ang(n1) = th1 ang(n2) = th2 ang(n3) = th3 else nerr = 2 go to 210 endif costh(n1) = cos (ang(n1)) costh(n2) = cos (ang(n2)) costh(n3) = cos (ang(n3)) sinth(n1) = sin (ang(n1)) sinth(n2) = sin (ang(n2)) sinth(n3) = sin (ang(n3)) c.... Form the components of the rotation matrix. rotm(n1,n1) = costh(n2) * costh(n3) rotm(n1,n2) = sinth(n1) * sinth(n2) * costh(n3) - & par * costh(n1) * sinth(n3) rotm(n1,n3) = sinth(n1) * sinth(n3) + & par * costh(n1) * sinth(n2) * costh(n3) rotm(n2,n1) = par * costh(n2) * sinth(n3) rotm(n2,n2) = costh(n1) * costh(n3) + & par * sinth(n1) * sinth(n2) * sinth(n3) rotm(n2,n3) = costh(n1) * sinth(n2) * sinth(n3) - & par * sinth(n1) * costh(n3) rotm(n3,n1) = -par * sinth(n2) rotm(n3,n2) = par * sinth(n1) * costh(n2) rotm(n3,n3) = costh(n1) * costh(n2) 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 return c.... End of subroutine aptrots. (+1 line.) end UCRL-WEB-209832