subroutine aptrotp (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, tol, rotm, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTROTP
c
c call aptrotp (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, tol, rotm, nerr)
c
c Version: aptrotp Updated 1990 October 29 17:45
c aptrotp 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 the
c plane containing the vectors a = (ax, ay, az) and
c b = (bx, by, bz) to be parallel to the plane containing the
c vectors c = (cx, cy, cz) and d = (dx, dy, dz), around an axis
c parallel to both planes. Any components of rotm within tol
c of -1.0, 0.0, or 1.0, 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 subroutine. Allows small magnitudes.
c Protected costh from absolute values exceeding 1.0.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz.
c
c Output: rotm, nerr.
c
c Calls: aptvaxb, aptvuna, aptvxun
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z components of a vector.
c
c bx,by,bz Input The x, y, z components of a vector.
c
c cx,cy,cz Input The x, y, z components of a vector.
c
c dx,dy,dz Input The x, y, z components of a vector.
c
c nerr Output Indicates an input error, if not 0.
c 1 if the magnitude of any input vector is too small,
c or the two vectors in a plane are almost parallel.
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 /laptrotp/ costh
c---- 1.0 - cosine of rotation angle.
common /laptrotp/ costhm
c---- Index of matrix row.
common /laptrotp/ i
c---- Index of matrix column.
common /laptrotp/ j
c---- Sine of rotation angle.
common /laptrotp/ sinth
c---- Component x of normal vector.
common /laptrotp/ uabx
c---- Component y of normal vector.
common /laptrotp/ uaby
c---- Component z of normal vector.
common /laptrotp/ uabz
c---- Component x of normal vector.
common /laptrotp/ ucdx
c---- Component y of normal vector.
common /laptrotp/ ucdy
c---- Component z of normal vector.
common /laptrotp/ ucdz
c---- Component x of unit axis vector.
common /laptrotp/ ux
c---- Component y of unit axis vector.
common /laptrotp/ uy
c---- Component z of unit axis vector.
common /laptrotp/ uz
c---- Magnitude of a vector.
common /laptrotp/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrotp. Rotating plane containing vectors' /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' to be parallel to plane containing vectors' /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' dx,dy,dz=',1p3e22.14)
cbug write ( 3, 9901) ax, ay, az, bx, by, bz, cx, cy, cz,
cbug & dx, dy, dz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Find the unit vector normal to the first plane.
call aptvxun (ax, ay, az, bx, by, bz, 1, tol,
& uabx, uaby, uabz, vlen, nerr)
c---- Vector magnitude too small.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
c.... Find the unit vector normal to the second plane.
call aptvxun (cx, cy, cz, dx, dy, dz, 1, tol,
& ucdx, ucdy, ucdz, 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 aptvaxb (uabx, uaby, uabz, ucdx, ucdy, ucdz, 1, tol,
& ux, uy, uz, vlen, nerr)
c---- Planes are parallel.
if (vlen .le. tol) then
vlen = 1.0
ux = 1.0
uy = 0.0
uz = 0.0
else
call aptvuna (ux, uy, uz, 1, 0., vlen, nerr)
endif
c.... Find the needed functions of the angle of rotation.
costh = uabx * ucdx + uaby * ucdy + uabz * ucdz
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
cbugc***DEBUG begins.
cbug 9902 format (' 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 aptrotp. (+1 line.)
end
UCRL-WEB-209832