ME 596 Spacecraft Attitude Dynamics and Control

Attitude Kinematics

Professor Christopher D. Hall

Mechanical Engineering

University of New Mexico

October 18, 2021

Dynamics = Kinematics + Kinetics

Translational dynamics (Newton’s 2nd Law)

— 2nd order includes kinematics and kinetics
m~¨r = ~f
— Rewriting as 1st order equations separates the two

kinematics
~r˙ = ~p/m
kinetics
~p˙ = ~f

In applying Newton’s 2nd Law, the reference origin must be

inertial, meaning non-accelerating.
2 / 63
Dynamics = Kinematics + Kinetics

Rotational dynamics (Euler’s Law)

— Includes both kinematics and
kinetics
~h˙ = ~g
— ~h is the angular momentum
vector
— ~g is the applied torque vector
I One of the more difficult
— Rewriting as 1st order equations elements of modeling rotational
separates kinematics and kinetics motion is the connection
between orientation of body and
kinematics angular momentum
To Be Developed I Rotational kinematics is
kinetics
sufficiently important that we
ω̇ = −I−1 ω× Iω + I−1 g
cover it in some detail before
discussing rigid body motion

3 / 63
Translational vs Rotational

This table relates the kinetic and kinematic relationships for

translational and rotational motion.
Linear momentum Angular momentum
= =
mass × velocity inertia × angular velocity

d/dt (linear momentum) d/dt (angular momentum)

= =
applied forces applied torques

d/dt (position) d/dt (angular momentum)

= =
linear momentum/mass angular momentum “/” inertia

4 / 63
Back to Reference Frames
Denote frames as triads of mutually orthogonal unit vectors:
ECI Orbital Body-Fixed
{î1 ,î2 ,î3 } {ô1 ,ô2 ,ô3 } {b̂1 ,b̂2 ,b̂3 }
Fi Fo Fb

Orthogonality means unit vectors are perpendicular to each other:

î1 · î1 = 1 î1 · î2 = 0 î1 · î3 = 0
î2 · î1 = 0 î2 · î2 = 1 î2 · î3 = 0
î3 · î1 = 0 î3 · î2 = 0 î3 · î3 = 1
may be written more concisely as

1 if i=j
îi · îj =
0 if i 6= j
or even more concisely as
îi · îj = δij Kronecker delta
5 / 63
Right-Handedness

The right-handedness of the unit vectors defines the order by:

î1 × î1 = ~0 î1 × î2 = î3 î1 × î3 = −î2

î2 × î1 = −î3 î2 × î2 = ~0 î2 × î3 = î1
î3 × î1 = î2 î3 × î2 = −î1 î3 × î3 = ~0

This set of rules may be written more concisely as

îi × îj = εijk îk

where εijk is the permutation symbol, defined as


 1 for i, j, k an even permutation of 1,2,3
εijk = −1 for i, j, k an odd permutation of 1,2,3
0 otherwise (i.e., if any repetitions occur)


The right-handedness can also be written as

 
 ~0 î3 −î2 
  T  ×
î × î = −î3 ~0 î1 = − î

î2 −î1 ~0 

6 / 63
Vectors
A vector, ~v, is an abstract mathematical object with two
properties: length or magnitude, ||~v||, and direction
~v = v1 î1 + v2 î2 + v3 î3
Scalars, v1 , v2 , and v3 , are the components of ~v expressed in Fi .

The components are dot products of ~v with the base vectors of Fi .

Specifically, v1 = ~v · î1 , v2 = ~v · î2 , v3 = ~v · î3

Since î vectors are unit vectors, components may be written as

v1 = v cos α1 , v2 = v cos α2 , v3 = v cos α3
where v = k~vk, αj is angle between ~v and îj for j = 1, 2, 3.

We frequently collect components of ~v into a column matrix v:


v1


v =  v2 
v3
7 / 63
Vector Notation Summary

Throughout the course, we will work with vectors and with their
representations in several reference frames.

I ~v is a vector, v = ||~v|| is its magnitude

I v1 = v cos α1 = ~v · î1 is the “1” component
I v is a 3 × 1 matrix of components in unspecified reference
frame
I vi is a 3 × 1 matrix of components in Fi
I vb is a 3 × 1 matrix of components in Fb
I vo is a 3 × 1 matrix of components in Fo
no n o
I ~v = viT î = vbT b̂ = voT {ô}

The vector is an abstract mathematical object with direction and

magnitude; its matrix representation is a 3 × 1 matrix of scalars.

These two mathematical objects are not the same thing.

8 / 63
Rotations

I Suppose we know components in body frame, vb

I We want to know components in inertial frame, vi
no n o
~v = viT î = vbT b̂
I Frames are related by a 3 × 3 Rotation Matrix
no n o
î = R b̂
n o no
I We want to find R, so substitute R b̂ for î
n o n o
I Substitution leads to viT R b̂ = vbT b̂
n o
I Since the b̂ vectors are mutually orthogonal,

viT R = vbT
I Transpose both sides to obtain RT vi = vb
I So RT transforms vectors from Fi to Fb
9 / 63
Rotations (continued)

I The equation RT vi = vb is a linear system of the form

Ax = b
I Thus, to determine the components in the inertial frame, we
need to determine R and solve the linear system
I Can write components of R as
î1 = R11 b̂1 + R12 b̂2 + R13 b̂3
î2 = R21 b̂1 + R22 b̂2 + R23 b̂3
î3 = R31 b̂1 + R32 b̂2 + R33 b̂3
I By inspection, the elements Rij are the “direction cosines” of
the unit vectors of Fi with respect to Fb
I That is, Rij = îi · b̂j = cos αij
I This result is why Rotation Matrices are also called Direction
Cosine Matrices (DCM)

10 / 63
Rotations (continued)

I Can also write the rotation matrix as the dot product of two
“vectrices”: n o n oT
R = î · b̂
I Can show that the inverse of a rotation matrix is simply its
transpose
R−1 = RT
I Thus a rotation matrix is an orthonormal matrix: its rows and
columns are components of mutually orthogonal unit vectors
I Introduce superscripts bi to denote rotation from Fi to Fb
vb = Rbi vi
Think of the “i” superscript on R as cancelling with the “i”
subscript on v
T
I Similarly vi = Rib vb , where Rib = Rbi

11 / 63
Rotation Notation

As dot product of “vectrices”:

n o n oT n o n oT
Rbi = b̂ · î Rib = î · b̂

As matrix transforming vectors from one frame to another

vb = Rbi vi vi = Rib vb

As matrix with rows and columns being unit vectors of one frame
expressed in the other
 
îT
 1b
h i
ib T T T
R =  îT = b̂ b̂ b̂

2b  1i 2i 3i
îT
3b

The matrix Rbi is the rotation from Fi to Fb

12 / 63
An Example

The reference frame Fb has the following unit vectors in Fi . What are the
rotation matrices Rbi and Rib ?

b̂1 = 0.8218 î1 − 0.4063 î2 − 0.3994 î3

b̂2 = 0.2490 î1 + 0.8867 î2 − 0.3896 î3
b̂3 = 0.5125 î1 + 0.2207 î2 + 0.8299 î3

The direction cosines can be identified by inspection. For Rbi , the components
of the b̂ vectors in Fi go in the rows.

−0.4063 −0.3994
" #
0.8218
bi
R = 0.2490 0.8867 −0.3896
0.5125 0.2207 0.8299

For Rib , the components of the b̂ vectors in Fi go in the columns. That is,
T
Rib = Rbi " #
0.8218 0.2490 0.5125
ib
R = −0.4063 0.8867 0.2207
−0.3994 −0.3896 0.8299

13 / 63
Attitude Kinematics Representations

I The rotation matrix represents the attitude

I A rotation matrix has 9 numbers, but they are NOT
independent
I There are 6 constraints on the 9 elements of a rotation
matrix (what are they?)
I Thus rotation has 3 degrees of freedom
I There are many different sets of parameters that can be
used to represent or parameterize rotations
I Euler angles, Euler parameters (aka quaternions),
Rodrigues parameters (aka Gibbs vectors), Modified
Rodrigues parameters,

14 / 63
Aircraft Roll, Pitch, and Yaw

15 / 63
Spacecraft Roll, Pitch, and Yaw

Note: RPY Axes can vary significantly from spacecraft to

spacecraft, depending on specifics of a spacecraft’s mission.

16 / 63
Ship Roll, Pitch, and Yaw

Illustration from on-line notes of T. I. Fossen, author of Guidance and Control

of Ocean Vehicles, Wiley, 1994
17 / 63
 Euler Angles

Leonhard Euler (1707-1783) Let’s consider the rotation from

reasoned that rotation from one Fi to Fb using three Euler angles
frame to another can be θ1 , θ2 , θ3
visualized as a sequence of three
We can make the first rotation
simple rotations about base
about any of the three î axes;
vectors
choose î3
Each rotation is through an angle
The first rotation about the î3
(Euler angle) about a specified
axis is through angle θ1
axis
The resulting intermediate frame
Thus a simple rotation is n o
identified by the angle θ and the is denoted Fi0 or î0
axis of rotation
The rotation matrix from Fi to Fi0 is denoted
0 0
R3 (θ1 ) = Ri i ⇒ vi0 = Ri i vi
In the next few slides, we continue the development of Rbi
18 / 63
Fi looking down the î3 axis
The first rotation is R3 (θ1 ): subscript 3 denotes rotation about
the “3” axis, and subscript 1 on θ denotes that this is the first
rotation in the three-Euler Angle rotation sequence.

19 / 63
Rotation from Fi to Fi0

Remember the right-hand rule: counterclockwise rotation about î3 .

20 / 63
Rotation from Fi to Fi0 Direction Cosines
By inspection, we can write down the Fi0 vectors in terms of θ1
and the Fi vectors, noting that î03 = î3

21 / 63
Rotation from Fi to Fi0 Rotation Matrix

Since î3 and î03 are in the same direction, we can easily write the
0
rotation matrix R3 (θ1 ) = Ri i

Note that the minus sign in R3 (θ1 ) is in the row above the row
with a 1 in it.

Next we have to perform two more rotations to get from Fi0 to Fb

22 / 63
 Second Rotation: Fi0 to Fi00

Second rotation can be about As with first rotation, we can

either “1” or “2” axis write down R2 (θ2 ) by inspection
Choose “2” axis, î02 , through θ2 cos θ2 0 − sin θ2
 

R2 (θ2 ) =  0 1 0 
Denote resultingnintermediate
o sin θ2 0 cos θ2
frame as Fi00 or î00
Note that the minus sign in
Rotation matrix is R2 (θ2 ) is in the row above the
R2 (θ2 ) = Ri
00 i0
⇒ vi00 = Ri
00 i0
vi0 row with a 1 in it.
Now we can rotate from Fi to
Note well: Rotation matrix Fi00
notation for “simple”rotations is
Ri (θj ), where subscript i denotes vi00 = R2 (θ2 )vi0 = R2 (θ2 )R3 (θ1 )vi
ith axis and subscript j denotes
jth rotation in three-Euler Angle
Keep in mind that we are
sequence
building Rbi .

23 / 63
 Third Rotation: Fi00 to Fb
Third rotation can be about Note that minus sign in R1 (θ3 ) is
either “1” or “3” axis in bottom row, which should be
viewed as being “above” row with 1.
Choose “1” axis, î001 , through θ3
Note also that minus sign in all
Denoten resulting
o “final” frame as cases is on sine.
Fb or b̂
The pattern is:
Rotation matrix is 1) 1 on diagonal corresponding to
axis;
00 00
R1 (θ3 ) = Rbi ⇒ vb = Rbi vi00 2) 0 on off-diagonal elements of row
and column with 1;
By inspection, the rotation 3) cosines on remaining diagonal
matrix is elements;
  4) sines on remaining off-diagonal
1 0 0
elements;
R1 (θ3 ) =  0 cos θ 3 sin θ3 
5) minus sign on sine in row above
0 − sin θ3 cos θ3 1.

vb = R1 (θ3 )vi00 = R1 (θ3 )R2 (θ2 )R3 (θ1 )vi = Rbi vi

24 / 63
An Example
Reference frame Fa is transformed into reference frame Fb by performing a
3-1-2 Euler rotation sequence through the angles θ1 , θ2 , and θ3 .

a) Determine the direction cosine matrix Rba (ie, from Fa to Fb ).

b) Derive expressions that allow you to determine θi for Rba .

a) The direction cosine matrix Rba is

Rba = R2 (θ3 )R1 (θ2 )R3 (θ1 )
− sin θ3
"
cos θ3
#" #"
cos θ1 sin θ1
#
0 1 0 0 0
= 0 1 0 0 cos θ2 sin θ2 − sin θ1 cos θ1 0
sin θ3 0 cos θ3 0 − sin θ2 cos θ2 0 0 1
cθ3 cθ1 − sθ3 sθ2 sθ1 −sθ3 cθ2
"
cθ3 sθ1 + sθ3 sθ2 cθ1
#
= −cθ2 sθ1 cθ2 cθ1 sθ2
sθ3 cθ1 + cθ3 sθ2 sθ1 sθ3 sθ1 − cθ3 sθ2 cθ1 cθ3 cθ2
where we have used the shorthand c = cos and s = sin.

Note that the order of simple rotation matrices in an Euler angle sequence is
from left to right:
Third rotation × Second rotation × First rotation
25 / 63
An Example (continued): Rotation Matrix to Euler Angles

Part b) is about how we “extract” the Euler angles from a given

rotation matrix
The primary reason for going through this exercise is to illustrate
how tedious it is, as compared with other approaches that we cover
next.

b) Denoting ij element (row i, column j) of Rba as Rij , we can

determine θi as
R21
 
θ1 = tan−1 −
R22
−1
θ2 = sin (R23 )
R13
 
θ3 = tan−1 −
R33

Quadrant checks are imperative. Use atan2(y,x).

26 / 63
An Example (continued): Some numbers
For θ1 = 30◦ , θ2 = −60◦ , θ3 = 45◦ , we have
−0.1768 −0.3536
" #
0.9186
ba
R = R2 (θ3 )R1 (θ2 )R3 (θ1 ) = −0.2500 0.4330 −0.8660
0.3062 0.8839 0.3536
Using the equations for backing out θi , we obtain directly
θ1 = 30◦ , θ2 = −60◦ , and θ3 = 45◦

% Matlab code to do calculations

% note that R1, R2, and R3 are functions that return the
% appropriate 3\times3 rotation matrix (easy to write!)
Rba=R2(45*pi/180)*R1(-60*pi/180)*R3(30*pi/180)
th1=atan2(-Rba(2,1),Rba(2,2))*180/pi
th2=asin(Rba(2,3))*180/pi
th3=atan2(-Rba(1,3),Rba(3,3))*180/pi

Example of the rotation matrix functions:

% R = R1(angle)
% returns elementary "1" rotation matrix
function Rotmat = R1(angle)
c = cos(angle); s = sin(angle);
Rotmat = [ 1 0 0; 0 c s; 0 -s c ];
27 / 63
An Example (continued): Different numbers
For θ1 = 30◦ , θ2 = 120◦ , θ3 = 45◦ , we have
 
0.3062 0.8839 0.3536
Rba = R2 (θ3 )R1 (θ2 )R3 (θ1 ) =  0.2500 −0.4330 0.8660 
0.9186 −0.1768 −0.3536
Using the equations for backing out θi , we obtain directly
θ1 = −150◦ , θ2 = 60◦ , and θ3 = −135◦

Note that these three angles are NOT the same as the three
angles we used to compute Rba

Check by calculating R2 (θ3 )R1 (θ2 )R3 (θ1 ) using these three
angles, and you will find that you obtain the same Rba

What, if anything, went wrong?

Conclusion: Euler angles extracted in this way are “correct,” but

not unique.
28 / 63
How many Euler angle sequences are possible?

I We have considered a couple of different Euler angle sequences

I Answering the question posed above:
— Three choices for first rotation
— Two choices for second rotation, because we cannot rotate
about the same axis
— Two choices for third rotation, because we can rotate about
the same “number” axis as for the first rotation since it is not
really the same axis
— Thus, there are 3 × 2 × 2 = 12 possible Euler angle
sequences:
• 1-2-3, 1-3-2, 2-1-3, 2-3-1, 3-1-2, 3-2-1 (“asymmetric”)
• 1-2-1, 1-3-1, 2-1-2, 2-3-2, 3-1-3, 3-2-3 (“symmetric”)
I The asymmetric set uses all three possible simple rotations
I The symmetric set uses the same rotation for the 3rd rotation
as for the 1st

You need to be able to derive R for any of 12 different sequences.

29 / 63
Roll, Pitch, and Yaw

Roll, pitch and yaw are Euler angles and are sometimes defined as
a 3-2-1 sequence and sometimes defined as a 1-2-3 sequence
What’s the difference?
The 3-2-1 sequence (we did earlier) leads to
−sθ2

cθ1 cθ2 sθ1 cθ2

bi
R =  −cθ3 sθ1 + cθ1 sθ2 sθ3 cθ1 cθ3 + sθ1 sθ2 sθ3 cθ2 sθ3 
cθ1 sθ2 cθ3 + sθ1 sθ3 sθ1 sθ2 cθ3 − cθ1 sθ3 cθ2 cθ3
where θ1 is yaw angle, θ2 is pitch angle, and θ3 is roll angle

The 1-2-3 sequence leads to

cθ2 cθ3 sθ1 sθ2 cθ3 + cθ1 sθ3 sθ1 sθ3 − cθ1 sθ2 cθ3
 

Rbi =  cθ2 sθ3 cθ1 cθ3 − sθ1 sθ2 sθ3 sθ1 cθ3 + cθ1 sθ2 sθ3 
sθ2 −sθ1 cθ2 cθ1 cθ2
where θ1 is roll angle, θ2 is pitch angle, and θ3 is yaw angle

30 / 63
 Roll, Pitch, and Yaw (continued)

Note that the two matrices differ: Rotations do not commute

However, if we assume that angles are small (appropriate for many vehicle
dynamics problems), then approximations of the two matrices are equal
3-2-1 Sequence 1-2-3 Sequence
cos θ ≈ 1 and sin θ ≈ θ cos θ ≈ 1 and sin θ ≈ θ
   
1 θ1 −θ2 1 θ3 −θ2
⇒ Rbi ≈  −θ1 1 θ3  ⇒ Rbi ≈  −θ3 1 θ1 
θ2 −θ3 1 θ2 −θ1 1
θ1 is yaw, θ2 is pitch, and θ3 is roll θ1 is roll, θ2 is pitch, and θ3 is yaw

 ×  ×  ×
θ3 roll θ1
Rbi ≈ 1 −  θ2  ⇒ Rbi ≈ 1 −  pitch  ⇐ Rbi ≈ 1 −  θ2 
θ1 yaw θ3

Thus, roll, pitch, and yaw are the same for 1-2-3 and 3-2-1 sequences, if
the angles are small and a linear approximation is valid.
31 / 63
What You Need To Know About Euler Angles

I Given a sequence, say “3-1-2” for example, derive the

Euler angle representation of R
— Be sure to get the order correct
I Given a sequence and some values for the angles, compute
the numerical values of R
— Be sure to know the difference between degrees and
radians
I Given the numerical values of R, extract numerical values
of the Euler angles associated with a specified sequence
— Be sure to make appropriate quadrant checks, and to
check your answer

32 / 63
Euler’s Theorem

The complexity of using many trig and inverse trig functions leads
us to Euler’s Theorem:

The most general motion of a rigid body with a fixed point is

a rotation about a fixed axis.

The axis, denoted a, is called the eigenaxis or Euler axis.

The angle of rotation, Φ is called the Euler angle or the principal

Euler angle.

Note that simple rotations used in developing Euler angles satisfy

Euler’s Theorem.

Rbi is a rotation matrix, so has Euler axis a and Euler angle Φ.

Rbi a = a Rbi = cos Φ1 + (1 − cos Φ)aaT − sin Φa×
(1 is the identity matrix)

33 / 63
Observations Regarding a and Φ

Since the rotation is about a,

Rbi a = a
from which we can deduce that a is an eigenvector of the rotation
matrix, and that the associated eigenvalue is +1.

Thus every rotation matrix has an eigenvalue equal to +1.

This fact justifies the term eigenaxis for the Euler axis.

This parameterization requires four parameters: the three elements

of a and the scalar Φ.

Given (a, Φ), we can calculate R:

R = cos Φ1 + (1 − cos Φ)aaT − sin Φa×
Next we need to be able to calculate (a, Φ) from a given R.

34 / 63
Extracting a and Φ from R

Just as we need to be able to compute Euler angles from a given

rotation matrix, we need to be able to compute the Euler axis and
Euler angle from R.

The following relationships are fairly straightforward to derive (but

you are not responsible for the derivation):
1
 
−1
Φ = cos (trace R − 1)
2
1  T 
a× = R −R
2 sin Φ

Exercise: Use Euler angles to calculate a rotation matrix, then

extract (a, Φ), and use (a, Φ) to calculate R(a, Φ).

What can you say about the sin Φ = 0 case?

35 / 63
Another Four-Parameter Set: Quaternions

The (a, Φ) representation also suffers from the complexity

assoicated with its trig function terms.

Perhaps the most common rotation representation used for

spacecraft attitude dynamics is the quaternion.

The Euler parameter set, or quaternion, is a four-parameter set

with some advantages over the Euler axis/angle set.

This set of variables can be written in terms of (a, Φ):

Φ
q = a sin
2
Φ
q4 = cos
2
The vector component, q is 3 × 1, and q4 is the scalar component.
h iT
The quaternion is denoted by q̄ = qT q4 , a unit 4 × 1 matrix.
36 / 63
Quaternions: Calculating R(q̄) and q̄(R)
We know how to calculate R from Euler angles and vice versa.

We know how to calculate R from Euler axis and angle and vice
versa.

Here are the relationships between R and q̄

 
R = q42 − qT q 1 + 2qqT − 2q4 q×

1√
q4 = ± 1 + trace R
2
R − R32
 
1  23
q = R31 − R13 
4q4 R − R
12 21

Exercise: Calculate R(a, Φ), extract q̄ from R, recalculate R(q̄).

Check that q̄ agrees with the equations relating q̄ to (a, Φ) .

These calculations provide excellent methods of checking your

answers.
37 / 63
Summary of Kinematics Notation

We have several equivalent methods of describing attitude or

orientation:

I Rotation matrix = DCM = vectors of one frame expressed

in the other = dot products of vectors of one frame with
those of the other
I Euler angles: 3 × 2 × 2 = 12 different sets
I Euler axis/angle: unit vector and angle
I Euler parameters = quaternions: unit 4 × 1

You must be able to compute one from the other for any given
representation
Next: How does attitude vary with time?

38 / 63
Dynamics = Kinematics + Kinetics

Translational dynamics (Newton’s 2nd Law)

— 2nd order includes kinematics and kinetics
m~¨r = ~f
— Rewriting as 1st order equations separates the two

kinematics
~r˙ = ~p/m
kinetics
~p˙ = ~f

In applying Newton’s 2nd Law, the reference origin must be

inertial, meaning non-accelerating.
39 / 63
Dynamics = Kinematics + Kinetics

Rotational dynamics (Euler’s Law)

— Includes both kinematics and
kinetics
~h˙ = ~g
— ~h is the angular momentum
vector
— ~g is the applied torque vector
I One of the more difficult
— Rewriting as 1st order equations elements of modeling rotational
separates kinematics and kinetics motion is the connection
between orientation of body and
kinematics angular momentum
To Be Developed I Rotational kinematics is
kinetics
sufficiently important that we
ω̇ = −I−1 ω× Iω + I−1 g
cover it in some detail before
discussing rigid body motion

40 / 63
Translational vs Rotational

This table relates the kinetic and kinematic relationships for

translational and rotational motion.
Linear momentum Angular momentum
= =
mass × velocity inertia angular velocity

d/dt (linear momentum) d/dt (angular momentum)

= =
applied forces applied torques

d/dt (position) d/dt (angular momentum)

= =
linear momentum/mass angular momentum “/” inertia

41 / 63
Differential Equations of Kinematics

Given the velocity of a point mass and initial conditions for its
position, we can compute its position as a function of time by
integrating the differential equation:
~r˙ = ~v
We now need to develop the equivalent differential equations for
the attitude when the angular velocity is known.
Sneak Preview
  
0 sin θ3 / cos θ2 cos θ3 / cos θ2 ω1
θ̇ =  0 cos θ3 − sin θ3   ω2  = S−1 ω
1 sin θ3 sin θ2 / cos θ2 cos θ3 sin θ2 / cos θ2 ω3

q× + q4 1
 
Φ̇ = aTω 1
q̄˙ = ω = Q(q̄)ω
2 −qT

1 × Φ × ×
ȧ = a − cot a a ω
2 2
42 / 63
Euler Angles and Angular Velocity

We develop θ̇ = S−1 ω one frame at a time, just as we developed

rotation matrices in terms of Euler angles. (3-2-1 here) 
0

0 0 0
~ i i = θ̇1 î03 = θ̇1 î3
ω ωii0 i = ωii i =  0 
θ̇1

0
 
00 i0 i00 i0 i00 i0
~i
ω = θ̇2 î002 = θ̇2 î02 ω i00 =ω i0 =  θ̇2 
0
 
θ̇3
bi00 bi00 bi00
~
ω = θ̇3 b̂1 = θ̇3 î001 ωb =ω i00 = 0 
0
00 00 i0 0
~ bi = ω
Add the angular velocity vectors: ω ~ bi + ω
~i ~i i
+ω

In words: ω of Fb with respect to Fi = ω of Fb w.r.t. Fi00

+ ω of Fi00 w.r.t. Fi0 + ω of Fi0 w.r.t. Fi
43 / 63
Adding the Angular Velocities
The three angular velocities are expressed in different frames:
00 00 0 0
ωbi
b ωii00 i ωii0 i
00
Remember the notation: ωbi b is the angular velocity of Fb with
respect to Fi00 (superscripts), expressed in Fb (subscript).

Rotate all into same frame in order to add them.

Typically, we use Fb , but we can use any reference frame.

00
h iT
We already have ωbi
b in Fb : θ̇3 0 0

00 0 0
Rotate ωii00 i from Fi00 to Fb , and ωii0 i from Fi0 to Fi00 to Fb .
00 i0 00 00 0
ωib = Rbi ωii00 i
0 00 00 i0 0
ωib i = Rbi Ri ωii0 i
Carry out matrix multiplications and complete addition, keeping in
mind 3-2-1 Euler angle sequence.
44 / 63
Completing the Operation
Carry out the matrix multiplications and add the three results:
    
1 0 0 0 0
00 0
ωib i =  0 cos θ3 sin θ3   θ̇2  =  cos θ3 θ̇2 
0 − sin θ3 cos θ3 0 − sin θ3 θ̇2
   
0
1 0 0 cos θ2 0 − sin θ2 0
ωib i =  0 cos θ3 sin θ3   0 1 0  0 
0 − sin θ3 cos θ3 sin θ2 0 cos θ2 θ̇1
 
− sin θ2 θ̇1
=  cos θ2 sin θ3 θ̇1 
cos θ2 cos θ3 θ̇1
Completing the addition and writing as a matrix expression, we obtain
  
− sin θ2 0 1 θ̇1
ωbi
b =
 cos θ2 sin θ3 cos θ3 0   θ̇2 
cos θ2 cos θ3 − sin θ3 0 θ̇3
This result is the angular velocity of Fb with respect to Fi , expressed in
Fb , also written as ωbi
b = S(θ)θ̇.
45 / 63
Completing the Operation (continued)

We have the angular velocity of Fb with respect to Fi , expressed in Fb :

  
− sin θ2 0 1 θ̇1
ωbi
b =
 cos θ2 sin θ3 cos θ3 0   θ̇2  = S(θ)θ̇
cos θ2 cos θ3 − sin θ3 0 θ̇3

We can invert S(θ) to obtain

θ̇ = = S−1 (θ)ω
  
0 sin θ3 / cos θ2 cos θ3 / cos θ2 ω1
=  0 cos θ3 − sin θ3   ω2 
1 sin θ3 sin θ2 / cos θ2 cos θ3 sin θ2 / cos θ2 ω3
This result is the kinematic differential equation for rotational motion in
terms of 3-2-1 Euler angles.

Note that the right-hand side of the equation goes to infinity as θ2

approaches π/2 or 3π/2.

This kinematic singularity is an important topic.

46 / 63
Kinematic Singularity in Euler Angles Differential Equation

I Note that for 3-2-1 Euler angle set, the Euler rates go to
infinity when cos θ2 → 0
I The reason is that at θ2 = π/2 (or 3π/2) the first and second
rotations are indistinguishable
I For the “symmetric” Euler angle sequences (3-1-3, 2-1-2,
1-3-1, etc) the singularity occurs when θ2 = 0 or π
I For the “asymmetric” Euler angle sequences (3-2-1, 2-3-1,
1-3-2, etc) the singularity occurs when θ2 = π/2 or 3π/2
I This kinematic singularity is a major disadvantage of using
Euler angles for large-angle motion
I For example, for Roll-Pitch-Yaw Euler angles, the singularity is
when the Pitch angle is π/2 or 3π/2; for an airplane, this
condition corresponds to a vertical ascent or descent

47 / 63
What You Need To Know ...

... about kinematics differential equations using Euler Angles

I Given a sequence, such as “3-1-2,” derive the expression

for angular velocity vector of Fb with respect to Fi
expressed in Fb ; i.e., derive S(θ), which is different for
each sequence
I Invert S(θ) to obtain the Euler angle differential equations
I Given a sequence and some initial conditions for the
angles, and an expression for the angular velocity as a
function of time ω(t) integrate the differential equations
numerically to obtain the attitude θ(t)
I Given the numerically integrated values of θ(t), compute
the rotation matrix R to get the attitude of the body
frame with respect to the inertial frame

48 / 63
Euler Axis/Angle and Euler Parameters

Euler axis and angle differential equations

Φ̇ = aTω
1 × Φ × ×

ȧ = a − cot a a ω
2 2
Kinematic singularity occurs when Φ = 0 or 2π
Euler parameter differential equations
1 q × + q4 1
 
q̄˙ = ω = Q(q̄)ω
2 −qT
No singularity!

Later we will learn how to get the angular velocity from the
kinetics differential equation:
ω̇ = −I−1 ω× Iω + I−1 g

49 / 63
Typical Problem Involving Angular Velocity and Attitude

Given initial conditions for the attitude (in any form), and a time
history of angular velocity, compute R(t) or any other attitude
representation as a function of time

Requires integration of one of the sets of differential equations

involving angular velocity

Examples:

Given Euler angle sequence, θ(t0 ), and ω(t), integrate θ̇ = S−1 (θ)
and compute R(θ(t))

Given Euler angle sequence, θ(t0 ), and ω(t), compute R(t0 ),

q̄(t0 ), integrate q̄˙ = Q(q̄)ω, and compute R(q̄(t))

50 / 63
Problem 1

I The orientations of two spacecraft A and B relative to an

inertial frame are given through the 3-2-1 Euler angle sequence:
— θA = [30 − 45 60]T and θB = [10 25 − 15]T degrees
I What is the relative orientation of spacecraft A relative to B
in terms of the 3-2-1 Euler angles?
I In other words, what are θ1 , θ2 , and θ3 for the rotation matrix
RAB ?
I The answer is (359.1, 287.7, 79.96) degrees
I Need to do the matrix multiplication to get the rotation matrix:
T
RAB = RAi RiB , where RiB = RBi
I Easy to do the calculations in Matlab ....

51 / 63
Problem 1 (Matlab Code)

% Kinematics Example Problem 1

% The orientations of two spacecraft A and B relative to an inertial fra
% are given for the 3-2-1 Euler angle sequence:
% theta_A=(30, -45, 60) and theta_B=(10, 25, -15) degrees
%
% What is the relative orientation of spacecraft A relative to B in term
% of the 3-2-1 Euler angles?
% In other words, what are theta_1, theta_2, and theta_3
% for the rotation matrix R^{AB}?
%
% Define the 3x1 matrices of the Euler angles for spacecraft A and B
% Convert from degrees to radians
thetaA = [30; -45; 60]*pi/180;
thetaB = [10; 25; -15]*pi/180;

% Compute the two rotation matrices

RAi = R1(thetaA(3))*R2(thetaA(2))*R3(thetaA(1));
RBi = R1(thetaB(3))*R2(thetaB(2))*R3(thetaB(1));

% Compute RAB = RAi*RiB = RAi*inv(RBi) = RAi*RBi’

RAB = RAi*RBi’;
52 / 63
Problem 1 (Matlab Code continued)

% Kinematics Example Problem 1 CONTINUED

% Extract the Euler angles from RAB

th2 = asin(-RAB(1,3))
th3 = atan2(RAB(2,3),RAB(3,3))
th1 = atan2(RAB(1,2),RAB(1,1))

% Check that these three angles give RAB

RABcheck = R1(th3)*R2(th2)*R3(th1)

% RAB*RABcheck’ should be identity

identitycheck = RAB*RABcheck’

% It works, so no need to go back and do further quadrant checks.

% also nice to see the angles in degrees:
thetaABdeg = [th1; th2; th3]*180/pi

% Since some angles are negative, we could also shift them to be

% between 0 and 360 degrees by simply adding 360 to the negative ones.
% More generally, though, the function mod can do it for us

thetaABdeg = mod(thetaAB,360)
53 / 63
Problem 2

I The attitude of FB relative to the inertial frame is given

through the 3-2-1 Euler angle sequence
— θB = [10 25 − 15]T degrees
I Find the principal rotation axis a and angle Φ.
I Calculate the rotation matrix RBi , then
1
 
Φ = cos−1 (trace R − 1)
2
1  T 
a× = R −R
2 sin Φ
— Φ = 0.5546 rad = 31.78◦
— a = [−0.5320 0.7403 0.4110]T
I Easy to do in Matlab ....

54 / 63
Problem 2 (Matlab code)

% Kinematics Example Problem 2

% Define the 3x1 matrix of the Euler angles for spacecraft B
% Convert from degrees to radians
thetaB = [10; 25; -15]*pi/180;

% Compute the rotation matrix

RBi = R1(thetaB(3))*R2(thetaB(2))*R3(thetaB(1));

% Extract the Euler angle and Euler axis from this matrix
cosPhi = (trace(RBi)-1)/2;
Phi = acos(cosPhi);
sinPhi = sin(Phi);
askew = (RBi’-RBi)/(2*sinPhi);
a = [askew(3,2); askew(1,3); askew(2,1)]

% Check that this Euler axis/angle set gives RBi

RBicheck = cosPhi*eye(3)+(1-cosPhi)*a*a’-sinPhi*askew;

% RBi*RBicheck’ should be identity

identitycheck = RBi*RBicheck’

% It works, so let’s show the results, with angle in rad and deg
Phirad = Phi
Phideg = Phi*180/pi

55 / 63
Problem 3

I The orientation of a spacecraft B relative to an inertial frame

is given through the 3-2-1 Euler angle sequence, with initial
conditions
— θB (0) = [30 − 45 60]T degrees
I The attitude control system causes the spacecraft to rotate
with angular velocity ωBi = [1 1 1]T rad/s expressed in FB
I Integrate the equations of motion to obtain θ(t), for t = 0 to
60 seconds
I Easy to do the calculations in Matlab ....
I Matlab includes several functions that are useful for integration
of ordinary differential equations, particularly ode45
— Code right-hand side of differential equations in one file rhs.m
— Code parameters, initial conditions, ode45 call, and post-processing
into driver.m
— Name rhs.m and driver.m with descriptive names for particular
problem
56 / 63
Problem 3 (Plot 1)

Matlab produces the following plot using the default settings. Clearly this plot
is not as useful as it could be. Fortunately, Matlab allows the user to make
modifications to the plots.

57 / 63
Problem 3 (Plot 2)
Changing the linewidth and fontsize parameters, adding labels and legends,
produces a more satisfying plot of the Euler angles vs time.

Observations: 1) Euler angles do not approach kinematic singularity;

2) θ1 exhibits “whirling” behavior, continuous increase from 0 to 2π;
3) θ2,3 exhibit oscillatory behavior

Exercise: Can you explain these observations in terms of the kinematic

differential equation?
58 / 63
Problem 3 (Matlab Code)
% th321driver.m
th0=[30;-45;60]*pi/180; %initial conditions
tspan=linspace(0,60,1000); %integration timespan
options=odeset(’abstol’,1e-9,’reltol’,1e-9); %set tolerances
[t,th]=ode45(’th321rhs’,tspan,th0); %call ode45
th(:,1:2)=mod(th(:,1:2),2*pi); %make angles [0,2pi]
plot(t,th); %plot 1
figure; hg=plot(t,th); %plot 2
set(hg,’linewidth’,3)
set(gca,’fontsize’,48)
hg=xlabel(’t (seconds)’);
set(hg,’fontsize’,48);
hg=ylabel(’\theta(t) (radians)’);
set(hg,’fontsize’,48);
hg=legend(’\theta_1’,’\theta_2’,’\theta_3’);
set(hg,’fontsize’,48);

% th321rhs.m right hand side function called by ode45

function thdot = th321rhs(t,th)
s2=sin(th(2)); %calculate trig f’s
c2=cos(th(2));
s3=sin(th(3));
c3=cos(th(3));
Sinv=[0 s3/c2 c3/c2; 0 c3 -s3; 1 s3*s2/c2 c3*s2/c2]; %Sinverse
thdot=Sinv*[1;1;1]; %dtheta/dt
end
59 / 63
Problem 4

I The orientation of a spacecraft B relative to an inertial frame

is given through the 3-2-1 Euler angle sequence, with initial
conditions:
— θB = [30 − 45 60]T degrees
I The attitude control system makes the spacecraft rotate with
angular velocity ωBi = [1 1 1]T rad/s expressed in FB
I Integrate the quaternion equations of motion to obtain q̄(t),
for t = 0 to 10 seconds
I Easy to do the calculations in Matlab

60 / 63
Problem 4 (Plot)

Observations: The components of the quaternion, (also called quaternions)

exhibit simple oscillatory behavior, always in the interval [−1, 1].

Exercise: Compare the results of Euler angle integration to quaternion

integration, both of which provide equivalent results in terms of the attitude.

61 / 63
Summary
Learning Objective 1. Describe attitude kinematics using reference frames,
rotation matrices, Euler parameters, Euler angles, and quaternions.

I Define reference frame. Define inertial origin. Define inertial

reference frame.
I Derive rotation matrices using direction cosines. ... dot products. ...
Euler angles. ... Euler parameters.
I Derive attitude kinematics differential equations for Euler angles.
I Identify and describe kinematic singularities.

Learning Objective 2. Accurately perform attitude kinematics calculations.

I Calculate rotation matrices using direction cosines. ... dot products.

... Euler angles. ... Euler parameters.
I Given attitude in one set of variables (e.g., Euler angles), calculate
the rotation in a different set of variables (e.g., quaternion).
I Calculate attitude as function of time by integrating the attitude
kinematics differential equations.
62 / 63
Looking Forward

I We need to develop the kinetic differential equations for rigid bodies,

in order to obtain the angular velocity, ωbi (t).
I Knowing the angular velocity, we can determine the attitude as a
function of time.
I Before beginning the Rigid Body Dynamics material, we will cover
the topic of Attitude Determination, which requires
— attitude sensors
— mathematical models of subjects of attitude sensors
— algorithms that use attitude sensor measurements and
mathematical models to determine attitude
I So, the remainder of the course includes Exam 1, Attitude
Determination, Rigid Body Dynamics, Exam 2, Satellite Dynamics,
and Final Exam.

63 / 63