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
getCurrentTime
2016-03-21 01:07:24.199926 UTC
In [2]:
:set -XFlexibleContexts

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

-- Scalar multiplication of a vector.
scale :: Num t => Vec t -> t -> Vec t
scale (Vec3 x y z) s = Vec3 (s*x) (s*y) (s*z)

scale (Vec3 1 2 3) 3

div :: Fractional t => Vec t -> t -> Vec t
div (Vec3 x y z) d = Vec3 (x/d) (y/d) (z/d)

Vec3 3.0 3.0 3.0 `div` 3.0

extend :: Num t => Vec t -> Vec t -> Vec t
extend (Vec3 x1 y1 z1) (Vec3 x2 y2 z2) = Vec3 (x1 + x2) (y1 + y2) (z1 + z2)

extend (Vec3 1.0 1.0 1.0) (Vec3 4.0 5.0 6.0)

-- Vector cross product. This can be verified using the right hand rule
-- and some examples.
cross :: Num t => Vec t -> Vec t -> Vec t
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 0 0) (Vec3 0 1 0)

-- Vector magnitude/absolute value.
magnitude :: (Floating t) => Vec t -> t
magnitude (Vec3 x y z) = sqrt $ x^2 + y^2 + z^2

magnitude (Vec3 1 0 0)
magnitude (Vec3 1 1 0)
magnitude (Vec3 1 1 1)
Vec3 3 6 9
Vec3 1.0 1.0 1.0
Vec3 5.0 6.0 7.0
Vec3 0 0 1
1.0
1.4142135623730951
1.7320508075688772

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 [4]:
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
Vec3 3.0 0.0 0.0
Vec3 15 0 0

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 [5]:
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 [6]:
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 [7]:
-- Mass in kg and distance in meters, result in Newtons
grav 1.989e30 5.972e24 149.6e9
3.541706104156539e22

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 [8]:
-- Some helpers for angles
degrees :: Floating x => x -> x
degrees x = x/pi*180

degreesV :: Floating x => Vec x -> Vec x
degreesV (Vec3 x y z) = Vec3 (degrees x) (degrees y) (degrees z)

radians :: Floating x => x -> x
radians x = x/180*pi

radiansV :: Floating x => Vec x -> Vec x
radiansV (Vec3 x y z) = Vec3 (radians x) (radians y) (radians z)

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 -- Solid cylinder with radius
                | 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) m = (m * r^2) / 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 (kg*m^2)
moment :: Num a => a -> a -> a
moment r m = r^2 * m

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

-- Angular momentum of a body with moment of inertia ia and angular velocity omega.
-- Omega is a vector quantity describing the rotation about the 3 principle axes in rad/s.
-- Units: kg*m^2/s^-1
angMomentum :: Num a => a -> Vec a -> Vec a
angMomentum ia omega = scale omega ia

-- Momentum of a 20kg hollow cylinder with radius 2m spinning at 10rad/s about the Z axis
angMomentum (momentV (HollowCyl 2.0) 20) (Vec3 0 0 10)
Vec3 0.0 (-7.5) 0.0
Vec3 0.0 0.0 800.0

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 [9]:
-- The angular momentum of a body after applying a torque for a period of time t.
torqueMoment :: Floating a => Vec a -> Vec a -> a -> Vec a
torqueMoment torque am t = extend am (scale torque t)

-- 3N force applied 5m from the center of a 20kg object
-- with radius 2m spinning at 10rad/s around the Y axis for 1s
cylTMoment = torqueMoment trq cylMomentum 1
    where
    trq = torque (Vec3 0.0 0.0 3) (Vec3 5.0 0.0 0.0)
    cylMomentum = angMomentum (momentV (HollowCyl 2.0) 20.0) $ Vec3 0.0 10.0 0.0

:t cylTMoment
cylTMoment
cylTMoment :: forall a. Floating a => Vec a
Vec3 0.0 785.0 0.0

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 zero forces and torques, these quantities will remain the same. We can use this to answer the Cassini-Huygens problem in the section 1.2 review.

Cassini

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

Solution

We know that the spacecraft is rotating 5600kg of body and 100kg of antenna with a rate of 6deg/s. Thus, the initial angular momentum is:

In [10]:
cassiniBodyIa = (momentV (SolidCyl 1.0) 5600) -- Main body of the spacecraft
cassiniRetAntIa = moment 1 50 -- Moment of inertia from a single arm when retracted
cassiniStartIa = cassiniBodyIa + (2 * cassiniRetAntIa)

cassiniL = angMomentum cassiniStartIa $ Vec3 (radians 6.0) 0.0 0.0

cassiniL
Vec3 303.68728984701335 0.0 0.0

After the arms are extended, we have two bodies: the main fuselage of the spacecraft, and two point masses extended 10m out. This changes our moment of inertia, and will change the rotational velocity as well; let's define a function to help us figure out what the resulting velocity will be.

In [11]:
angVelocity :: Floating a => a -> Vec a -> Vec a
angVelocity ia am = div am ia

cassiniExtAntIa = moment (1 + 10) 50.0
cassiniEndIa = cassiniBodyIa + (2 * cassiniExtAntIa)
cassiniEndVelocity = angVelocity cassiniEndIa cassiniL

degreesV cassiniEndVelocity -- Convert back to degrees for EdX
Vec3 1.167785234899329 0.0 0.0

So, approximately 1.17 deg/s.