Abstract

Introduction

This document covers the CSPICE routines that deal with rotations and
the mathematical ideas behind the routines.

There are three chapters:

-- CSPICE routines

The ``CSPICE routines'' chapter tells you what routines are available, what
they do, and how to call them.

-- Tutorial introduction to rotations

If rotations are new to you, you'll profit most by reading the ``tutorial''
chapter first. This chapter lists most of the facts about rotations used in
the CSPICE code. The emphasis is on building intuition about rotations;
proofs of any noticeable difficulty or length are deferred to the
``Mathematical Road Map'' chapter.

-- Mathematical road map

This chapter contains detailed explanations of a number of the mathematical
ideas used in the CSPICE rotation routines. And as we've said, proofs of
some of assertions made in the tutorial chapter are stowed here.

Using this document

For many readers, the first chapter will be the only one of interest.
The mechanics of using the rotation routines in your code are covered
here.

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.

Categories of functions

The discussions of the functions are categorized according to the type
of problem that the functions solve. This chapter contains one section
for each category. The categories are:

-- Euler angles

-- Quaternions

-- Rotating vectors and matrices

-- Rotation axis and angle

-- Rotation derivatives

-- Validating rotation matrices

The rotation functions constitute a ``family'' of functions insofar as
they deal with related problems, but they do not constitute a
``system.'' So, there is no logical interdependence between the
discussions of functions in different categories.

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

There are several pitfalls associated with converting matrices to Euler
angles. First, any mapping from matrices to Euler angles has
singularities. These come in two flavors: some matrices don't map to a
unique set of Euler angles, and some matrices have the property that a
small change in the matrix can result in a large change in the
corresponding Euler angles.

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

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.

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.

Multiplying quaternions

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'.

Rotating vectors and matrices

This set of routines are frequently used as utilities from which more
complicated routines may be constructed. For example, the routine
eul2m_c constructs a rotation matrix from a sequence of three rotations
about specified coordinate axis, so there is no need to call rotmat_c
directly to accomplish this.

A word of warning

Some care is needed when dealing with signs of rotation angles: a
rotation of a vector by an angle `theta' can be viewed as rotating the
coordinate system by -`theta'. We try to avoid confusion here by
referring to routines as ``coordinate system rotations'' or ``vector
rotations,'' depending on whether a positive rotation angle corresponds
to rotating the coordinate system by a positive angle, or to rotating a
vector by a positive angle. The same criterion applies to matrix
rotations. According to this classification, rotvec_c and rotmat_c
perform coordinate rotations, and vrotv_c performs a vector rotation.

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

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.

raxisa_c and axisar_c can be used to perform linear interpolation
between two rotation matrices. Here's a code fragment illustrating the
procedure:

/*
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

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 rotation derivative routines are utilities that simplify finding
derivatives of time-varying coordinate transformations. In particular,
these routines are used to transform state vectors between non-inertial
reference frames.

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

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

Validating a rotation matrix

isrot_c is a logical function that indicates whether a matrix is a valid
rotation matrix. The criteria for validity are:

-- The columns of the matrix are unit vectors, within a specified tolerance.

-- The determinant of the matrix formed by unitizing the columns of the input
matrix is 1, within a specified tolerance. This criterion ensures that the
columns of the matrix are nearly orthogonal, and that they form a
right-handed basis.

Tutorial introduction to rotations

This tutorial is intended to bridge the gap between knowing linear
algebra and understanding rotations. If you haven't reached the gap yet,
consult reference [2], or any reasonable textbook on the subject.

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

If you're going to read the rest of this tutorial, you're going to see a
lot of definitions, symbol manipulation, and proofs. This information
may be more accessible if you have some kind of reasonable mental model
of rotations; then you can test our claims against your model.

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.

Definition of ``rotation''

Actually, we present three definitions of the term ``rotation.'' What
for? Having more than one way of knowing that a mapping is a rotation
makes it easier to check whether any particular mapping is a rotation or
not. Some properties of rotations are easier to derive from one
definition than from another.

Uses of the definitions

Definition (1) is useful for checking that a mapping is a rotation,
because there's not much to check.

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

A ``rotation matrix'' M is a 3 by 3 matrix whose columns form an
orthonormal set, and whose third column is the cross product of the
first two. Given a rotation R and an orthonormal basis B, the matrix
representation of R relative to B is a rotation matrix.

Any matrix whose columns form an orthonormal set is called an
``orthogonal'' matrix.

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.

Composition of rotations

If we multiply a sequence of rotation matrices, is the result a rotation
matrix? Equivalently, if we compose a sequence of rotations, is the
resulting mapping a rotation?

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:

Coordinate transformations

Change-of-basis transformations between right-handed, orthonormal bases
are rotations. You can verify this using our first definition of
rotations. In particular, the inertial coordinate system transformations
available in CSPICE are all rotations.

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

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

How do we rotate a two-dimensional vector v by theta radians? Trivial as
this problem may seem, it's probably worth your while to get a firm grip
on its solution. When you've understood it, you've almost understood
three-dimensional rotations.

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

Suppose we have a rotation R, defined in three-dimensional space, that
rotates vectors by angle theta about unit vector n. For an arbitrary
vector r, we can break r up into two orthogonal components: one parallel
to n and one perpendicular to n. We can call these components rParallel
and rPerp.

The two-dimensional diagram below shows this decomposition. All of the
vectors lie in the plane containing r and n.

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.

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:

-- The trace of the above matrix is 1 + 2 cos(theta). Since the trace of a
matrix is invariant under similarity transformations, every rotation matrix
has trace equal to 1 + 2 cos(theta). So we can easily find the rotation
angle of any rotation matrix.

-- Every rotation has 1 as one of its eigenvalues. We already knew that, but
there it is, sitting alone in its one-dimensional diagonal block. The other
eigenvalues are complex unless the rotation angle is a multiple of pi.

-- Every rotation is an orthogonal mapping, that is, orthogonal vectors map to
orthogonal vectors. This has to be true because the canonical form is an
orthogonal matrix, represented relative to an orthonormal basis.

Rotation axis and angle

Our soccer ball model shows that a rotation has a fixed vector, called
the ``axis.'' Now if the vector n is fixed by R, then -n is fixed as
well, so the direction of the rotation axis is not unique.

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.

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

Given a rotation matrix M (usually representing a coordinate
transformation), we occasionally have the need to express it as the
product

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.

Quaternion arithmetic

CSPICE does not currently contain any routines that make use of
quaternion arithmetic. Eventually, CSPICE may use quaternion
multiplication to compose rotations, since quaternion arithmetic is more
efficient than matrix multiplication. At present, this section is merely
for the curious.

Definitions

There are two binary operations defined on quaternions: addition and
multiplication.

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.

Deducing the multiplication formula

One really interesting fact about the product formula is that it is a
sum of binary operations, each of which is linear (where the
coefficients are scalars, not quaternions) in both operands. This
implies that the product formula itself is linear in q1 and q2.

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

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

There is one last assertion to check: we've said that you can carry out
composition of rotations using quaternion multiplication. Let's examine
what that means:

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 purpose of this chapter is to familiarize you with the mathematical
ideas essential to dealing with rotations. If you understand the
relevant mathematics, you are in a position to judge the merits of
alternative software designs based on CSPICE routines. If you don't
understand the mathematics, you can still build programs that work by
paying careful attention to function interface specifications, but the
design process is more error-prone, and you're unlikely to hit upon
efficient and elegant solutions.

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.

Formation of a rotation matrix from axis and angle

In this section, we derive an expression for a rotation matrix that
explicitly relates the matrix to the rotation axis and angle. This
expression is valuable for understanding how to find the rotation axis
and angle of a rotation matrix, as well as how to build a rotation
matrix having a given rotation axis and angle. The problem of finding a
quaternion corresponding to a specified rotation matrix is also solved
by the expression derived here.

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

Formation of a rotation matrix from a quaternion

Since the quaternion gives us a rotation's axis and angle, an earlier
discussion in this chapter gives us one way of recovering the rotation
matrix: twice the arccosine of the first component of the quaternion
gives us the rotation angle, and the rest of the quaternion is the
rotation axis, so AXISAR can be used to form the matrix. In this
approach, we may want to treat small rotation angles as a special case,
since the arccosine function is very inaccurate when the argument is
close to 1. We would use the norm of the ``vector'' portion of the
quaternion to give us the sine of half the rotation angle instead, and
recover the rotation angle from this.

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.

Equivalence of rotation definitions

The idea discussed here is used implicitly throughout the CSPICE
rotation routines.

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

We observe that the characteristic polynomial of a rotation is of degree
three, and so has either zero or two complex roots, hence at least one
real root. Because rotations preserve norms, the magnitudes of all of
the roots (eigenvalues), real or complex, are equal to one. So the real
roots are 1 or -1. The determinant of any rotation matrix is 1, since
the determinant of any 3 by 3 matrix is the dot product of the third
column with the cross product of the first and second columns, and for
rotations, this cross product is the third column. But the determinant
is also the product of the eigenvalues. In the case where all three
roots are 1 or -1, we cannot get a product of 1 unless at least one
eigenvalue is equal to 1. If there are complex roots, they are complex
conjugates, so their product is 1, which implies that the real root must
be 1 as well, if the product of all three roots is 1.

A geometric approach

Again, we assume that our rotation R satisfies definition (1), and we
prove that R has a fixed axis.

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.

Recovery of Euler angles from a rotation matrix

Here's the problem: Given a rotation matrix M, and a set of coordinate
axes indexed by i1, i2, i3, find angles w1, w2, w3 such that

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:

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

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.

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