Rotation 
Table of ContentsRotation Abstract Introduction Using this document References Notation SPICELIB Routines Categories of routines Euler angle routines Constructing a matrix from Euler angles Finding Euler angles that represent a matrix Programming hazards Working with RA, Dec and Twist Quaternions Finding a quaternion that represents a matrix Finding the matrix represented by a quaternion Q2M and M2Q are approximate inverses of each other Multiplying quaternions Obtaining angular velocity from quaternions Rotating vectors and matrices A word of warning Rotating a vector about a coordinate axis Rotating a matrix about a coordinate axis Rotating a vector about an arbitrary axis Rotation axis and angle Constructing a matrix from a rotation axis and angle Finding the axis and angle of a rotation matrix AXISAR and RAXISA are approximate inverses Using RAXISA and AXISAR Constructing a coordinate axis rotation matrix Rotation derivatives Differentiating rotations State transformations Validating a rotation matrix Tutorial introduction to rotations A comment of the heuristic variety Definition of ``rotation'' Definition 1 Definition 2 Definition 3 Uses of the definitions Definition of ``rotation'' and ``orthogonal'' matrix Rotations preserve inner products Inverses of rotation matrices Composition of rotations Coordinate transformations Rotation of vectors in the plane A canonical representation for rotations Rotation axis and angle Timedependent coordinate transformations Euler angles Quaternions Quaternion arithmetic Definitions Basic properties of multiplication Deducing the multiplication formula Composing rotations using quaternions Mathematical road map Rotation of a vector about an axis Formation of a rotation matrix from axis and angle Finding the axis and angle of a rotation matrix Formation of a rotation matrix from a quaternion Equivalence of rotation definitions An algebraic approach A geometric approach Quaternion multiplication Assertion 1 Assertion 2 Recovery of Euler angles from a rotation matrix Euler angle recovery: aba case Euler angle recovery: abc case Appendix A: Document Revision History March 9, 2017 May 27, 2010 November 17, 2005 January 10, 2005 February 2, 2004 December 2, 2002 April 26, 1999 Rotation
Abstract
Introduction
There are three chapters:
Using this document
The rest of the document covers the ideas behind the code. This material is meant to be used as a reference rather than to be read from start to finish; the topics are ordered loosely according to logical dependence, but there is no narrative progression from section to section. The purpose of the tutorial and ``road map'' chapters is to make it easier to be certain that you're using the code correctly. In our experience, thinking about this category of code only in terms of ``inputs'' and ``outputs'' is a tricky and errorprone approach; really understanding the mathematics helps you to verify that you're putting the pieces together in ways that make sense. Because some of the ideas required to understand the code seem to exist as lore and are rarely written down, we've collected them here for your convenience. References
Notation
Symbol Meaning   E3 Threedimensional Euclidean space. E2 Twodimensional Euclidean space. < u, v > Inner product of vectors u and v. T M Transpose of the matrix M. 1 M Inverse of the matrix M. 1 f Inverse of the function f. u x v Cross product of vectors u and v. M v Product of matrix M and vector v. M N Product of matrix M and matrix N. Trace(M) Sum of elements on the main diagonal of M. [w] Matrix that rotates a coordinate system by i w radians about the ith coordinate axis (and rotates vectors by w radians about the same axis). We also use this notation to refer to the linear transformation corresponding to this matrix.There are only three types of rotation matrices representing rotations about coordinate axes:
+ +  1 0 0  [w] =  0 cos(w) sin(w)  1  0 sin(w) cos(w)  + + + +  cos(w) 0 sin(w)  [w] =  0 1 0  2  sin(w) 0 cos(w)  + + + +  cos(w) sin(w) 0  [w] =  sin(w) cos(w) 0  3  0 0 1  + +The composition of rotations
[w1] [w2] [w3] a b cis sometimes referred to as a ``abc'' rotation. For example, we may talk about a ``313'' rotation. About angles: all angles in this document are measured in radians. About directions: the ``right hand rule'' is in effect at all times in this document, so counterclockwise rotations about an axis have positive measure. SPICELIB Routines
The SPICELIB routines that deal with general rotations are:
Categories of routines
Euler angle routines
Constructing a matrix from Euler angles
The call
CALL EUL2M ( ANG(3), ANG(2), ANG(1), . I(3), I(2), I(1), M )returns the matrix M, where
M = [ANG(3)] [ANG(2)] [ANG(1)] . I(3) I(2) I(1)The indices I(3), I(2), I(1) must belong to the set
{1, 2, 3}. Finding Euler angles that represent a matrix
M = [ANG(3)] [ANG(2)] [ANG(1)] . I(3) I(2) I(1)The call
CALL M2EUL ( M, I(3), I(2), I(1), . ANG(3), ANG(2), ANG(1) )returns the desired angles. M2EUL restricts the ranges of the output angles so as to guarantee that the Euler angle representation is unique. The output angles ANG(3) and ANG(1) are always in the range (pi,pi]. The range of ANG(2) is determined by the set of rotation axes. When I(3) equals I(1), ANG(2) is in the range [0, pi]. Otherwise, ANG(2) is in the range [pi/2, pi/2]. These ranges make unique determinations of Euler angles possible, except in degenerate cases. In cases where the Euler angles are not uniquely determined, EUL2M sets the first angle (called ANG(3) above) to zero. The other two angles are then uniquely determined. Again, the indices I(3), I(2), I(1) are members of
{1, 2, 3}.There is a restriction on the allowed set of coordinate axes: I(2) must not equal I(3) or I(1). If this constraint is not met, the desired representation of M by Euler angles may not exist; M2EUL signals an error in this case. Programming hazards
The first category of singularity occurs with matrices that represent rotations about the first or third axis in the sequence of rotation axes (for example, axis 3 for a 213 rotation). In practical terms, if EUL2M encounters one of these special matrices, EUL2M must choose the Euler angles. Immediately the possibility arises that EUL2M will disagree with any other code performing the same task. The second kind of singularity occurs when any of the Euler angles corresponding to a matrix is at one of the endpoints of its range, for example, when the first angle has the value pi. If the matrix is perturbed slightly, the first angle may jump from pi to a value close to pi. Again, two different pieces of code may give different results in such a case, merely because of roundoff error. Euler angles near the limits of their ranges should be regarded with suspicion. The existence of singularities in the matrixtoEuler angle mapping prevents EUL2M and M2EUL from being exact inverses: most of the time, the code fragment
CALL EUL2M ( ANG(3), ANG(2), ANG(1), . AXIS(1), AXIS(2), AXIS(3), M ) CALL M2EUL ( M, AXIS(1), AXIS(2), AXIS(3), . ANG(3), ANG(2), ANG(1) )leaves the angles ANG(3), ANG(2), ANG(1) unchanged, except for roundoff error, but in some cases, the angles may change drastically. If we reverse the order of the routine calls in the last code fragment, the matrix M should be preserved, except for errors due to loss of precision. The loss of precision can be considerable, though, for matrices whose entries are nearly those of any degenerate case matrix. For more details on this topic, consult the ``Mathematical road map'' section. Working with RA, Dec and Twist
C = [ Twist ] [ pi/2  Dec ] [ pi/2 + RA ] 3 1 3so the code fragments
CALL EUL2M ( TWIST, HALFPI()  DEC, HALFPI() + RA . 3, 1, 3, C )and
CALL M2EUL ( C, 3, 1, 3, . ANG(3), ANG(2), ANG(1) ) TWIST = ANG(3) DEC = HALFPI()  ANG(2) RA = ANG(1)  HALFPI()carry out the conversion from RA, Dec, and Twist to a C matrix, and back. Note that definitions of ``RA, Dec, and Twist'' vary: on the Galileo project, the C matrix is related to the angles ``RA, Dec, and Twist'' by the equation
C = [ Twist ] [ pi/2  Dec ] [ RA ] . 3 2 3 Quaternions
Finding a quaternion that represents a matrix
The code fragment
CALL M2Q ( M, Q )returns a unit quaternion Q that represents the rotation matrix M. If you really want to know about it, the elements of the quaternion are defined as follows: Let the unit vector A be a choice of rotation axis for M, and let ANGLE be the rotation angle, where ANGLE is in the interval
[0, pi]Then the elements of the unit quaternion Q returned by M2Q are
Q(0) = cos( ANGLE/2 ) Q(1) = sin( ANGLE/2 ) * A(1) Q(2) = sin( ANGLE/2 ) * A(2) Q(3) = sin( ANGLE/2 ) * A(3)M2Q considers the rotation angle of M to lie in [0, pi]. Therefore, half the rotation angle lies in [0, pi/2], so Q(0) is always in [0, 1]. For a given rotation matrix M, the corresponding quaternion is uniquely determined except if the rotation angle is pi. Finding the matrix represented by a quaternion
CALL Q2M ( Q, M ) Q2M and M2Q are approximate inverses of each other
CALL M2Q ( M, Q ) CALL Q2M ( Q, M )always preserve M, except for roundoff error. However, since there are two quaternions that represent each rotation, the sequence of calls
CALL Q2M ( Q, M ) CALL M2Q ( M, Q )do not necessarily preserve Q. Multiplying quaternions
CALL QXQ ( Q1, Q2, QOUT )The resulting product QOUT is computed using the multiplication formula given in the section ``Quaternion Arithmetic'' below. QOUT represents the rotation formed by composing the rotations represented by Q1 and Q2. Obtaining angular velocity from quaternions
CALL QDQ2AV ( Q, DQ, AV )The resulting angular velocity vector AV has units of radians/T, where 1/T is the time unit associated with DQ. Rotating vectors and matrices
A word of warning
Rotating a vector about a coordinate axis
[angle] ito a vector, use ROTVEC. The following code fragment applies
[ ANGLE ] 3to the vector V, yielding VOUT.
CALL ROTVEC ( V, ANGLE, 3, VOUT )The components of VOUT are the coordinates of the vector V in a system rotated by ANGLE radians about the third coordinate axis. We can also regard VOUT as V, rotated by
ANGLEradians about the third coordinate axis. Rotating a matrix about a coordinate axis
[angle] ito a matrix using ROTMAT, as follows:
CALL ROTMAT ( M, ANGLE, 3, MOUT )After this routine call, MOUT is equal to
[ANGLE] * M 3 Rotating a vector about an arbitrary axis
CALL VROTV ( V, AXIS, ANGLE, VOUT )rotates the vector V about the vector AXIS by ANGLE radians, yielding VOUT. Rotation axis and angle
Constructing a matrix from a rotation axis and angle
CALL AXISAR ( AXIS, ANGLE, M )produces M, the desired rotation matrix. What if we want generate a coordinate system rotation about an arbitrary axis, as opposed to a coordinate axis? We can use AXISAR for this. Let AXIS be the coordinate system rotation axis and ANGLE be the rotation angle; then the code fragment
CALL AXISAR ( AXIS, ANGLE, M )produces the desired coordinate system rotation matrix. Note that the input angle is the NEGATIVE of that associated with a vector rotation. AXISAR is designed this way for compatibility with RAXISA, which is an inverse routine for AXISAR. Finding the axis and angle of a rotation matrix
CALL RAXISA ( M, AXIS, ANGLE )AXIS and ANGLE have the property that for any vector V,
M Vyields V, rotated by ANGLE radians about the vector AXIS. If M is viewed as a coordinate transformation, we can say that M rotates the initial coordinate system by
 ANGLEradians about AXIS. AXISAR and RAXISA are approximate inverses
CALL RAXISA ( M, AXIS, ANGLE ) CALL AXISAR ( AXIS, ANGLE, M )leaves M unchanged, except for roundoff error. The code fragment
CALL AXISAR ( AXIS, ANGLE, M ) CALL RAXISA ( M, AXIS, ANGLE )usually leaves AXIS and ANGLE unchanged, except for roundoff error, provided that two conditions are met:
Using RAXISA and AXISAR
C C Let R(t) be a timevarying rotation matrix; R could be C a Cmatrix describing the orientation of a spacecraft C structure. Given two points in time t1 and t2 at which C R(t) is known, and given a third time t3, where C C t1 < t3 < t2, C C we can estimate R(t3) by linear interpolation. In other C words, we approximate the motion of R by pretending that C R rotates about a fixed axis at a uniform angular rate C during the time interval [t1, t2]. More specifically, we C assume that each column vector of R rotates in this C fashion. This procedure will not work if R rotates through C an angle of pi radians or more during the time interval C [t1, t2]; an aliasing effect would occur in that case. C C If we let C C R1 = R(t1) C R2 = R(t2), and C C 1 C Q = R2 * R1 , C C then the rotation axis and angle of Q define the rotation C that each column of R(t) undergoes from time t1 to time C t2. Since R(t) is orthogonal, we can find Q using the C transpose of R1. We find the rotation axis and angle via C RAXISA. C CALL MXMT ( R2, R1, Q ) CALL RAXISA ( Q, AXIS, ANGLE ) C C Find the fraction of the total rotation angle that R C rotates through in the time interval [t1, t3]. C FRAC = ( T3  T1 ) / ( T2  T1 ) C C Finally, find the rotation DELTA that R(t) undergoes C during the time interval [t1, t3], and apply that rotation C to R1, yielding R(t3), which we'll call R3. C CALL AXISAR ( AXIS, FRAC * ANGLE, DELTA ) CALL MXM ( DELTA, R1, R3 ) Constructing a coordinate axis rotation matrix
[ANGLE] , Iwhich corresponds to a rotation about a coordinate axis. This is a special case of the problem solved by AXISAR. Note however that the matrix produced by ROTATE is the inverse of that produced by AXISAR, if both routines are provided with the same input angle, and AXISAR is given the Ith coordinate basis vector as the rotation axis. The call
CALL ROTATE ( ANGLE, I, M )produces M, the desired matrix. Rotation derivatives
Differentiating rotations
[ANGLE] 1which has the matrix representation
+ +  1 0 0   0 cos(ANGLE) sin(ANGLE)   0 sin(ANGLE) cos(ANGLE)  + +yields, when differentiated with respect to ANGLE, the matrix
+ +  0 0 0   0 sin(ANGLE) cos(ANGLE)   0 cos(ANGLE) sin(ANGLE)  + +The routine DROTAT is useful for differentiating rotations that are defined by timevarying Euler angles. For example, if the rotation R is defined by
R = [ Twist ] [ pi/2  Dec ] [ pi/2 + RA ] 3 1 3where RA, Dec, and Twist are timedependent, then if we make the abbreviations
A(Twist) = [ Twist ] 3 B(Dec) = [ pi/2  Dec ] 1 C(RA) = [ pi/2 + RA ] 3we can write
d(R) d(A) d(Twist)  =  *  * B * C dt d(Twist) dt d(B)  d(Dec) + A *  *  * C d(Dec) dt d(C) d(RA) + A * B *  *  d(RA) dtThe derivatives of A, B, and C can be found using DROTAT. State transformations
P (t) Ibe a position vector referenced to an inertial frame ``I,'' and let
P (t) Nbe the equivalent position vector referenced to a noninertial frame ``N.'' If R(t) is the transformation from frame I to N at time t, then the two vectors are related as follows:
P (t) = R(t) P (t) N ITherefore, the derivatives of the position vectors satisfy
d [ P (t) ] d [ P (t) ] N I d [ R(t) ]  = R(t) *  +  * P (t) dt dt dt IIt's well to note that although R(t) may vary slowly, the second term in the above equation is not necessarily insignificant. For example, if R(t) describes a transformation between an inertial frame and a bodycentered frame that uses a bodycentertoSun vector to define one of its coordinate axes, then for any point that is fixed on this axis, the two addends above have equal and opposite magnitude. In particular, if the fixed point is the location of the Sun, the magnitude of the second addend is (ignoring the velocity of the Sun with respect to the inertial frame) that of the inertially referenced velocity of the body used to define the bodycentered frame. SPICELIB provides routines to transform states between inertial frames and bodyfixed planetocentric frames. The routine TISBOD returns the 6x6 transformation matrix required to transform inertially referenced state vectors to bodyfixed planetocentric coordinates. If REF is the name of the inertial frame of interest, BODY is the NAIF integer code of a body defining a bodyfixed planetocentric frame, and ET is ephemeris time used to define the bodyfixed frame, then the call
CALL TISBOD ( REF, BODY, ET, TSIPM )returns TSIPM, the desired 6x6 state transformation matrix. A state vector S can be transformed to the bodyfixed state vector SBFIXD by the routine call
CALL MXVG ( TSIPM, S, 6, 6, SBFIXD )Since the inverse of a state transformation matrix is not simply its transpose, SPICELIB provides the utility routine INVSTM to perform the inversion. If M is a state transformation matrix, the inverse matrix MINV can be obtained via the routine call
CALL INVSTM ( M, MINV ) Validating a rotation matrix
C C Set values for the column norm and determinant tolerances C NTOL and DTOL: C NTOL = 1.D7 DTOL = 1.D7 IF ( .NOT. ISROT ( M, NTOL, DTOL ) ) THEN [perform error handling] ELSE CALL M2Q ( M, Q ) . . . END IF Tutorial introduction to rotations
In this section, we make some assertions that we don't prove. Our goal is to supply you with the most important information first, and fill in the details later. Proofs are supplied only when they're instructive and not too distracting. The longer or more difficult proofs are deferred to the ``Mathematical road map'' chapter. A comment of the heuristic variety
We offer the following model: Take a soccer ball, put two fingers on diametrically opposed points on the ball, and rotate the ball through some angle, keeping your fingers in place. What you use to rotate the ball is up to you. Well, that's it. That's the effect of a rotation on a soccer ball. Now you're equipped to answer some questions about rotations. Do rotations preserve inner products of vectors? That is, is it true that for vectors u and v, and a rotation R,
< R u, R v > = < u, v > ?Well, presume that your soccer ball is centered at the origin, and mark the ball where u and v, or extensions of them, intercept the surface (perhaps you could hold a marker pen between your teeth). Does rotating the ball change the angular separation of the marks? No. So rotations preserve angular separation. Does rotating the ball change the norm of u or of v? No. So rotations preserve both angular separation and norms, and hence inner products. Do rotations preserve cross products? For vectors u and v, is it true that
( R u ) x ( R v ) = R ( u x v )?Mark the intercepts of u, v, and u x v on the soccer ball. After you rotate the ball, does the intercept mark of u x v still lie at the right place, relative to the u and v intercept marks? Yes. Since we already know that rotations preserve norms, we can conclude that they preserve cross products as well. The soccer ball model shows that rotations preserve geometrical relationships between vectors. Definition of ``rotation''
Definition 1
 R(v)  = v
R ( a x b ) = ( R a ) x ( R b ). Definition 2
R(n) = n.
Definition 3
Uses of the definitions
Definition (2) obviously implies definition (1). Less obviously, definition (1) implies definition (2). This had better be true if definition (1) is valid, since we expect rotations to have rotation axes. In the ``Mathematical road map'' chapter, we prove that the two definitions are equivalent. Definition (3) is a mathematical paraphrase of our soccer ball model of rotations. Definition of ``rotation'' and ``orthogonal'' matrix
Any matrix whose columns form an orthonormal set is called an ``orthogonal'' matrix. Rotations preserve inner products
Our definition of ``rotation'' says that rotations preserve norms of vectors. That is, if R is a rotation and v is a vector, then
 R(v)  =  v .Preserving norms also implies the seemingly stronger property of preserving inner products: if R preserves norms and u, v are vectors, then
2 2  R ( u  v )  =  u  v  = < u  v, u  v > 2 2 =  u   2 < u, v > +  v ,and also
2  R ( u  v )  = < R ( u  v ), R ( u  v ) > = < R(u)  R(v), R(u)  R(v) > 2 2 =  R(u)  +  R(v)   2 < R u, R v > 2 2 =  u  +  v   2 < R u, R v >so
< R(u), R(v) > = < u, v >.So rotations really do preserve inner products. In particular, for any orthonormal basis, the images of the basis vectors under a rotation are also an orthonormal set. Then rotation matrices, expressed relative to an orthonormal basis, are in fact orthogonal, as claimed. Inverses of rotation matrices
T R Ris the inner product of the ith column of R and the jth column of R, all entries of the product are zero except for those on the main diagonal, and the entries on the main diagonal are all 1. So
T R R = I.If A and B are square matrices with real or complex entries, it's a fact that if
A B = Ithen
B A = I.We won't prove this (the ``Nullity and Rank'' theorem is useful, if you wish to do so). But this result implies that
T R R = I,given the previous result. This shows that the rows of R form an orthonormal set as well, and that the transpose of R is also a rotation matrix. Composition of rotations
What does our soccer ball model say? We can rotate the ball as many times as we like, changing the rotation axis each time, without changing the distance of any surface point from the center, so norms are preserved. Similarly, after marking the intercepts of u, v, and u x v on the surface, we can perform any sequence of rotations without changing the position of the u x v intercept mark relative to those of u and v. So cross products are preserved. That's all we need to verify that the composition of a sequence of rotations is a rotation. It follows that the product of a sequence of rotation matrices is a rotation matrix. If you don't agree that that's all we need, we can present the same argument using the usual symbols: Suppose R1 and R2 are rotation mappings, and for any vector v,
R3(v) = R2 ( R1(v) ).Then for any vector v, we have
 R3(v)  =  R2 ( R1 v )  =  R1( v )  =  v .So R3 preserves norms. If u and v are vectors, then
R3 ( u x v ) = R2 ( R1 ( u x v ) ) = R2 ( R1(u) x R1(v) ) = ( R2 ( R1(u) ) x R2 ( R1(v) ) ) = R3(u) x R3(v),so R3 preserves cross products. We've used only the definition of R3 and the fact that R2 and R1 preserve cross products in this proof. We conclude that R3 is a rotation. We can extend the result to the product of a finite number of rotations by mathematical induction; the argument we've made is almost identical to the induction step. Coordinate transformations
How do we transform a vector v from one coordinate system to another? This result really belongs to linear algebra, but we'll state it here because it seems to come up a lot. Given two vector space bases,
B1 = { e1, e2, e3 }, B2 = { u1, u2, u3 },and a vector v that has components ( v1, v2, v3 ) relative to B1, we wish to express v relative to B2. We can say that
v = x1 u1 + x2 u2 + x3 u3,where the x's are unknowns. Let M be the matrix whose columns are u1, u2, and u3, represented relative to basis B1. M represents the linear transformation T defined by
T (e1) = u1, T(e2) = u2, T(e3) = u3.Then since
1 1 T (v) = T ( x1 u1 + x2 u2 + x3 u3 ) = ( x1 e1 + x2 e2 + x3 e3 ),we see that
1 M v = ( x1, x2, x3 ).So we've found the components of v, relative to basis B2. In the case where B1 and B2 are orthonormal bases, the matrix M is orthogonal. So we have
T M v = ( x1, x2, x3 ).Conversely, if M is the matrix that transforms vectors from orthonormal basis B1 to orthonormal basis B2, then the rows of M are the basis vectors of B2. For example, if M is the matrix that transforms vectors from J2000 coordinates to body equator and prime meridian coordinates, then the first row is the vector, expressed in J2000 coordinates, that points from the body center to the intersection of the prime meridian and body equator. The third row is the vector, expressed in J2000 coordinates, that points from the body center toward the body's north pole. Rotation of vectors in the plane
We can assume that v is a unit vector; since rotations are linear, it's easy to extend the result to vectors of any length. Now, if v is (1,0), the result of the rotation will be
( cos(theta), sin(theta) ).How does this help us if v is an arbitrary unit vector? Given a unit vector v, let v' be the vector perpendicular to v, obtained by rotating v by pi/2. Now v and v' form an orthonormal basis, and relative to this basis, v has coordinates (1,0). But we've already found out what we get by rotating v by theta radians: relative to our new basis, the result must be
( cos(theta), sin(theta) ).Relative to our original basis, this vector is
cos(theta) v + sin(theta) v'This is the result we're looking for: ``If you rotate a vector v by theta radians, you end up with cos(theta) v plus sin(theta) v','' where v' is v, rotated by pi/2. Scaling v does not affect this result. A consequence of this result is that the mapping R that rotates vectors by theta radians is represented by the matrix
+ +  cos(theta) sin(theta)   .  sin(theta) cos(theta)  + +It is useful to note that R has this exact representation relative to any orthonormal basis where the second vector is obtained from the first by a rotation of pi/2. A canonical representation for rotations
The twodimensional diagram below shows this decomposition. All of the vectors lie in the plane containing r and n.
\ . rParallel \ . \ . r \ . \ . \ . \ . \ ^ \  \  n \  ..........\ rPerpNow, what does R do to vectors that are perpendicular to n? Since R rotates each vector about n, if a vector v is perpendicular to n, then R(v) is perpendicular to n as well (remember that rotations preserve inner products, and orthogonality in particular). So rPerp just gets rotated in the plane perpendicular to n. We know from the last section how to find R(rPerp): if we let rPerp' be the vector obtained by rotating rPerp by pi/2 about n, then
R(rPerp) = cos(theta) rPerp + sin(theta) rPerp'We will also want to know what R(rPerp') is. Since rotating rPerp' by pi/2 about n yields rPerp, applying our familiar formula to rPerp' gives us
R(rPerp') = cos(theta) rPerp'  sin(theta) rPerp.Now, since n, rPerp, and rPerp' are mutually orthogonal, these vectors form a basis. Since we can scale r so that rPerp has norm 1, and since rPerp' has the same norm as rPerp, we may assume that the basis is actually orthonormal. The matrix of R relative to this basis is
+ +  1 0 0     0 cos(theta) sin(theta) .    0 sin(theta) cos(theta)  + +Since the rotation we're representing is arbitrary, we've shown that every rotation can be represented by a matrix of the above form. Equivalently, every rotation matrix is similar to one of the above form. This fact justifies the use of the term ``canonical form.'' The canonical form we've found shows why threedimensional rotations are very much like twodimensional rotations: The effect of a threedimensional rotation on any vector is to rotate the component of that vector that is normal to the rotation axis, and leave the component parallel to the rotation axis fixed. This rotation matrix is a useful ``model'' to keep in mind when dealing with rotations because of its particularly simple form. It's easy to read off some types of information directly from this matrix. Some examples:
Rotation axis and angle
Given a rotation R and a vector v, normal to the rotation axis n of R, the angle between v and R(v), measured counterclockwise around n, is the rotation angle of R. We see that the rotation angle depends on the direction of the axis: if we pick n as the axis, we change the sign of the angle. Note that while the rotation axis and angle of a rotation are not uniquely defined, a choice of axis and angle do determine a unique rotation. How do we find the rotation matrix R that rotates vectors by angle theta about the unit vector n? If n is
n = (n1, n2, n3),then
2 R = I + ( 1  cos(theta) ) N + sin(theta) N,where
+ +  0 n3 n2    N =  n3 0 n1 .    n2 n1 0  + +How do we recover the rotation angle and axis of a rotation R from a corresponding rotation matrix, M? We've already seen in the ``canonical form'' section that the rotation angle is
ACOS ( ( Trace(M)  1 ) / 2 ).If the rotation angle is not zero or pi, then the relation
T M  M = 2 sin(theta) Nallows us to recover the rotation axis n from M, while if the rotation angle is pi, we have
2 M = I + 2 N,again determining n. In the ``Mathematical road map'' chapter, we'll verify these assertions. Timedependent coordinate transformations
B2 = { v1(t), v2(t), v3(t) }.An example of a timedependent coordinate transformation is the transformation from J2000 to body equator and prime meridian coordinates. If R(t) transforms vectors from basis B1 to basis B2, the basis vectors of B2 are the rows of the matrix R(t). Let p(t) and p'(t) be position and velocity vectors expressed relative to B1. What is the corresponding velocity, expressed relative to B2? We know that p(t) has coordinates
R(t) p(t)relative to B2, so the time derivative of R(p(t)) is
R(t) p'(t) + R'(t) p(t),relative to B2. If R(t) is expressed as a product of the form
R(t) = [ w1(t) ] [ w2(t) ] [ w3(t) ] , i j kthen
R'(t) = [ w1(t) ]' [ w2(t) ] [ w3(t) ] i j k + [ w1(t) ] [ w2(t) ]' [ w3(t) ] i j k + [ w1(t) ] [ w2(t) ] [ w3(t) ]' i j kSince we know the explicit form of the factors (given in the ``Notation'' section), we can compute R'(t). We must take care when converting velocity vectors between systems whose bases are related in a timedependent way. If R(t) varies extremely slowly, it is tempting to ignore the R' term, and in fact this is a valid approximation in some cases. However, since the magnitude of this term is proportional to the magnitude of p, the term can be large when R is quite slowly varying. An example: Let B1 be the basis vectors of the J2000 system, and let
B2 = { v1(t), v2(t), v3(t) }be defined as follows: v1(t) is the geometric JupiterSun vector at ephemeris time t, v3(t) is orthogonal to v1(t) and lies in the plane containing v1(t) and Jupiter's pole at time t, and v2(t) is the cross product of v3(t) and v1(t). Let R(t) be the transformation matrix from basis B1 to B2. Then the period of R(t) is 1 Jovian year (we're ignoring movement of Jupiter's pole). Now if p(t) is the JupiterSun vector in J2000 coordinates, then p'(t) is the negative of Jupiter's velocity in J2000 coordinates. But in B2 coordinates, R(t) ( p(t) ) always lies along the xaxis, and if we approximate Jupiter's motion as a circle, then R(p(t))' is the zero vector. So we have the equation
R(t) p'(t) + R'(t) p(t) = [ R(t)( p(t) ) ]' = 0,which implies
R'(t)p(t) =  R(t) p'(t).So in this case, the term involving R' has the same magnitude as the term involving R, even though R is slowly varying. Euler angles
M = [w1] [w2] [w3] . i1 i2 i3The angles w1, w2, and w3 are called ``Euler angles.'' It is not necessarily obvious that this ``factorization'' is possible. It turns out that as long as i2 does not equal i1 or i3, it is possible, for any rotation matrix M. In the ``Mathematical road map'' chapter, we exhibit the formulas for calculating w1, w2, and w3, given M and i1, i2, and i3. Quaternions
Unit quaternions may be associated with rotations in the following way: if a rotation R has unit vector n = (n1, n2, n3) as an axis and w as a rotation angle, then we represent R by
Q = ( cos(w/2), sin(w/2) n1, sin(w/2) n2, sin(w/2) n3 ).As you might suspect, this association is not unique: substituting (w + 2*pi) for w, we see that Q is also a representation for R. If we choose the rotation axis and angle of R so that the angle lies in [0, pi], then there is a unique quaternion representing R, except in the case where R is a rotation by pi radians. Quaternion arithmetic
Definitions
The main interest of quaternion multiplication is that we can actually carry out composition of rotations using the multiplication defined on the quaternions. If quaternions Q1 and Q2 represent rotations R1 and R2, then Q2*Q1 represents R2(R1). So the mapping from unit quaternions to rotations is a group homomorphism, where the ``multiplication'' operation on the rotations is functional composition. Quaternion addition is simple vector addition. Multiplication is a little more complicated. Before defining it, we're going to introduce a new notation for quaternions that makes it easier to deal with products. The quaternion
Q = ( Q0, Q1, Q2, Q3 )can be represented as
Q0 + ( Q1, Q2, Q3 ),or
s + v,where s represents the ``scalar'' Q0 and v represents the ``vector''
( Q1, Q2, Q3 ).We define the ``conjugate'' of the quaternion
q = s + vas
* q = s  v.Given two quaternions,
q1 = s1 + v1, q2 = s2 + v2,we define the product q1 * q2 as
( s1 * s2  < v1, v2 > ) + ( s1 * v2 + s2 * v1 + v1 x v2 ).We've grouped the ``scalar'' and ``vector'' portions of the product. Basic properties of multiplication
Is multiplication commutative? No; if s1 and s2 above are zero, then the product is
 < v1, v2 > + v1 x v2,which is not commutative. However, multiplication is associative: given three quaternions q1, q2, and q3, we have
q3 * ( q2 * q1 ) = ( q3 * q2 ) * q1.We'll forgo checking this; it's messy but straightforward. If you do check it, the vector identities
A x ( B x C ) = < A, C > B  < A, B > C ( A x B ) x C = < C, A > B  < C, B > Awill be useful. What's the product of q and its conjugate? It comes out to
2  q .What's the conjugate of the product
q1 * q2,where q1 and q2 are as defined above? The product formula allows us to verify that the answer is
* * q2 * q1. Deducing the multiplication formula
You can check this: if we scale q1 by x, the product gets scaled by x. The same thing happens if we scale q2. If we replace q1 by the sum of two quaternions, say
q1 = q + q' = ( s + s' ) + ( v + v' ),the product is
q * q2 + q' * q2.The analogous result occurs when we replace q2 by a sum of two quaternions. Because of this linearity property, we can define multiplication on a small set of quaternions, and then define multiplication on the whole set of quaternions by insisting that the multiplication operator is linear in both operands. This gives us an equivalent definition of multiplication. To carry out this definition, we first define multiplication on the four quaternions
1 + ( 0, 0, 0 ), which we call ``1,'' 0 + ( 1, 0, 0 ), which we call ``i,'' 0 + ( 0, 1, 0 ), which we call ``j,'' 0 + ( 0, 0, 1 ), which we call ``k.''We treat ``1'' as a scalar and i, j, and k as vectors, and define the products
1 * v = v, for v = i, j, k; v * v = 1, for v = i, j, k; v1 * v2 =  v2 * v1, for v1, v2 = i, j, k; i * j = k; j * k = i; k * i = j.Multiplication of i, j, and k works just like taking cross products. If we now proclaim that multiplication is linear in both operands, then since all quaternions can be expressed as linear combinations of ``1,'' i, j, and k, we've defined multiplication on the entire set of quaternions. You can check that this definition of multiplication is consistent with our formula above. Composing rotations using quaternions
We've defined a mapping from quaternions to rotations, since the relation
Q = ( cos(w/2), sin(w/2) n1, sin(w/2) n2, sin(w/2) n3 )allows us to recover w and the axis ( n1, n2, n3 ), hence the corresponding rotation. Now suppose we have two quaternions Q1 and Q2 that represent rotations R1 and R2, respectively. We're claiming that the product Q2 * Q1 represents R2(R1). So, we should be able to recover the rotation axis and angle of R2(R1) from the quaternion Q2 * Q1. In the ``Mathematical road map' chapter, we will verify this claim. Mathematical road map
The difference between the two perspectives is a bit like the difference between having a set of directions to get from point A to point B, and having a road map of the entire area. This chapter is not organized sequentially, since there is little logical dependence of one section on another. It is simply a collection of discussions. Rotation of a vector about an axis
As in the tutorial discussion of the canonical form for rotations, we can express r as the sum of two orthogonal components:
r = rParallel + rPerp.Let's give the name rPerp' to the vector obtained by rotating rPerp by pi/2 radians about n. We know, from the results of the ``canonical form'' section, that applying our rotation to r will yield
rParallel + cos(theta) rPerp + sin(theta) rPerp'So all we have to do is find rPerp and rPerp' in terms of r, n, and theta. It turns out that rPerp' is precisely n x r, since n x r is parallel to rPerp' and has the same magnitude as rPerp, namely
r sin(phi),where phi is the angle between r and n. Rotating rPerp' by another pi/2 radians yields rPerp, so
rPerp = n x ( n x r ).In the picture below,
n x ( n x r)
 n x ( n x r)
\ . rParallel \ phi . \ / . r \ . \ . \ . \ . \ ^ \  \  n \ _  n x ( n x r ) ..........._......... n x ( n x r ) /_/ = rPerp . =  rPerp . . . n x r = rPerp'Now we're ready to compute the image of r under the rotation. It is:
rParallel + cos(theta) rPerp + sin(theta) rPerp' = ( r  rPerp ) + cos(theta) rPerp + sin(theta) rPerp' = r + ( cos(theta)  1 ) rPerp + sin(theta) rPerp' = r + ( 1  cos(theta) ) ( n x ( n x r ) ) + sin(theta) ( n x r ).This is what we were after: an expression for the image of r under the rotation, given in terms of r, theta, and n. Formation of a rotation matrix from axis and angle
What's the rotation matrix R that rotates vectors by theta radians about the vector n? If n is a unit vector, then the result of the last section implies that
R * r = r + ( 1  cos(theta) ) ( n x ( n x r ) ) + sin(theta) ( n x r ),for any vector r. Now, let
n = (n1, n2, n3),and define the matrix N by
+ +  0 n3 n2    N =  n3 0 n1 ;    n2 n1 0  + +this definition implies that
N * r = n x rfor all r. So we can rewrite the above expression as
R * r = r + ( 1  cos(theta) ) ( N * ( N * r ) ) + sin(theta) ( N * r ),or
2 R * r = [ I + ( 1  cos(theta) ) N + sin(theta) N ] r.Since r is arbitrary, we must have
2 R = I + ( 1  cos(theta) ) N + sin(theta) N.R is the desired matrix. Finding the axis and angle of a rotation matrix
There are many ways to recover the rotation axis. The most elegant method we know of is presented in [1]. The idea is based on the observation that any rotation matrix R can be expressed by
2 R = I + ( 1  cos(theta) ) N + sin(theta) N,where N is derived from the rotation axis, as in the last section. Now N is skewsymmetric and N squared is symmetric, so
T R  R = 2 sin(theta) N.As long as sin(theta) is nonzero, we've found N and hence the axis itself. If theta is pi, we have
2 R = I + 2 N,which still allows us to recover the axis. In the tutorial section, we showed that the rotation angle can be recovered from the trace of a rotation matrix:
angle = ACOS ( ( trace  1 ) / 2 ).If the angle is very small, we will determine it more accurately from the relation
T R  R = 2 sin(theta) N. Formation of a rotation matrix from a quaternion
There is a fast, accurate solution available. It depends on the formula relating a rotation matrix to its axis and angle, which we derived earlier in the chapter. In this approach, we compute the matrix corresponding to a quaternion, component by component. Define
c = cos(theta/2), s = sin(theta/2),and let the quaternion
q = c + s n = q0 + s ( q1, q2, q3 )represent a rotation R having unit axis vector n and rotation angle theta. If n = ( n1, n2, n3 ), and we define the matrix N by
+ +  0 n3 n2    N =  n3 0 n1 ,    n2 n1 0  + +then the matrix M representing R is
2 M = I + ( 1  cos(theta) ) N + sin(theta) N.Now we can make the substitutions
sin(theta) = 2 c s, 2 ( 1  cos(theta) ) = 2 sto obtain
2 M = I + 2 ( s N ) + 2 c ( s N ).Substituting the elements of our quaternion into s N, we find
+ +  2 2  + +  (q2 + q3 ) q1 q2 q1 q3   1 0 0       2 2  M =  0 1 0  + 2  q1 q2 (q1 + q3 ) q2 q3       0 0 1   2 2  + +  q1 q3 q2 q3 (q1 + q2 )  + + + +    0 q0 q3 q0 q2      + 2  q0 q3 0 q0 q1 ,      q0 q2 q0 q1 0  + +so
+ +  2 2   1  2 ( q2 + q3 ) 2 (q1 q2  q0 q3) 2 (q1 q3 + q0 q2)       2 2  M =  2 (q1 q2 + q0 q3) 1  2 ( q1 + q3 ) 2 (q2 q3  q0 q1) .      2 2   2 (q1 q3  q0 q2) 2 (q2 q3 + q0 q1) 1  2 ( q1 + q2 )    + + Equivalence of rotation definitions
We wish to prove that definitions (1) and (2) from the ``Definition of rotations'' section of the tutorial are equivalent. To do this, we need to show that a mapping R that satisfies definition (1) also satisfies definition (2). This amounts to showing that R has a fixed subspace of dimension 1, or equivalently, that R has 1 as one of its eigenvalues. An algebraic approach
A geometric approach
We're going to look at the effect of R on the unit sphere, and demonstrate that two points on the sphere are fixed. We'll assume that the rotation is not the identity and does not map any vector v to v. This last case corresponds to a rotation of pi radians. Our first observation is that R maps great circles to great circles. This follows from the fact that a great circle is a set of unit vectors, all orthogonal to some particular vector v. Since R preserves inner products, the image of the great circle is a set of unit vectors, all orthogonal to R(v). Now, consider the distances that vectors on the unit sphere move when the rotation R is applied; there is some vector v, not necessarily unique, that moves the maximum distance. Let C1 be a great circle passing through v and R(v), and let C2 be a great circle that passes through v and intersects C1 at right angles. Now R(C2) passes through R(v), and if we can show that it passes through at right angles to C1, then C2 and R(C2) intersect at vectors p and p, both of which are normal to v and R(v). So R(p) is either p or p. But we've assumed that R does not map any vector to its inverse, so R(p) = p, and we have a fixed vector. So, we must show that R(C2) passes through R(v) at right angles to C1. If it did not, there would be some point w on C2, close to v, such that
 R(w)  w  >  R(v)  v ,contradicting our hypothesis that no vector moves farther that v. We will leave the rigorous proof of this last assertion to the energetic reader. Quaternion multiplication
There are two assertions that we need to prove:
* R(v) = q * v * q,
q2 * q1. Assertion 1
R(v) = v + sin(theta) n x v + ( 1  cos(theta) ) n x ( n x v ),where n is a unit axis vector and theta is the corresponding rotation angle. We also define the constants C and S by
C = cos(theta/2), S = sin(theta/2).The quaternion
q = C + S nrepresents R. To check the assertion, we compute
( C + S n ) * v * ( C  S n ) = ( C + S n ) * [ ( S <v,n> ) + ( C v  S v x n ) ] 2 = [ C S <v, n>  S C < n, v > + S < n, v x n > ] 2 2 + [ C v  C S v x n + S <v, n> n + S C n x v 2  S n x ( v x n ) ].Since n is normal to v x n, the scalar part of the last line is zero, which leaves us with
2 2 2 C v + S <v, n> n + 2 S C n x v  S n x ( v x n ).We can rewrite this again as
2 2 2 C v + S <v, n> n + 2 S C n x v + 2 S n x ( n x v ) 2  S n x ( n x v ),and using the vector identity
A x ( B x C ) = < A, C > B  < A, B > C,we can modify the final term to arrive at
2 2 2 C v + S <v, n> n + 2 S C n x v + 2 S n x ( n x v ) 2  S ( <n, v> n  <n, n> v ).Since n is a unit vector, the entire expression reduces to
2 v + 2 S C n x v + 2 S n x ( n x v ) = v + 2 sin(theta/2) cos(theta/2) n x v 2 + 2 sin (theta/2) n x ( n x v) = v + sin(theta) n x v + ( 1  cos(theta) ) n x ( n x v ) = R(v). Assertion 2
* q1 * v * q1 ,we can express R2(R1(v)) by
* * q2 * ( q1 * v * q1 ) * q2 * = ( q2 * q1 ) * v * ( q2 * q1 ).Now let q be a quaternion that represents R2(R1); then
* q * v * q = R2(R1(v))for all v. We'll be done if we can show that, in general, for unit quaternions x and y, if
* * x * v * x = y * v * yfor all vectors v, then x equals y or y. But this equation implies that
* * y * x * v = v * y * x,for all v, which in turn implies that
* y * xis a scalar, since only scalar quaternions commute with every vector quaternion (due to the cross product term in the product formula). Since y and x are unit quaternions, either
* y * x = 1or
* y * x = 1,so x = y or y. We conclude that q2 * q1 = q or q, so q2 * q1 does represent R2(R1). Recovery of Euler angles from a rotation matrix
M = [w1] [w2] [w3] . (1) i1 i2 i3There are a couple of reasons why we might want to solve this problem: first, the representation of a rotation by three Euler angles is a common one, so it is convenient to be able to convert the matrix representation to this form. Also, the three angles on the right hand side of equation (1) often allow you to visualize a rotation more readily than does the matrix representation M. This ``factorization'' is possible if i2 does not equal i1 or i3. For each valid sequence (i1i2i3) of axes, there is a set of functions that give us w1, w2, and w3 as a function of M:
w1 = f1 ( M ), i1i2i3 w2 = f2 ( M ), i1i2i3 w3 = f3 ( M ). i1i2i3How can we derive the functions
f1 , f2 , f3 ? i1i2i3 i1i2i3 i1i2i3One approach is to multiply the matrices on the right hand side of equation (1); this yields a matrix whose entries are sums of products of sines and cosines of w1, w2, and w3. We can then equate the entries of this matrix to those of M, and find formulas for w1, w2, and w3 that arise from the componentwise correspondence. In subsequent sections, we actually carry out this procedure for 313 and 123 factorizations. There are twelve sets of axes to consider, so there are potentially twelve sets of functions to compute. However, the procedure we've just described is not enjoyable enough to justify doing it twelve times. We'd like to find a slicker way of solving the problem. One approach is to find a way of ``recycling'' the formulas we derived for one particular axis sequence. Here's an example of how we might do this: Suppose that we already have functions
f1 , f2 , f3 313 313 313that allow us to factor rotation matrix M as a 313 product: If
M = [w1] [w2] [w3] . (2) 3 1 3then
w1 = f1 ( M ), 313 w2 = f2 ( M ), 313 w3 = f3 ( M ). 313We'd like to somehow use the functions we've already got to factor M as a 232 product: we want to find functions
f1 , f2 , f3 232 232 232such that
M = [y1] [y2] [y3] , (3) 2 3 2and
y1 = f1 ( M ), 232 y2 = f2 ( M ), 232 y3 = f3 ( M ) 232without having to derive
f1 , f2 , f3 232 232 232from scratch. We'll start out by using a new basis, relative to which the right hand side of (3) is not a 232, but rather a 313 rotation. It is important to note here that bases are ordered sets of vectors; changing the order changes the basis. Let the basis B1 be the ordered set of vectors
{e(1), e(2), e(3)},and let the basis B2 be the ordered set of vectors
{e(3), e(1), e(2)}.Now the rotation matrix
[y] 2expressed relative to B1 represents the same rotation as the matrix
[y] 3expressed relative to B2. Both matrices represent a rotation of y radians about the vector e(2). Similarly, the matrix
M = [y1] [y2] [y3] 2 3 2expressed relative to B1 represents the same rotation as the matrix
M' = [y1] [y2] [y3] 3 1 3expressed relative to B2. So if C is the matrix whose columns are the elements of B2, expressed relative to B1, namely
+ +  0 1 0  C =  0 0 1 ,  1 0 0  + +then
1 C M C = M' (4)We can use the functions
f1 , f2 , f3 313 313 313to factor M' as a 313 product: applying (4), we have
1 y1 = f1 ( C M C ), (5) 1 y2 = f2 ( C M C ), (6) 1 y3 = f3 ( C M C ), (7)so we've found functions that yield the angles y1, y2 and y3 that we sought. ``No muss, no fuss.'' How much mileage can we get out of our 313 factorization functions? Looking at our example, we see that the main ``trick'' is to find a basis so that the factorization we want is a 313 factorization with respect to that basis. It is important the new basis be righthanded; otherwise the form of the matrices
[w] iis not preserved. It turns out that for any axis sequence of the form aba, we can find a righthanded basis such that the factorization we want is a 313 factorization with respect to that basis. There are two cases: if we define a successor function s on the integers 1, 2, 3 such that
s(1) = 2, s(2) = 3, s(3) = 1,we either have b = s(a) or a = s(b). In the first case, b = s(a), and if our original ordered basis is
B1 = { e(1), e(2), e(3) },then
B2 = { e(b), e( s(b) ), e(a) }is the righthanded basis we're looking for. You can check that
e(a) = e(b) x e( s(b) ).We recall that the transformation matrix C we require has the elements of B2 as columns. For example, if a is 2 and b is 3, then B2 is
{ e(3), e(1), e(2) },and the matrix C is
+ +  0 1 0  C =  0 0 1 ,  1 0 0  + +as we have seen previously. The axis sequences that can be handled by the above procedure are 121, 232, and 313. In the second case, a = s(b), and if our original ordered basis is
B1 = { e(1), e(2), e(3) },then
B2 = { e(b), e( s(a) ), e(a) }is the righthanded basis we're looking for. Again, you can verify this by taking cross products. The transformation matrix C we require has the elements of B2 as columns. For example, if a = 2 and b = 1, then B2 is
{ e(1), e(3), e(2)) }and the matrix C is
+ +  1 0 0  C =  0 0 1 .  0 1 0  + +The axis sequences that can be handled by the above procedure are 131, 212, and 323. So we can use our 313 formula to handle all of the aba factorizations, just by computing the correct transformation matrix C. What about abc factorizations? As you might guess, the procedure we've described also applies to these, with very little modification. Suppose we have the formulas we need to carry out a 123 factorization. We'd like to find a basis that allows us to represent the abc product
[w1] [w2] [w3] a b cas the product
[w1] [w2] [w3] . 1 2 3Again, there are two cases, depending on whether b is the successor of a or a is the successor of b, according to our cyclic ordering. In the case where b is the successor of a, the righthanded basis we want is
B2 = { e(a), e(b), e(c) }.With respect to the basis B2, our abc factorization is a 123 factorization. Again, we can form the transformation matrix C by letting its columns be the elements of B2. In the second case, a is the successor of b. Our new basis is
B2 = { e(a), e(b), e(c)}.In this case, there is a slight twist: the change of basis we use negates the third rotation angle. This is not a serious problem; the change of basis converts the product
[w1] [w2] [w3] a b cto
[w1] [w2] [w3] , 1 2 3so we can still recover the angles w1, w2, and w3 easily. So our 123 factorization formula allows us to handle all the abc factorizations. Having shown that we can perform all of the aba and abc factorizations using just one formula for each type of factorization, we now proceed to derive those formulas. This is not a particularly instructive procedure, but the derivations ought to be written down somewhere, and this is as good a place as any. Euler angle recovery: aba case
In this case, the right hand side of (1) is
+ + + +  cos(w1) sin(w1) 0   1 0 0   sin(w1) cos(w1) 0  *  0 cos(w2) sin(w2)  *  0 0 1   0 sin(w2) cos(w2)  + + + + + +  cos(w3) sin(w3) 0   sin(w3) cos(w3) 0   0 0 1  + +which equals
+ +  cos(w1) sin(w1) 0   sin(w1) cos(w1) 0  *  0 0 1  + + + +  cos(w3) sin(w3) 0   cos(w2)sin(w3) cos(w2)cos(w3) sin(w2)   sin(w2)sin(w3) sin(w2)cos(w3) cos(w2)  + +which comes out to
+ +  cos(w1)cos(w3) cos(w1)sin(w3) sin(w1)sin(w2)   sin(w1)cos(w2)sin(w3) +sin(w1)cos(w2)cos(w3)     sin(w1)cos(w3) sin(w1)sin(w3) cos(w1)sin(w2) .  cos(w1)cos(w2)sin(w3) +cos(w1)cos(w2)cos(w3)     sin(w2)sin(w3) sin(w2)cos(w3) cos(w2)  + +At this point, we can recover w1, w2, and w3 from the elements of M. The inverse trigonometric functions used below are borrowed from Fortran. We find w2 from the relation
w2 = ACOS( M(3,3) ).(So w2 is in [0, pi].) If w2 is not equal to 0 or pi, then we can recover w1 as follows:
M(1,3) sin(w1)sin(w2)  =  = tan(w1), M(2,3) cos(w1)sin(w2) w1 = ATAN2 ( M(1,3), M(2,3) ).We find w3 in an analogous fashion, again assuming w2 is not equal to 0 or pi. We find
M(3,1) sin(w2)sin(w3)  =  = tan(w3), M(3,2) sin(w2)cos(w3) w3 = ATAN2 ( M(3,1), M(3,2) ).Note the minus sign used in the second ATAN2 argument. For ATAN2 to determine the correct value, it is necessary that the first and second arguments have the same signs as sin(w3) and cos(w3), respectively. Now if w2 is equal to zero or pi, we have a degenerate case: M is the product of two rotations about the third coordinate axis. The angles of the rotations are not determined uniquely, only the sum of the angles is. One way of finding a factorization is to set w3 to zero, and solve for w1. The matrix M then is equal to
+ +  cos(w1) cos(w2)sin(w1) 0   sin(w1) cos(w2)cos(w1) 0 ,  0 0 cos(w2)  + +so we can recover w1 by computing
w1 = ATAN2( M(2,1), M(1,1) ). Euler angle recovery: abc case
In this case, the right hand side of (1) is
+ + + +  1 0 0   cos(w2) 0 sin(w2)   0 cos(w1) sin(w1)  *  0 1 0  *  0 sin(w1) cos(w1)   sin(w2) 0 cos(w2)  + + + + + +  cos(w3) sin(w3) 0   sin(w3) cos(w3) 0   0 0 1  + +which equals
+ +  1 0 0   0 cos(w1) sin(w1)  *  0 sin(w1) cos(w1)  + + + +  cos(w2)cos(w3) cos(w2)sin(w3) sin(w2)     sin(w3) cos(w3) 0     sin(w2)cos(w3) sin(w2)sin(w3) cos(w2)  + +which comes out to
+ +  cos(w2)cos(w3) cos(w2)sin(w3) sin(w2)     cos(w1)sin(w3) cos(w1)cos(w3) sin(w1)cos(w2)   +sin(w1)sin(w2)cos(w3) +sin(w1)sin(w2)sin(w3) .    sin(w1)sin(w3) sin(w1)cos(w3) cos(w1)cos(w2)   +cos(w1)sin(w2)cos(w3) +cos(w1)sin(w2)sin(w3)  + +We recover w2 by
w2 = ASIN ( M(1,3) ),so w2 is in the interval [pi/2, pi/2]. As long as w2 does not equal pi/2 or pi/2, we can find w1 by the formula
w1 = ATAN2 ( M(2,3), M(3,3) ),and w3 from the formula
w3 = ATAN2 ( M(1,2), M(1,1) ).If w2 is pi/2 or pi/2, we have a degenerate case. The sum of w1 and w3 is determined, but w1 and w3 are not determined individually. We can set w3 to zero, which reduces our right hand side to
+ +  0 0 sin(w2)   sin(w1)sin(w2) cos(w1) 0 ,  cos(w1)sin(w2) sin(w1) 0  + +so we can recover w1 from the formula
w1 = ATAN2 ( M(3,2), M(2,2) ). Appendix A: Document Revision HistoryMarch 9, 2017
May 27, 2010
November 17, 2005
January 10, 2005
February 2, 2004
December 2, 2002
April 26, 1999
