subroutine aptrota (theta, ku, ax, ay, az, tol, rotm, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTROTA
c
c call aptrota (theta, ku, ax, ay, az, tol, rotm, nerr)
c
c Version: aptrota Updated 1990 March 14 16:00.
c aptrota 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
c counterclockwise rotation by the angle theta, around the
c axis parallel to the vector a = (ax, ay, az).
c (Counterclockwise: with the axis pointed at the observer.)
c The angle may be in degrees (ku = 0) or radians (ku = 1).
c Flag nerr indicates any error in vector "a" (too small)
c or ku (not 0, 1).
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c
c Input: theta, ku, ax, ay, az.
c
c Output: rotm, nerr.
c
c Calls: aptvunb
c
c Glossary:
c
c ax,ay,az Input The x, y, z components of the vector parallel to the
c rotation axis.
c
c ku Input Indicates theta units are degrees (0) or radians (1).
c
c nerr Output Indicates an input error, if not 0.
c 1 if magnitude of (ax, ay, az) too small.
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 theta Input Angle of rotation around axis "a", counterclockwise
c when rotation axis is pointed at observer. Units
c are degrees (ku = 0) or radians (ku = 1).
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 in radians.
common /laptrota/ angle
c---- Cosine of rotation angle.
common /laptrota/ costh
c---- 1.0 - cosine of rotation angle.
common /laptrota/ costhm
c---- Index of matrix row.
common /laptrota/ i
c---- Index of matrix column.
common /laptrota/ j
c---- Mathematical constant, pi.
common /laptrota/ pi
c---- Sine of rotation angle.
common /laptrota/ sinth
c---- Component x of unit axis vector.
common /laptrota/ ux
c---- Component y of unit axis vector.
common /laptrota/ uy
c---- Component z of unit axis vector.
common /laptrota/ uz
c---- Magnitude of a vector.
common /laptrota/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrota. Rotate by angle ',f12.6,' (ku=',i2,')' /
cbug & ' around axis (',1p3e13.5,').')
cbug write ( 3, 9901) theta, ku, ax, ay, az
cbugc***DEBUG ends.
c.... Initialize.
c---- Mathematical constant, pi.
pi = 3.14159265358979323
nerr = 0
c.... Normalize rotation axis vector.
call aptvunb (ax, ay, az, 1, 0., ux, uy, uz, vlen, nerr)
if (vlen .le. tol) then
nerr = 1
cbugc***DEBUG begins.
cbug write ( 3, '(/ "aptrota fatal. bad a length=",1pe22.14)') vlen
cbugc***DEBUG ends.
go to 210
endif
c.... Get rotation angle in radians, and its sine and cosine.
if (ku .eq. 0) then
angle = theta * pi / 180.0
elseif (ku .eq. 1) then
angle = theta
else
nerr = 2
cbugc***DEBUG begins.
cbug write ( 3, '(/ "aptrota fatal. bad ku=",i3)') ku
cbugc***DEBUG ends.
go to 210
endif
sinth = sin (angle)
costh = cos (angle)
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
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptrota results. rotm=',3(/ 11x,1p3e22.14))
cbug write ( 3, 9902) ((rotm(i,j), j = 1, 3), i = 1, 3)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptrota. (+1 line.)
end
UCRL-WEB-209832