| Rotation |
Table of ContentsRotation Abstract Introduction Using this document References Notation CSPICE Functions Categories of functions Euler angle functions 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_c and m2q_c 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_c and raxisa_c are approximate inverses Using raxisa_c and axisar_c 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 Time-dependent 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: a-b-a case Euler angle recovery: a-b-c 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 error-prone 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 Three-dimensional Euclidean space.
E2 Two-dimensional 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 c
is sometimes referred to as a ``a-b-c'' rotation. For example, we may
talk about a ``3-1-3'' 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. CSPICE Functions
The CSPICE functions that deal with general rotations are:
Categories of functions
Euler angle functions
Constructing a matrix from Euler angles
The call
eul2m_c ( ang[2], ang[1], ang[0], i[2], i[1], i[0], m );returns the matrix `m', where
m = [ang[2]] [ang[1]] [ang[0]] .
i[2] i[1] i[0]
The indices i[2], i[1], i[0] must belong to the set
{1, 2, 3}.
Finding Euler angles that represent a matrix
m = [ang[2]] [ang[1]] [ang[0]] .
i[2] i[1] i[0]
The call
m2eul_c ( m, i[2], i[1], i[0], ang+2, ang+1, ang );returns the desired angles. m2eul_c restricts the ranges of the output angles so as to guarantee that the Euler angle representation is unique. The output angles ang[2] and ang[0] are always in the range (-pi,pi]. The range of ang[1] is determined by the set of rotation axes. When i[2] equals i[0], ang[1] is in the range [0, pi]. Otherwise, ang[1] 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_c sets the first angle (called ang[2] above) to zero. The other two angles are then uniquely determined. Again, the indices i[2], i[1], i[0] are members of
{1, 2, 3}.
There is a restriction on the allowed set of coordinate axes: i[1] must
not equal i[2] or i[0]. If this constraint is not met, the desired
representation of `m' by Euler angles may not exist; m2eul_c 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 2-1-3 rotation). In practical terms, if eul2m_c encounters one of these special matrices, eul2m_c must choose the Euler angles. Immediately the possibility arises that eul2m_c 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 round-off error. Euler angles near the limits of their ranges should be regarded with suspicion. The existence of singularities in the matrix-to-Euler angle mapping prevents eul2m_c and m2eul_c from being exact inverses: most of the time, the code fragment
eul2m_c ( ang[2], ang[1], ang[0], axis[2], axis[1], axis[0], m ); m2eul_c ( m, axis[2], axis[1], axis[0], ang+2, ang+1, ang );leaves the angles ang[2], ang[1], ang[0] unchanged, except for round-off error, but in some cases, the angles may change drastically. If we reverse the order of the function 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 3
so the code fragments
eul2m_c ( twist, halfpi_c() - dec, halfpi_c() + ra, 3, 1, 3, c );and
m2eul_c ( c, 3, 1, 3, ang+2, ang+1, ang ); twist = ang[2]; dec = halfpi_c() - ang[1]; ra = ang[0] - halfpi_c();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
m2q_c ( 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_c are
q[0] = cos( angle/2 ) q[1] = sin( angle/2 ) * a[0] q[2] = sin( angle/2 ) * a[1] q[3] = sin( angle/2 ) * a[2]m2q_c 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
q2m_c ( q, m ); q2m_c and m2q_c are approximate inverses of each other
m2q_c ( m, q ); q2m_c ( q, m );always preserve `m', except for round-off error. However, since there are two quaternions that represent each rotation, the sequence of calls
q2m_c ( q, m ); m2q_c ( m, q );do not necessarily preserve `q'. Multiplying quaternions
qxq_c ( 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
qdq2av_c ( 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]
i
to a vector, use rotvec_c. The following code fragment applies
[ angle ]
3
to the vector `v', yielding `vout'.
rotvec_c ( 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]
i
to a matrix using rotmat_c, as follows:
rotmat_c ( m, angle, 3, mout );After this function call, `mout' is equal to
[ angle ] * m
3
Rotating a vector about an arbitrary axis
vrotv_c ( 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
axisar_c ( 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_c for this. Let `axis' be the coordinate system rotation axis and `angle' be the rotation angle; then the code fragment
axisar_c ( 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_c is designed this way for compatibility with raxisa_c, which is an inverse routine for axisar_c. Finding the axis and angle of a rotation matrix
raxisa_c ( 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_c and raxisa_c are approximate inverses
raxisa_c ( m, axis, &angle ); axisar_c ( axist, angle, m );leaves `m' unchanged, except for round-off error. The code fragment
axisar_c ( axist, angle, m ); raxisa_c ( m, axis, &angle );usually leaves `axis' and `angle' unchanged, except for round-off error, provided that two conditions are met:
Using raxisa_c and axisar_c
/*
Let r(t) be a time-varying rotation matrix; r could be
a C-matrix describing the orientation of a spacecraft
structure. Given two points in time t1 and t2 at which
r(t) is known, and given a third time t3, where
t1 < t3 < t2,
we can estimate r(t3) by linear interpolation. In other
words, we approximate the motion of r by pretending that
r rotates about a fixed axis at a uniform angular rate
during the time interval [t1, t2]. More specifically, we
assume that each column vector of r rotates in this
fashion. This procedure will not work if r rotates through
an angle of pi radians or more during the time interval
[t1, t2]; an aliasing effect would occur in that case.
If we let
r1 = r(t1)
r2 = r(t2), and
-1
q = r2 * r1 ,
then the rotation axis and angle of q define the rotation
that each column of r(t) undergoes from time t1 to time
t2. Since r(t) is orthogonal, we can find q using the
transpose of r1. We find the rotation axis and angle via
raxisa_c.
*/
mxmt_c ( r2, r1, q );
raxisa_c ( q, axis, &angle );
/*
Find the fraction of the total rotation angle that r
rotates through in the time interval [t1, t3].
*/
frac = ( t3 - t1 ) / ( t2 - t1 );
/*
Finally, find the rotation delta that r(t) undergoes
during the time interval [t1, t3], and apply that rotation
to r1, yielding r(t3), which we'll call r3.
*/
axisar_c ( axis, frac * angle, delta );
mxm_c ( delta, r1, r3 );
Constructing a coordinate axis rotation matrix
[ angle ] ,
i
which corresponds to a rotation about a coordinate axis. This is a
special case of the problem solved by axisar_c. Note however that the
matrix produced by rotate_c is the inverse of that produced by axisar_c,
if both routines are provided with the same input angle, and axisar_c is
given the ith coordinate basis vector as the rotation axis.
The call
rotate_c ( angle, i, m );produces `m', the desired matrix. Rotation derivatives
Differentiating rotations
[ angle ]
1
which 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_c is useful for differentiating rotations that are defined by time-varying Euler angles. For example, if the rotation `r' is defined by
r = [ Twist ] [ pi/2 - Dec ] [ pi/2 + RA ]
3 1 3
where RA, Dec, and Twist are time-dependent, then if we make the
abbreviations
A(Twist) = [ Twist ]
3
B(Dec) = [ pi/2 - Dec ]
1
C(RA) = [ pi/2 + RA ]
3
we 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) dt
The derivatives of A, B, and C can be found using drotat_c.
State transformations
p (t)
I
be a position vector referenced to an inertial frame ``I,'' and let
p (t)
N
be the equivalent position vector referenced to a non-inertial 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 I
Therefore, 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 I
It'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
body-centered frame that uses a body-center-to-Sun 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 body-centered frame.
CSPICE provides routines to transform states between inertial frames and body-fixed planetocentric frames. The routine tisbod_c returns the 6x6 transformation matrix required to transform inertially referenced state vectors to body-fixed planetocentric coordinates. If `ref' is the name of the inertial frame of interest, `body' is the NAIF integer code of a body defining a body-fixed planetocentric frame, and `et' is ephemeris time used to define the body-fixed frame, then the call
tisbod_c ( ref, body, et, tsipm );returns `tsipm', the desired 6x6 state transformation matrix. A state vector `s' can be transformed to the body-fixed state vector `sbfixd' by the function call
mxvg_c ( tsipm, s, 6, 6, sbfixd );Since the inverse of a state transformation matrix is not simply its transpose, CSPICE provides the utility routine invstm_c to perform the inversion. If `m' is a state transformation matrix, the inverse matrix `minv' can be obtained via the function call
invstm_c ( m, minv ); Validating a rotation matrix
/*
Set values for the column norm and determinant tolerances
ntol and dtol:
*/
ntol = 1.e-7;
dtol = 1.e-7;
if ( ! isrot_c ( m, ntol, dtol ) )
{
[perform error handling]
}
else
{
m2q_c ( m, q );
.
.
.
}
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 R
is 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 two-dimensional diagram below shows this decomposition. All of the vectors lie in the plane containing r and n.
\ . rParallel
\ .
\ .
r \ .
\ .
\ .
\ .
\ ^
\ |
\ | n
\ |
..........\|
rPerp
Now, 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 three-dimensional rotations are very much like two-dimensional rotations: The effect of a three-dimensional 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) N
allows 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. Time-dependent coordinate transformations
B2 = { v1(t), v2(t), v3(t) }.
An example of a time-dependent 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 k
then
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 k
Since 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 time-dependent 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 Jupiter-Sun 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 Jupiter-Sun 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 x-axis, 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 i3
The 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 > A
will 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 skew-symmetric and N squared is symmetric, so
T
R - R = 2 sin(theta) N.
As long as sin(theta) is non-zero, 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 s
to 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 re-write 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 * y
for 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 * x
is 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 = 1
or
*
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 i3
There 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 (i1-i2-i3) of axes, there is a set of functions that give us w1, w2, and w3 as a function of M:
w1 = f1 ( M ),
i1-i2-i3
w2 = f2 ( M ),
i1-i2-i3
w3 = f3 ( M ).
i1-i2-i3
How can we derive the functions
f1 , f2 , f3 ?
i1-i2-i3 i1-i2-i3 i1-i2-i3
One 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 component-wise correspondence. In subsequent sections, we
actually carry out this procedure for 3-1-3 and 1-2-3 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
3-1-3 3-1-3 3-1-3
that allow us to factor rotation matrix M as a 3-1-3 product:
If
M = [w1] [w2] [w3] . (2)
3 1 3
then
w1 = f1 ( M ),
3-1-3
w2 = f2 ( M ),
3-1-3
w3 = f3 ( M ).
3-1-3
We'd like to somehow use the functions we've already got to factor M as
a 2-3-2 product: we want to find functions
f1 , f2 , f3
2-3-2 2-3-2 2-3-2
such that
M = [y1] [y2] [y3] , (3)
2 3 2
and
y1 = f1 ( M ),
2-3-2
y2 = f2 ( M ),
2-3-2
y3 = f3 ( M )
2-3-2
without having to derive
f1 , f2 , f3
2-3-2 2-3-2 2-3-2
from scratch.
We'll start out by using a new basis, relative to which the right hand side of (3) is not a 2-3-2, but rather a 3-1-3 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]
2
expressed relative to B1 represents the same rotation as the matrix
[y]
3
expressed 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 2
expressed relative to B1 represents the same rotation as the matrix
M' = [y1] [y2] [y3]
3 1 3
expressed 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
3-1-3 3-1-3 3-1-3
to factor M' as a 3-1-3 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 3-1-3 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 3-1-3 factorization with respect to that basis. It is important the new basis be right-handed; otherwise the form of the matrices
[w]
i
is not preserved. It turns out that for any axis sequence of the form
a-b-a, we can find a right-handed basis such that the factorization we
want is a 3-1-3 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 right-handed 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 1-2-1, 2-3-2, and 3-1-3. 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 right-handed 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 1-3-1,
2-1-2, and 3-2-3. So we can use our 3-1-3 formula to handle all of the
a-b-a factorizations, just by computing the correct transformation
matrix C.
What about a-b-c 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 1-2-3 factorization. We'd like to find a basis that allows us to represent the a-b-c product
[w1] [w2] [w3]
a b c
as the product
[w1] [w2] [w3] .
1 2 3
Again, 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 right-handed basis we want is
B2 = { e(a), e(b), e(c) }.
With respect to the basis B2, our a-b-c factorization is a 1-2-3
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 c
to
[w1] [w2] [-w3] ,
1 2 3
so we can still recover the angles w1, w2, and w3 easily. So our 1-2-3
factorization formula allows us to handle all the a-b-c factorizations.
Having shown that we can perform all of the a-b-a and a-b-c 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: a-b-a 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: a-b-c 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
|