subroutine aptrott (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, tol, rotm, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTROTT
c
c call aptrott (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, tol, rotm, nerr)
c
c Version: aptrott Updated 1990 March 14 16:00.
c aptrott 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 vector a = (ax, ay, az) and the plane containing vectors "a" and
c b = (bx, by, bz), to be parallel to vector c = (cx, cy, cz) and
c the plane containing vectors "c" and d = (dx, dy, dz).
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 If vectors "a" and "b" are the first two vectors of the positive
c vector triple (a, b, a x b), and c and d are the first two
c vectors of the positive vector triple (c, d, c x d), then rotm
c rotates (a, b, a x b) onto (c, d, c x d), or equivalently,
c redefines the coordinate axes to be (a, b, c x d) instead of
c (c, d, c x d). (a x b indicates the vector product of a and b.)
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, tol.
c
c Output: rotm, nerr.
c
c Calls: aptvxun, aptvunb
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---- Index of matrix row.
common /laptrott/ i
c---- Index of matrix column.
common /laptrott/ j
c---- Component x of unit vector "uab".
common /laptrott/ uabx
c---- Component y of unit vector "uab".
common /laptrott/ uaby
c---- Component z of unit vector "uab".
common /laptrott/ uabz
c---- Component x of unit vector "ua".
common /laptrott/ uax
c---- Component y of unit vector "ua".
common /laptrott/ uay
c---- Component z of unit vector "ua".
common /laptrott/ uaz
c---- Component x of unit vector "ub".
common /laptrott/ ubx
c---- Component y of unit vector "ub".
common /laptrott/ uby
c---- Component z of unit vector "ub".
common /laptrott/ ubz
c---- Component x of unit vector "ucd".
common /laptrott/ ucdx
c---- Component y of unit vector "ucd".
common /laptrott/ ucdy
c---- Component z of unit vector "ucd".
common /laptrott/ ucdz
c---- Component x of unit vector "uc".
common /laptrott/ ucx
c---- Component y of unit vector "uc".
common /laptrott/ ucy
c---- Component z of unit vector "uc".
common /laptrott/ ucz
c---- Component x of unit vector "ud".
common /laptrott/ udx
c---- Component y of unit vector "ud".
common /laptrott/ udy
c---- Component z of unit vector "ud".
common /laptrott/ udz
c---- Magnitude of a vector.
common /laptrott/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrott. 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 & ' with vector a parallel to vector c.')
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 positive unit vector triple "ua", "ub", "uab".
c.... Find the unit vector parallel to vector "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 vector "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 perpendicular to vectors "a" and "b".
call aptvxun (uax, uay, uaz, ubx, uby, ubz, 1, tol,
& uabx, uaby, uabz, vlen, nerr)
c---- Vectors almost parallel.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
c.... Replace "ub" with the unit vector perpendicular to vectors "uab" and "a".
call aptvxun (uabx, uaby, uabz, uax, uay, uaz, 1, tol,
& ubx, uby, ubz, vlen, nerr)
c.... Find the positive unit vector triple "uc", "ud", "ucd".
c.... Find the unit vector parallel to vector "c".
call aptvunb (cx, cy, cz, 1, 0., ucx, ucy, ucz, 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 vector "d".
call aptvunb (dx, dy, dz, 1, 0., udx, udy, udz, vlen, nerr)
c---- Vector magnitude too small.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
c.... Find the unit vector perpendicular to vectors "c" and "d".
call aptvxun (ucx, ucy, ucz, udx, udy, udz, 1, tol,
& ucdx, ucdy, ucdz, vlen, nerr)
c---- Vectors almost parallel.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
c.... Replace "ud" with the unit vector perpendicular to vectors "ucd" and "c".
call aptvxun (ucdx, ucdy, ucdz, ucx, ucy, ucz, 1, tol,
& udx, udy, udz, vlen, nerr)
c.... Find the components of the rotation matrix.
rotm(1,1) = uax * ucx + ubx * udx + uabx * ucdx
rotm(1,2) = uay * ucx + uby * udx + uaby * ucdx
rotm(1,3) = uaz * ucx + ubz * udx + uabz * ucdx
rotm(2,1) = uax * ucy + ubx * udy + uabx * ucdy
rotm(2,2) = uay * ucy + uby * udy + uaby * ucdy
rotm(2,3) = uaz * ucy + ubz * udy + uabz * ucdy
rotm(3,1) = uax * ucz + ubx * udz + uabx * ucdz
rotm(3,2) = uay * ucz + uby * udz + uaby * ucdz
rotm(3,3) = uaz * ucz + ubz * udz + uabz * ucdz
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 aptrott. (+1 line.)
end
UCRL-WEB-209832