Introduction to Space Operations

This iHaskell notebook contains notes on the EdX course "Introduction to Space Operations", taught by Claude Nicollier.

Last Run

In [1]:
import Data.Time
2016-03-21 02:46:31.843515 UTC

A note on units

Throughout this notebook I will use SI units. I will add annotations from the units Haskell package once I figure out how to use it.

In [2]:
-- TypeFamilies is needed when defining Unit instances
:set -XTypeFamilies

-- TypeOperators is needed for defining the EnergyFlux type
:set -XTypeOperators

-- DataKinds is needed for defining AngleDim and Radians
:set -XDataKinds

import Data.Metrology
import Data.Metrology.SI
import Data.Metrology.Show

Unit 1: Review of Mechanics

Unit 1 contains a brief review of the laws of mechanics as formulated by Newton. These are relatively easy to encode, so let's begin with some basic ones.

First, let's define some helper functions for working with vectors.

In [3]:
-- We'll primarily be working in 3-space. For problems with a rotational
-- component, we'll use two vectors: one for linear motion and one for
-- angular motion.
data Vec a = Vec3 a a a deriving Show

unitV :: Vec Length
unitV = Vec3 (1.0 % Meter) (1.0 % Meter) (1.0 % Meter)

:t unitV

-- Scalar multiplication of a vector.
scale :: (Num n) => Vec (Qu d l n) -> n -> Vec (Qu d l n)
scale (Vec3 x y z) s = Vec3 (s *| x) (s *| y) (s *| z)

scale (Vec3 (1 % Meter) (2 % Meter) (3 % Meter)) 3
unitV :: Vec (Qu '['F Length One] 'DefaultLCSU Double)
Vec3 3.0 m 6.0 m 9.0 m
In [4]:
div :: (Fractional n) => Vec (Qu d l n) -> n -> Vec (Qu d l n)
div (Vec3 x y z) d = Vec3 (x |/ d) (y |/ d) (z |/ d)

Vec3 (3.0 % Meter) (3.0 % Meter) (3.0 % Meter) `div` 3.0
Vec3 1.0 m 1.0 m 1.0 m
In [5]:
-- Vector cross product. This can be verified using the right hand rule
-- and some examples.
cross (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 i j k where
    i = (y1 |*| z2) |-| (z1 |*| y2)
    j = (x1 |*| z2) |-| (z1 |*| x2)
    k = (x1 |*| y2) |-| (y1 |*| x2)

cross (Vec3 (1 % Meter) (0 % Meter) (0 % Meter)) (Vec3 (0 % Meter) (1 % Meter) (0 % Meter))
Vec3 0.0 m^2 0.0 m^2 1.0 m^2
In [ ]:
answer1 = (2.0 % Meter) |^ sTwo

-- Vector magnitude/absolute value.
--magnitude :: Vec (Qu d l n) -> Qu d l n
magnitude (Vec3 x y z) = qSqrt ((qSq x) |+| (qSq y) |+| (qSq z))

:t magnitude

--magnitude $ Vec3 (1 % Meter) (0 % Meter) (0 % Meter)
--magnitude $ Vec3 (1 % Meter) (1 % Meter) (0 % Meter)
--magnitude $ Vec3 (1 % Meter) (1 % Meter) (1 % Meter)

<<< This is where units ends >>>

Imbalanced forces acting upon an object will cause it to accelerate. The direction and the magnitude of this acceleration is given by the vector sum of the set of forces acting on the object, which we can represent by:

$$\vec F = m\vec a$$

In Haskell, we can trivially calculate the magnitude of the acceleration imparted by a given force or the magnitude of a force required for a given acceleration using:

In [ ]:
accel :: (Fractional t) => Vec t -> t -> Vec t
accel f m = div f m

accel (Vec3 15.0 0.0 0.0) 5.0 -- 15N force applied to a 5kg object, result in m/s

force :: (Num t) => Vec t -> t -> Vec t
force a m = scale a m

force (Vec3 5 0 0) 3 -- 3kg mass to 5m/s^2

The generalization of this law when the mass of the body varies over time is given by:

$$ \vec F = \frac{d\vec p}{dt} $$

where the momentum $p$ is given by:

$$ \vec p = m \vec v $$
In [ ]:
momentum :: (Num t) => Vec t -> t -> Vec t
momentum v m = scale v m

The magnitude of the gravitational forces between two bodies is given by the equality:

$$ F = G \frac{Mm}{r^2}$$

Where $G$ is $6.673 \times 10^-11 \; \mathrm{m^3/kg/s^2}$, M is the mass of the larger body, m is the mass of the smaller body, and r is the distance between them.

In [ ]:
g = 6.673e-11
grav ma mb r = g * (ma * mb / r^2)

So, plugging in the values for the mass of the sun and earth and the distance between them we get:

In [ ]:
-- Mass in kg and distance in meters, result in Newtons
grav 1.989e30 5.972e24 149.6e9

For rotations, we consider a torque T which is generated by a force F acting on a body such that the line of action doesn't pass through the center of mass of the body or a fixed point. This force will tend to rotate the body around its CoM or fixed axis, changing its angular momentum L.

$$\vec T = \vec r \times \vec F$$$$\vec L = \vec r \times \vec p$$$$\vec T = \frac{d \vec{L}}{dt}$$

Torque vectors point along the axis of rotation, with the direction indicating the spin direction: when viewing the object along its axis in the direction of the torque vector, the object will rotate clockwise.

When we have a fixed axis of rotation, the angular momentum $L$ is related to the moment of inertia $I_a$ and the angular velocity $\omega$:

$$\vec L = I_a \vec \omega$$

Moment of inertia $I_a$ is given by the equation

$$ I_a = mr^2 $$

where $m$ is the mass of the body and $r^2$ is its distance from the axis of rotation.

If we have a composite body made of $i$ pieces each having a mass $m_i$ with its center of mass located a distance $r_i$ from the axis of rotation, then the moment if inertia is the sum:

$$I_a = \sum_i m_i r_i^2$$

Otherwise, for a continuous body, we apply the original equation continuously over the point masses comprising the body and sum them over the entire volume:

$$ I_a = \int\int_0^M r^2 dm $$

Note: In the videos they use this form of the equation:

$$I_a = \int\int\int_V r^2 \rho(r)dV$$

This is a more generalized form which accounts for varying density across the body. We'll be ignoring that factor and going with simpler forms. Specifically, we'll define the moments of inertia of some common volumes and use them to compose more complicated scenarios.

In [ ]:
torque :: (Num t) => Vec t -> Vec t -> Vec t
torque r f = r `cross` f

-- Force of 5N along the x axis acting on a point 1.5m from the center on the z axis.
-- Result in Nm, as a torque vector.
torque (Vec3 0.0 0.0 1.5) (Vec3 5.0 0.0 0.0)

data Volume a = SolidCyl a a -- Solid cylinder with radius and length
                | HollowCyl a -- Hollow cylinder with radius
                | SolidSphere a -- Solid sphere with radius
                | HollowSphere a -- Hollow sphere
                deriving Show

-- Moment of inertia of a volume along its axis of symmetry, in kg*m^2
momentV :: (Fractional a) => Volume a -> a -> a
momentV (SolidCyl r l) m = (m * r^2 * l) / 2
momentV (HollowCyl r) m = m * r^2
momentV (SolidSphere r) m = (2 * m * r^2) / 5
momentV (HollowSphere r) m = (2 * m * r^2) / 3

-- Moment of a point mass
moment :: (Num a) => a -> a -> a
moment r m = r^2 * m

-- Moment of inertia of a set of bodies
compMoment :: (Num a) => [(a, a)] -> a
compMoment bodies = foldr (+) 0 $ mmnts bodies
    where mmnts bodies = map (\x -> moment (fst x) (snd x)) bodies

angMomentum m r s = m * r * s

angMomentum 20 2 10 -- Momentum of a 20kg object with radius 2m spinning at 10deg/s

A torque applied to a rotating object will change its momentum. If it is applied in the direction of the rotation it will increase the rotational speed; otherwise it will decrease it.

The right hand rule can be used to determine the sign of an applied torque.

In [ ]:
torqueMoment t m s = m + (t * s)

-- 3N force applied 5m from the center of a 20kg object
-- with radius 2m spinning at 10deg/s for 1s
angMomentum 20 2 10
torqueMoment (torque 3 5) (angMomentum 20 2 10) 1

The two quantities p (linear momentum) and L (angular momentum) are conserved in a system - that is to say, given a balanced force diagram with net 0 forces, these quantities will remain the same. We can use this to answer the Cassini-Huygens problem in the section 1.2 review.


The Cassini-Huygens spacecraft was launched in 1997 towards Saturn. It was roughly of cylindrical shape with two straight arms along the cylinder whose purpose was to deploy antennas (considered as point mass) far from the main body some time after launch.

After deployment from the launcher, the spacecraft was turning along its main axis at a rate of 1 rpm or 6 degrees per second. After deployment of the antennas 90 degrees off the axis of the spacecraft, at the end of the arms supposed massless, what was the rotation rate of the spacecraft?

Data : Cassini total mass $M_c = 5600kg$ without the antennas, height $h_c = 5m$, diameter $d_c = 2m$. Mass of each antenna $m_a = 50kg$, length of each arm $l = 10m$.


We know that the spacecraft is rotating 5700kg (remember to add the mass of the two antennas) with a rate of 6deg/s. Thus, the initial angular momentum is:

In [ ]:
cassiniL = angMomentum 5700 1 6


After the arms are extended, we have two bodies: the main fuselage of the spacecraft, and two point masses extended 10m out.

In [ ]: