In [2]:

```
from IPython.display import Latex
from IPython.display import Image
from IPython.core.display import HTML
```

- SMAD 19

- Derivation of Transport Equation.
- Derivation of torque-free motion equations.
- Stability analysis for torque-free motion.
- Energy analysis for torque-free motion.
- Historical anecdote about energy dissipation.
- Derivation of equation of motion for spacecraft with reaction wheels/thrusters.
- Discussion of control moment gyroscopes.

- Review of TRIAD for attitude determination.

- Star trackers.
- Sun/Earth sensors.
- Horizon sensors.
- Magnetometers.
- CDGPS.
- MEMS gyros.
- Bugs!
- Fiber-optic gyros.
- Optical navigation.

I want to start this lecture by reviewing how we derive the equation of motion for a spacecraft that may be either uncontrolled or acted upon by reaction wheels, control moment gyroscopes, thrusters, dampers, and other actuators. I believe that you have all seen this content, either in part or in full, in your other classes (4060, and perhaps 3260). That is ok. These equations, and the equations which describe the spacecraft's trajectory, are the foundational equations upon which any spacecraft's attitude determination and control systems are built. It is a good idea to see them again.

Also, because you have seen some of this material before, I am going to review it by presenting it in a way that you perhaps have not seen before. I think this is the appropriate way to important material like this because, at worst, doing so will show you things that you've already seen. At best, presenting it in a new way will improve your understanding of these concepts. *You should get into the habit of looking back at old concepts and ideas as you learn new concepts and ideas. You almost always see them in a new light.*

My goal is to remove every shred of magic or mystery from these equations by deriving them from very, very basic principles. You likely have not seen that entire process. You will have taken classes where you get the mathematical tools for the basic principles (calculus, linear algebra, etc.) and you will have taken classes that use the final equations to draw conclusions about the system. But now that you've done both, you have the ability to look at the entire derivation. We can look at old concepts with new knowledge and (maybe) see them in new light.

**Note**: **There will be aspects of this document that will seem very elementary to you. I'm not including this content because I think you don't know it, I'm including it for completeness.**

There is one equation at the foundation of attitude dynamics, determination, and control. It is the famous Transport Theorem, shown below:

\begin{align}
\frac{{}^Nd}{dt}\textbf{v} &= \frac{{}^{B}d}{dt}\textbf{v} + \omega^{B/N} \times \textbf{v}
\end{align}

This equation tells us that the time derivative of some vector in an inertial (not moving) frame (N) is equal to the time derivative of that vector in a rotating frame (B) plus the cross product between the rotation rate between the two frames and the vector. This equation is fundamental enough that it is often handed to you as something foundational to be remembered (and it is), but let's revisit it with a senior-level background in calculus and physics to try to get a visceral understanding for why this equation must be what it is.

Suppose that we have two frames, one of which is rotating and the other of which is inertial. The inertial frame, $N$, is represented by the unit vectors $\hat{i}$, $\hat{j}$, and $\hat{k}$. To represent the rotating frame, $B$, we'll simply add primes to each unit vector: $\hat{i}'$, $\hat{j}'$, $\hat{k}'$. The $B$ frame is rotating at some constant rate $\omega^{B/N}$, as shown below. In the below image, the $B$ frame happens to be rotating along the $\hat{k}$ axis, but it could be spinning around any arbitrary axis.

In [7]:

```
Image("frames.png", width=500)
```

Out[7]:

Now suppose that we have some (potentially time-varying) vector $\textbf{v}$, which is represented in the rotating $B$ frame, as illustrated below. Perhaps this is the boresight of some instrument, which is known in the spacecraft body frame (which is spinning relative to the inertial frame).

In [6]:

```
Image("addv.png", width=500)
```

Out[6]:

What we would like to know is the time derivative of $\textbf{v}(t)$ in the *inertial* frame.

Let's just go through the steps that we would take to start doing that, and we'll discover that we arrive at the Transport Theorem. Start by simply writing down $\textbf{v}$ in its $B$ frame components, which we know.

\begin{align}
{}^{B}\textbf{v} = v_i' \hat{\textbf{i}}' + v_j' \hat{\textbf{j}}' + v_k' \hat{\textbf{k}}'
\end{align}

And now let's take the total derivative in order to get the expression for the rate of change of $\textbf{v}$ in the inertial frame. Using the Chain Rule:

\begin{align}
\frac{{}^Nd}{dt}\textbf{v} &= \left[ \left(\frac{d}{dt} v_i'\right) \hat{\textbf{i}}' + v_i' \left(\frac{d}{dt} \hat{\textbf{i}}'\right) \right] + \left[ \left(\frac{d}{dt} v_j'\right) \hat{\textbf{j}}' + v_j' \left(\frac{d}{dt} \hat{\textbf{j}}'\right) \right] + \left[ \left(\frac{d}{dt} v_k'\right) \hat{\textbf{k}}' + v_k' \left(\frac{d}{dt} \hat{\textbf{k}}'\right) \right]
\end{align}

Rearranging:

\begin{align}
\frac{{}^Nd}{dt}\textbf{v} &= \left[ \left(\frac{d}{dt} v_i'\right) \hat{\textbf{i}}' + \left(\frac{d}{dt} v_j'\right) \hat{\textbf{j}}' + \left(\frac{d}{dt} v_k'\right) \hat{\textbf{k}}' \right] + \left[ v_i' \left(\frac{d}{dt} \hat{\textbf{i}}'\right) + v_j' \left(\frac{d}{dt} \hat{\textbf{j}}'\right) + v_k' \left(\frac{d}{dt} \hat{\textbf{k}}'\right) \right]
\end{align}

Pause for a moment and think about the above expression. In the case that the $B$ and $N$ frames were not rotating relative to one another, what would happen to the second term? It would go to zero, because the unit vectors for the $B$ frame would be constant (derivatives=0) for the $N$ frame.

What is the first term? The first term is just the time derivative of $\textbf{v}$ in the $B$ frame! So we can rewrite this:

\begin{align}
\frac{{}^Nd}{dt}\textbf{v} &= \frac{{}^{B}d}{dt}\textbf{v} + \left[ v_i' \left(\frac{d}{dt} \hat{\textbf{i}}'\right) + v_j' \left(\frac{d}{dt} \hat{\textbf{j}}'\right) + v_k' \left(\frac{d}{dt} \hat{\textbf{k}}'\right) \right]
\end{align}

All that remains is to inerpret the time derivatives of the rotating $B$ basis vectors. In the most general case, the rate of change of some vector can be decomposed into a radial component (changing length) and a tangential component (changing direction). Since we are dealing with *fixed-length unit vectors*, they have no rate of change in the radial direction (i.e. they are not changing length). So, the derivatives of these vectors are entirely in the tangential direction. We could systematically prove this by differentiating the orthogonality rules, but I don't think that's necessary.

This reduces to phyics that you learned in high school. Think of each unit vector as a position vector (fixed length) undergoing circular motion. What is the equation for the tangential velocity of that vector? It is the angular velocity crossed with the position vector! $v = \omega \times \textbf{r}$. So, substituting into our above expression and simplifying:

\begin{align}
\frac{{}^Nd}{dt}\textbf{v} &= \frac{{}^{B}d}{dt}\textbf{v} + \left[ v_i' \left(\omega^{B/N} \times \hat{\textbf{i}}'\right) + v_j' \left(\omega^{B/N} \times \hat{\textbf{j}}'\right) + v_k' \left(\omega^{B/N} \times \hat{\textbf{k}}'\right) \right]\\
&= \frac{{}^{B}d}{dt}\textbf{v} + \left[ \left(\omega^{B/N} \times v_i'\hat{\textbf{i}}'\right) + \left(\omega^{B/N} \times v_j'\hat{\textbf{j}}'\right) + \left(\omega^{B/N} \times v_k'\hat{\textbf{k}}'\right) \right]\\
&= \frac{{}^{B}d}{dt}\textbf{v} + \omega^{B/N} \times \left[ v_i'\hat{\textbf{i}}' + v_j'\hat{\textbf{j}}' + v_k'\hat{\textbf{k}}' \right]\\
&= \frac{{}^{B}d}{dt}\textbf{v} + \omega^{B/N} \times \textbf{v}
\end{align}

There is nothing mysterious here. The Transport Theorem *must* be what it is, because nothing else would make sense.

In order to understand the rotational motion of a spacecraft, we are interested in the inertial derivative of one specific vector: *the angular momentum*. This is how we arrive at the torques (the rotational analogues for forces) acting on the spacecraft.

\begin{align}
\tau &= \frac{{}^{N}d}{dt}\textbf{H}
\end{align}

Where, you'll recall, the angular momentum is equal to the spacecraft inertial tensor dotted with the angular velocity. I'm using $I_c$ to designate the inertia tensor of the spacecraft, and $\textbf{H}_c$ to denote the angular momentum of the spacecraft.

\begin{align}
\textbf{H}_c &= I_c \cdot \omega^{B/N}
\end{align}

We can use the Transport Theorem that we've just derived in order to get an expression for this torque.

\begin{align}
\tau &= \frac{{}^{N}d}{dt}\textbf{H}_c\\
&= \frac{{}^Bd}{dt}\textbf{H}_c + \omega^{B/N} \times \textbf{H}_c\\
&= \frac{{}^Bd}{dt}\left(I_c \cdot \omega^{B/N}\right) + \omega^{B/N} \times \left(I_c \cdot \omega^{B/N}\right)\\
&= I_c \cdot \left(\frac{{}^Bd}{dt}\omega^{B/N}\right) + \omega^{B/N} \times \left(I_c \cdot \omega^{B/N}\right)
\end{align}

For clarity, I'm going to replace inertial derivatives with a superscript $N$, and body-frame derivatives with a superscript $B$. Rewriting the above:

\begin{align}
\tau = \overset{N}{\textbf{H}}_c = I_c \cdot \left(\overset{B}{\omega}^{B/N}\right) + \omega^{B/N} \times \left(I_c \cdot \omega^{B/N}\right)
\end{align}

We will start by considering the torque-free case (no thrusters, torque coils, etc). In this case, the above expression for the torque is set to 0.

\begin{align}
\tau = \overset{N}{\textbf{H}}_c = I_c \cdot \left(\overset{B}{\omega}^{B/N}\right) + \omega^{B/N} \times \left(I_c \cdot \omega^{B/N}\right) = 0
\end{align}

Recall that the inertia tensor describes the distribution of mass in an object, and that it has certain properties. They will be symmetric and (for physical systems) positive definite. So, the most general form of the above inertia tensor, $I_c$, can be written as:

\begin{align}
I_c &= \begin{bmatrix}
I_{11} & I_{12} & I_{13}\\
I_{12} & I_{22} & I_{23}\\
I_{13} & I_{23} & I_{33}\end{bmatrix}
\end{align}

where $I_{11}$, $I_{22}$ and $I_{33}$ are positive. However! It is *always* possible to diagonalize the inertia tensor. In other words, it is always possible to choose a $B$-frame coordinate system such that the mass of the object is evenly distributed around each axis. More mathematically put, we can align our $B$-frame coordinate system with the eigenvectors of the inertia tensor. Thus, if we assume that we've judiciously chosen our $B$ frame coordinates to align with the eigenvectors of the inertia tensor, then we can rewrite the above as below without loss of generality.

\begin{align}
I_c &= \begin{bmatrix}
I_{11} & 0 & 0\\
0 & I_{22} & 0\\
0 & 0 & I_{33}\end{bmatrix}
\end{align}

We call these the *principle axes*. When we're doing analysis, we almost always do so under the assumption that our non-inertial frames align with the principle axes.

We want to analyze each component of our spacecraft's angular velocity in a torque-free situation. To do so, we can rewrite our Transport Theorem equation in matrix form, so that we can then isolate each component of the angular velocity.

\begin{align}
\overset{N}{\textbf{H}}_c &= I_c \cdot \left(\overset{B}{\omega}^{B/N}\right) + \omega^{B/N} \times \left(I_c \cdot \omega^{B/N}\right)\\
&= \begin{bmatrix}
I_{11} & 0 & 0\\
0 & I_{22} & 0\\
0 & 0 & I_{33}\end{bmatrix} \begin{bmatrix}
\dot{\omega}_1\\
\dot{\omega}_2\\
\dot{\omega}_3
\end{bmatrix}^{B/N}+ \begin{bmatrix}
0 & -\omega_3 & \omega_2\\
\omega_3 & 0 & -\omega_1\\
-\omega_2 & \omega_1 & 0
\end{bmatrix}\left(\begin{bmatrix}
I_{11} & 0 & 0\\
0 & I_{22} & 0\\
0 & 0 & I_{33}\end{bmatrix} \begin{bmatrix}
\omega_1\\
\omega_2\\
\omega_3
\end{bmatrix}^{B/N}\right) = \begin{bmatrix}
0\\0\\0\end{bmatrix}
\end{align}

If we perform all of this matrix multiplication, simplify, and separate into the three equations associated with the 3x1 vector, then we arrive at the following set of equations:

\begin{align}
I_1 \dot{\omega}_1 &= \left(I_2 - I_3\right) \omega_2 \omega_3\\
I_2 \dot{\omega}_2 &= \left(I_3 - I_1\right) \omega_3 \omega_1\\
I_3 \dot{\omega}_3 &= \left(I_1 - I_2\right) \omega_1 \omega_2
\end{align}

**Maximum Axis**: Let's assume that $I_1 > I_2 > I_3$. Consider the situation in which the spacecraft (or whatever object) is rotating about the axis with moment of inertia $I_1$ (it's rotating about the axis with the greatest moment of inertia). Furthermore, assume very small angular velocities in the other two directions ($\omega_2$, $\omega_3$ small). Under the assumption of small $\omega_2$ and $\omega_3$, the first equation above tells us that $\dot{\omega}_1$ is very small and can be neglected. Let's consider what happens to the other two components of the angular velocity.

We'll start with $\omega_2$. Differentiate the second equation above:

\begin{align}
I_2 \ddot{\omega_2} &= (I_3 - I_1) \dot{\omega}_3 \omega_1 + (I_3-I_1) \omega_3 \dot{\omega}_1
\end{align}

Under the assumption of very small $\omega_2$ and $\omega_3$, we know that $\dot{\omega}_1$ is negligibly small. So, the second term may be eliminated. We can substitute the third equation above in for $\dot{\omega}_3$ in the second term:

\begin{align}
\ddot{\omega}_2 &= \frac{(I_3-I_1)(I_1 - I_2)}{I_2I_3} \omega_1^2 \omega_2\\
&= (\text{negative quantity}) \omega_2
\end{align}

What is this an equation for? It's Hooke's Law! This is the equation for a simple harmonic oscillator. Precisely the same analysis shows that we also get a simple harmonic oscillator for $\omega_3$. So, spinning around the maximium axis is *stable*.

**Minimum Axis**: Now the situation in which the spacecraft (or whatever object) is rotating about the axis with moment of inertia $I_3$ (it's rotating about the axis with the smallest moment of inertia). Furthermore, assume very small angular velocities in the other two directions ($\omega_1$, $\omega_2$ small). Under the assumption of small $\omega_1$ and $\omega_2$, the third equation above tells us that $\dot{\omega}_3$ is very small and can be neglected. Let's consider what happens to the other two components of the angular velocity.

We'll again start with $\omega_2$. Differentiate the second equation above:

Under the assumption of very small $\omega_1$ and $\omega_2$, we know that $\dot{\omega}_3$ is negligibly small. So, the first term may be eliminated. We can substitute the first equation above in for $\dot{\omega}_1$ in the second term:

\begin{align}
\ddot{\omega}_2 &= \frac{(I_3-I_1)(I_2 - I_3)}{I_2I_1} \omega_3^2 \omega_2\\
&= \text{(negative quantity}) \omega_2
\end{align}

And, again, we would get a similar expression for $\omega_1$. So, for spin about the minimum axis too, the system is stable. The transverse angular velocities will oscillate (i.e. we have precession), but it is stable.

**Intermediate Axis**: Now the situation in which the spacecraft (or whatever object) is rotating about the axis with moment of inertia $I_2$ (it's rotating about the axis with the intermediate moment of inertia). Furthermore, assume very small angular velocities in the other two directions ($\omega_1$, $\omega_3$ small). Under the assumption of small $\omega_1$ and $\omega_3$, the second equation above tells us that $\dot{\omega}_2$ is very small and can be neglected. Let's consider what happens to the other two components of the angular velocity.

Let's look at $\omega_1$. Differentiate the first equation above.

\begin{align}
I_1 \ddot{\omega_1} &= (I_2 - I_3) \dot{\omega_2}\omega_3 + (I_2 - I_3) \omega_2 \dot{\omega_3}
\end{align}

Eliminate the first term, and substitute for the second term:

\begin{align}
\ddot{\omega_1} &= \frac{(I_2 - I_3)(I_1 - I_2)}{I_1I_3} \omega_2^2 \omega_1 \\
&= (\textbf{positive quantity}) \omega_1
\end{align}

In the case of rotation about the intermediate axis, $\omega_1$ (and, by similar analysis, $\omega_3$) is **unstable** and will therefore grow. Even a small disturbance along other axes will cause the object to flip, as shown in the video below. You can try this yourself too. Try to flip your phone end over end without it rotating along any other axes.

Consider again these equations.

Notice what happens if we premultiply the first, second, and third equations by $\omega_1$, $\omega_2$, and $\omega_3$, respectively, and then sum all three equations:

\begin{align}
\omega_1I_1 \dot{\omega}_1 + \omega_2I_2 \dot{\omega}_2 + \omega_3I_3 \dot{\omega}_3 &= \omega_1\left(I_2 - I_3\right) \omega_2 \omega_3 + \omega_2\left(I_3 - I_1\right) \omega_3 \omega_1+\omega_3\left(I_1 - I_2\right) \omega_1 \omega_2\\
&= \omega_1 \omega_2 \omega_3 \left[I_2-I_3+I_3-I_1+I_1-I_2 \right]\\
&= 0
\end{align}

This is interesting, because it tells us something about the energy in the system. Consider the expression for the rotational kinetic energy of the spacecraft, shown below:

\begin{align}
T &= \frac{1}{2} \omega^{{B/N}^T} \cdot I_c \cdot \omega^{B/N}\\
&= \frac{1}{2}\left[I_2 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2\right]
\end{align}

Differentiate:

\begin{align}
\frac{d}{dt}T &= \frac{1}{2}\left[2I_2 \omega_1 \dot{\omega}_1 + 2 I_2 \omega_2 \dot{\omega}_2 + 2 I_3 \omega_3 \dot{\omega}_3\right]\\
&= \omega_1I_1 \dot{\omega}_1 + \omega_2I_2 \dot{\omega}_2 + \omega_3I_3 \dot{\omega}_3\\
&= 0 \text{ as shown above}
\end{align}

We've just shown that the rotational kinetic energy is a constant. Thus, the following equation must be true:

\begin{align}
T &= \frac{1}{2}\left[I_2 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2\right] = \text{constant}
\end{align}

*This is the equation for an ellipsoid.* Thus, the angular velocity of the spacecraft is constrained to move on the surface of an ellipsoid. But, in fact, there is another constraint. The magnitude of the angular momentum vector is also conserved.

\begin{align}
||\textbf{H}_c||^2 &= \left(I_c \cdot \omega^{B/N}\right)^T \left(I_c \cdot \omega^{B/N}\right)\\
&= \omega^{{B/N}^T} I_c^T I_c \omega^{B/N}\\
&= I_1^2 \omega_1^2 + I_2^2 \omega_2^2 + I_3^2 \omega_3^2
\end{align}

This is the expression for *another* ellipsoid. The angular velocity must satisfy both of these expressions, so it may only travel on the intersections of these two ellipsoids. These intersections are called *polhodes*, and they are either circular or taco-shaped, as shown below.

In [10]:

```
Image("poinsot.png", width=300)
```

Out[10]:

We can visually see the stability of the major and minor axes by looking at this ellipsoid. We see that the angular velocity vector will walk around the minimum/maximum axes, and that there is an infinitesimally small equilibrium point associated with the intermediate axis. We can use the three torque-free equations that we've been analyzing to draw the direction that the vector walk along each of these contours.

The precession of the angular velocity vector around these principle axes is called *nutation*. We see nutation when the spacecraft either has too much energy to spin about the maximum axis, or too little energy to spin about the minimum axis. Which raises the following question:

Suppose that our spacecraft has a small amount of non-rigidity that allows for some energy dissipation. This may come in the form of fuel slosh, floppy antennas, wiggling solar panel arms, etc. We assume that this non-rigidity is small enough that the inertia tensor for the spacecraft is well approximated as constant. How does the spacecraft behave in this case?

Remember, even when energy is being dissipated, the magnitude of the angular momentum vector remains contant in the body frame (again, with no external torques from thrusters, etc.), and the vector itself remains constant in the inertial frame.

\begin{align}
||\textbf{H}_c||^2 &= I_1^2 \omega_1^2 + I_2^2 \omega_2^2 + I_3^2 \omega_3^2
\end{align}

\begin{align}
2T &= I_2 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2
\end{align}

Let's define some new terms that allow us to consider both of the transverse axes as a pair:

\begin{align}
\omega_t &= \left(\omega_1^2 + \omega_2^2\right)^{\frac{1}{2}}\\
I_t &= I_1 = I_2
\end{align}

With the above quantities defined, we can rewrite the expressions for $||\textbf{H}_c||^2$ and $2T$:

\begin{align}
||\textbf{H}_c||^2 &= I_t^2 \omega_t^2 + I_3^2 \omega_3^2\\
2T &= I_t \omega_t^2 + I_3 \omega_3^2
\end{align}

The above expression for the kinetic energy can be rewritten in terms of the angular momentum magnitude:

\begin{align}
2T &= I_t \omega_t^2 - \left(\frac{I_t^2 \omega_t^2 - ||\textbf{H}_c||^2}{I_3}\right)
\end{align}

(If you substitute in for $||\textbf{H}_c||^2$ and simplify, you'll find that you get $2T = I_t \omega_t^2 + I_3 \omega_3^2$, as before). Rearranging:

\begin{align}
2T &= I_t\left(1 - \frac{I_t}{I_3}\right) \omega_t^2 + \frac{||\textbf{H}_c||^2}{I_3}
\end{align}

Now we take the total derivative. Recall that the magnitude of the angular momentum remains constant, even with energy dissipation, so the second term disappears when we differentiate. I will leave the derivative for $\omega_t^2$ unevaluated, since this is what we're interested in solving for:

\begin{align}
2\dot{T} &= I_t\left(1 - \frac{I_t}{I_3}\right) \frac{d}{dt}\left(\omega_t^2\right)
\end{align}

Solve for the derivative of the squared transverse angular velocity:

\begin{align}
\frac{d}{dt} \omega_t^2 &= \frac{2\dot{T}}{I_t \left(1 - \frac{I_t}{I_3}\right)}
\end{align}

Consider the above equation. In the case of energy dissipation, $\dot{T}<0$. If $I_3 > I_t$, then the denominator is positive, and consequently the derivative of the squared transverse angular velocity is negative. In other words, spin about the *maximum* axis is *stable* in the case of energy dissipation. A spacecraft shaped like a frisbee, with energy dissipation, will tend to lose all nutation and spin perfectly about its maximum axis.

In the case that $I_3 < I_t$ (i.e. a spacecraft shaped like a pencil), the rate of change for the transverse angular velocity is *positive*. So, spinning about the minimum axis is *unstable* in the case of energy dissipation. A pencil-shaped spacecraft that is spinning stably about its maximum axis will begin to nutate more and more as energy gets dissipated. Eventually, it will fall over and start spinning about its maximum axis.

We can visualize this on the Poinsot ellipsoid. As energy is dissipated, the angular velocity spirals away from the minimum axis, falls over the intermediate axis, and spirals toward the maximum axis.

**Question**: What happens if we start *adding* energy to the system?

The Explorer 1 spacecraft is picture below, launched 1958. It was designed to spin about its minimum axis. Look at this spacecraft, what do you think happened and why?

In [15]:

```
Image("explorer.jpeg", width=500)
```

Out[15]: