Introduction

Julia's mathematical syntax and support for linear algebra constructs makes it an ideal language for scientific computing.

In this problem, we will investigate the space flight of the failed Apollo 13 mission, using Julia to simulate the system of differential equations that govern the motion of the spacecraft.

Problem and code adapted from the online course, Differential Equations in Action, taught by professor Jörn Loviscach

Download the code

Background

For those of you who are unfamiliar with the story of Apollo 13, it was the seventh manned mission to the moon, and took place during April, 1970. Their mission was to land on the moon and explore a region of the moon's surface called the Fra Mauro formation. Two days and 200,000 miles after launch, disaster struck.

After turning on the hydrogen and oxygen tank stirring fans, the crew heard a loud explosion, which was the sound of the number-2 oxygen tank blowing up. Damage to the service module would have made a safe return from a moon landing impossible, so Houston decided to abort the mission.

There were two options for safe return. The quickest route would have been a direct trajectory, essentially slowing down and falling back to earth. However, this required jettisoning the Lunar Module, which would have prevented the crew from landing safely. Instead NASA chose a circumlunar route, which would "slingshot" the astronauts around the moon, using the moon's gravity to send the spacecraft back towards earth.

A plot of Apollo 13's trajectory.

Prior to the oxygen tank explosion, the mission crew had just completed a burn that would slow them down so they could enter the moon's atmosphere and land on the surface. With the landing aborted, the crew had to initiate another burn that would speed them up and carry them around the back of the moon.

Our task, as NASA engineers, is to determine the amount of boost to give the spaceship so that it will return home safely.

Physics

Before we can start programming the simulation, we must first understand the basic physics that govern the system.

Let's first classify the problem in terms of Newtonian mechanics. There are three bodies in the simulation: the Earth, the Moon, and the Command Module. The moon is in orbit around the earth, and the command module is a projectile between the two. During its flight, the command module will experience forces due to gravity from both bodies, which will induce an acceleration. Let's examine Newton's law of gravitation (and his second law of motion) to see how this will effect the motion of the command module.

This describes the gravitational force between any two bodies, where and are their masses. represents the distance between them and represents the unit vector that points from on body to the other. is the gravitational constant. This is the general expression. For our problem, we can use known masses for the earth and the moon to describe the gravitational force exerted by either body on the spacecraft:

Where , , and are the masses of the earth, moon, and spacecraft, respectively, and and are the distances (with associated unit vectors) between the spacecraft and the earth and the moon.

The total force on the spacecraft at any instant is the sum of all forces acting on it, which we will assume to be the sum of the gravity from the earth and the moon:

We know from Newton's second law that the force experienced by a body is equivalent to the product of its mass and the acceleration on the body.

This allows us to express the acceleration experienced by the spacecraft in terms of known constants and one dependent variable, namely, the position of the spacecraft relative to the moon and the earth.

We will assume that the position of the earth is constant, at the origin. The moon, however, moves periodically around the earth in the following fashion (using cartesian coordinates):

Where is the radius of the moon's orbit, which we will treat as constant, is the period of the moon's orbit, and is the phase angle of the orbit, or the initial angular position of the earth relative to the positive x-axis. is in seconds.

We can now express a system of linear differential equations:

Constants

Let's get started. Open Julia Studio and create a new folder, called Apollo. Create a new file within the project folder called main.jl and then create another new file called constants.jl

It's convention to name constants with all uppercase letters in Julia. It's also convention to name variables with underscores between logical separation points. We precede the definition of a constant with the keyword const, which instructs the compiler that this variable is immutable. It also treats the constant as a global variable, without the performance hit you incur by using globals. We will define the constants in a module and export them. We can then use them by using the module name.

Let's enter the following constants:

# constants.jl

module constants

export ME, RE, G, MM, RM, MCM, DISTANCE_TO_MOON, MOON_PERIOD, MOON_INITIAL_ANGLE, ORIGIN
export TOTAL_DURATION, MARKER_TIME, TOLERANCE, INITIAL_POSITION, INITIAL_VELOCITY

const ME = 5.97e24 # Mass of earth in kg
const RE = 6.378e6 # Radius of earth in m (at equator)
const G = 6.67e-11 # Gravitational constant in m3 / kg s2
const MM = 7.35e22 # Mass of moon in kg
const RM = 1.74e6 # Radius of moon in m
const MCM = 5000. # Mass of command module in kg
const DISTANCE_TO_MOON = 400.5e6 # m (actually, not at all a constant)
const MOON_PERIOD = 27.3 * 24.0 * 3600. # s
const MOON_INITIAL_ANGLE = pi / 180. * -61. # Radians
const ORIGIN = [0., 0.] # Vector that represents the origin

const TOTAL_DURATION = 12. * 24. * 3600.# s
const MARKER_TIME = 0.5 * 3600. # (used in error correction) s
const TOLERANCE = 100000. # (used in error correction) m

const INITIAL_POSITION = [-6.701e6, 0.] # Initial vector for the spacecraft in m
const INITIAL_VELOCITY = [0., -10.818e3] # Initial velocity for spacecraft in m/s

end

In main.jl, load in these constants:

# main.jl

using constants

Julia's Type System

Let's take a moment to talk about types in Julia. For a detailed description, we recommend reading the material in the Julia Manual, but we'll address some of it here.

In addition to the concrete types that the language defines for you, Julia allows you to express composite types. You can think of a composite type as a struct in C, or an object, or a template, depending on what languages you are familiar with. The key difference between concrete types in Julia and objects in other languages is that concrete types do not have methods associated with them. This functionality can be achieved, as we will see later, but not by defining methods on types.

Let's start with an example. First, create a file called types.jl. This is where we will define all of our types, again in a module. This simulation involves three bodies, so it makes sense to create a Body type:

# types.jl

module types

type Body{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
end

Let's talk about what's going on here. On the first line, we use the keyword type, followed by Body{T}. This is an example of a parametric type, where we can pass in a type, e.g. Float64 as a parameter in the type instantiation, and that type will be applied to each of the named fields. Here, we specify that each body will have 4 such named fields: mass, which will inherit the type passed in, velocity, which will be a vector (a one 1xN array) with values of type, T, radius, which will be of type, T, and position, another vector whose values are of type, T. The :: operater is a type assertion, asserting that those fields must be of the specified type.

Here's how we might instantiate a new Body:

# With explicit type parameter
earth = Body{Float64}(ME, [0., 0.], RE, ORIGIN)

# Without explicit type parameter
earth = Body(ME, [0., 0.], RE, ORIGIN)

We can then reference each of its attributes using the normal dot notation: earth.mass, etc.

Let's define a few more types:

# types.jl

module types

export Body, Moon, Command_Module, EarthMoonSystem

type Body{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
end

typealias Moon Body # Create a new type, Moon, that has all the
                    # attributes of Body, but is treated as its own type

type Command_Module{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
end

type EarthMoonSystem # This is a non-parametric type, since it doesn't take any arguments
    time::Float64
    earth::Body # See what we did there?
    moon::Moon # Now we can do type assertion with our composite types
    command_module::Command_Module
end

end

In these lines, we use the typealias keyword to create two new types, one for the moon, and one for the command module, that have the same named attributes and parameters as Body, but have different names. We will see why this is important later. Use this module in main.jl

# main.jl

using constants
using types

Building a simulation

Let's instantiate some new bodies in main.jl, and create a function, simulate, that will be our interface for running the simulation.

# main.jl

using constants
using types

# Initialization of our bodies
earth = Body(ME, [0.0, 0.0], RE, ORIGIN)
moon = Moon(MM, [0., 0.], RM, moon_position(0.0))
command_module = Command_Module(MCM, INITIAL_VELOCITY, 5.0, INITIAL_POSITION)
world = EarthMoonSystem(0.0, earth, moon, command_module)

function simulate()
    current_time = 1.

    while current_time <= TOTAL_DURATION
        update(world, current_time)

        current_time += 1.
    end
end

In this function, we initialize the simulation time at 1, and update the world once a second until we reach the end of the simulation.

Now we must define the update method. Create a new file called system.jl which will contain all the functions for instances of EarthMoonSystem types. What do we want this function to do? We want it to update the status of all the bodies within it.

# system.jl

function update(me::EarthMoonSystem, time::Float64)
    # `me` is a reference to this instance of EarthMoonSystem

    me.time = time # Explicitly set the time on the system

    update(me.moon, time) # Update the moon
    update(me.command_module, time) # Update the command_module

    return me # Return the system
end

Wait a second. Are we calling update again? Does that mean this is a recursive function? Shouldn't we name those functions something like update_moon and update_command_module? No, because of Julia's adherence to a pattern called multiple dispatch.

Multiple dispatch allows us to define functions as collections of methods, where each method is defined for a specific combination of types. Julia calls whichever function method best applies to the types of the arguments. Here is an example. Let's say we define a function, call_person, that calls a specific person:

type Person
    name::String
end

function call_person(person)
    print("Hello $person.name")
end

mom = Person("mom")

call_person(mom)

In this example, we define a Person type, and a function for calling a person, which says hello to that person. Now let's take it another step further to illustrate how mutliple dispatch works:

type Person
    name::String
end

typealias Mom Person
typealias Boss Person
typealias Doctor Person

function call_person(person)
    print("Hello $person.name")
end

function call_person(person::Mom)
    print("I love you, mom")
end

function call_person(person::Boss)
    print("Sorry I'm late, $person.name")
end

function call_person(person::Doctor)
    print("What's up, doc?")
end

doc = Doctor("Rich")

call_person(doc)

Now we've extended the call_person function with three more methods. The correct version of call_person will be called, depending on the type of arguments. If the type of the arguments do not match any of those specified, the first method will be called, since no type assertions are used there. In this case, we create a new Doctor and then call call_person. The output of this code should be, "What's up, Doc?".

Now that it's clear why we can have three different update methods, let's write them. Create two new files, command-module.jl and moon.jl.

In each, add the following method:

# command-module.jl
function update(me::Command_Module, time::Float64)
    return me
end

# moon.jl
function update(me::Moon, time::Float64)
    return me
end

Now let's focus on the moon. In moon.jl, enter the following.

# moon.jl
function moon_position(time::Float64)
    moon_angle = MOON_INITIAL_ANGLE + 2.0pi * time / MOON_PERIOD
    x::Float64 = DISTANCE_TO_MOON * cos(moon_angle)
    y::Float64 = DISTANCE_TO_MOON * sin(moon_angle)

    return [x, y]
end

function update(me::Moon, time)
    me.position = moon_position(time)

    me
end

Now the update method calls the moon_position function, which returns the x and y coordinates of the moon as a vector.

Now, in command-module.jl:

# command-module.jl

function acceleration(time::Float64, pos::Vector{Float64})
    moon_pos = moon_position(time) # Get the moon's position at this time

    distance_from_earth = pos # Earth is at the origin
    distance_to_moon = pos - moon_pos # Distance between command module and moon
    mag_e = norm(distance_from_earth) # Get the magnitude of this vector
    mag_m = norm(distance_to_moon)
    return -G * (ME * distance_from_earth / mag_e^3 + MM * distance_to_moon / mag_m^3)
end

function update(me::Command_Module, time::Float64)
    a = acceleration(time, me.position) # Calculate the acceleration on the module

    # Using Euler's method
    me.position += me.velocity # Increment position vector by velocity vector
    me.velocity += a # Add acceleration vector to velocity vector

    me
end

This code use's Euler's method to solve the system of differential equations. First, we calculate the acceleration on the command module based on its distance from the earth and the moon. Since velocity the the rate of change in position, we can add it to current position to get a new position. Since acceleration is the rate of change of velocity, we add it to velocity for use in the next step. If we use small enough steps, this numerical solution approximates the analytical solution to the system of equations. Notice how easy it is to add vectors in Julia. All normal mathematical operations work on arrays.

If we include these files in main.jl, and add a call to our simulate() function, we should have a running simulation!

Viewing the Results

We want to see the results, right? Let's create some variables to hold the results. Somewhere in the function body of the simulate function, before the start of the loop, create a new vector: position_list = Vector{Float64}[]. This variable is a vector that will hold all the position vectors. Now, in the body of the loop, after calling update, add the following line: push!(position_list, copy(world.command_module.position)).

Now, after the conclusion of the loop, return the value of position_list. We will use Julia's writecsv method to output our results to a csv file:

# in main.jl, after simulate() function

@time pos = simulate() # get the results, show elapsed time
writecsv("output.csv", pos) # write the output to a csv file

The entire file should look like this:

# main.jl

using constants
using types
include("physics.jl")
include("moon.jl")
include("command-module.jl")
include("system.jl")

# Initialization of our bodies
earth = Body(ME, [0.0, 0.0], RE, ORIGIN)
moon = Moon(MM, [0., 0.], RM, moon_position(0.0))
command_module = Command_Module(MCM, INITIAL_VELOCITY, 5.0, INITIAL_POSITION)
world = EarthMoonSystem(0.0, earth, moon, command_module)

function simulate()
    position_list = Vector{Float64}[] # m
    current_time = 1.

    while current_time <= TOTAL_DURATION
        update(world, current_time)

        push!(position_list, copy(world.command_module.position))

        current_time += 1
    end

    return position_list
end

@time pos = simulate()
writecsv("output.csv", pos)

If we run this, notice that it takes a long time and generates a lot of data. If we graph the data in Excel the output looks wrong, because Euler's method leads to a large amount of error in the calculations. If we substitute another method called Heun's method, and use an adaptive step size instead of 1, then we can greatly reduce the error in our calculations. Here's how we can implement it.

First, change update for command_module to look like this:

# in command_module.jl

function update(me::Command_Module, time::Float64, h::Float64)
    acceleration0 = acceleration(time, me.position) # Get acceleration at current time
    velocityE = me.velocity + h * acceleration0 # Calculate Euler's velocity,
                                                # Accounting for the adaptive step size
    positionE = me.position + h * me.velocity   # Calculater Euler's position

    # Heun's method looks ahead to the next step and averages the two values
    velocityH = me.velocity + h * 0.5 * (acceleration0 + acceleration(time + h, positionE))
    positionH = me.position + h * 0.5 * (me.velocity + velocityH)

    # We will use the Heun's values
    me.velocity = velocityH
    me.position = positionH

    # Store these on the command_module so the simulation can access them
    me.positionH = positionH
    me.velocityH = velocityH
    me.positionE = positionE
    me.velocityE= velocityE

    me
end

We now have to modify the Command_Module type to include these new named fields:

# in types.jl

type Command_Module{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
    positionE::Vector{T}
    positionH::Vector{T}
    velocityE::Vector{T}
    velocityH::Vector{T}
end

And we modify the initialization code in main.jl:

# in main.jl

command_module = Command_Module(MCM, INITIAL_VELOCITY, 5.0, INITIAL_POSITION, INITIAL_POSITION, INITIAL_POSITION, INITIAL_VELOCITY, INITIAL_VELOCITY)

You'll also notice that the update method in command_module.jl now takes an additional argument, h, which we need to pass in. This is the step size to use. The simulation will adjust the step size depending on the error we generate between Heun's values and Euler's values, hence the name adaptive step size. Smaller values of h give more accurate results, but take longer to run. Larger values generate less accurate results, but in fewer steps. We have arbitrarily set a TOLERANCE in constants.jl to 100,000 m, but you can experiment with different values.

Because we've modified the update method in command-module.jl, make sure to modify each update call to include h as the third argument. This happens twice: in main.jl: update(world, current_time, h), and also in system.jl: update(me.command_module, time, h).

Now, back in main.jl, let's initialize step size with a value of 0.1 right before the start of the loop. Also create a variable called h_new to store the adapted step size.

# in main.jl, in simulate() function

position_list = Vector{Float64}[]
current_time = 1.
h = 0.1 # s, set as initial step size right now but will store current step size
h_new = h # s, will store the adaptive step size of the next step

In place of current_time += 1, we will insert the following:

# in main.jl, in simulate() function

positionE = world.command_module.positionE
positionH = world.command_module.positionH
velocityE = world.command_module.velocityE
velocityH = world.command_module.velocityH

error_amt = norm(positionE - positionH) + TOTAL_DURATION * norm(velocityE - velocityH)
h_new = min(0.5 * MARKER_TIME, max(0.1, h * sqrt(TOLERANCE / error_amt))) # Restrict step size to reasonable range

current_time += h
h = h_new

Let's examine more closely what we're doing. When we calculate error_amt, we are getting a normalized error vector in our current position and adding in our long term error, which is our current normalized error vector in position, multiplied by the total duration of the simulation. This vector explains how much error we expect in our calculations given the current step size.

The following line picks a new step size, h_new, based on the ratio of our tolerance to our projected error.

We then increment current time by step size, and set the step size for the next iteration to be our new step size. Here's a snapshot of what main.jl should look like now:

# main.jl

using constants
using types
include("physics.jl")
include("moon.jl")
include("command-module.jl")
include("system.jl")

# initialization of our bodies
earth = Body(ME, [0.0, 0.0], RE, ORIGIN)
moon = Moon(MM, [0., 0.], RM, moon_position(0.0))
command_module = Command_Module(MCM, INITIAL_VELOCITY, 5.0, INITIAL_POSITION, INITIAL_POSITION, INITIAL_POSITION, INITIAL_VELOCITY, INITIAL_VELOCITY)
world = EarthMoonSystem(0.0, earth, moon, command_module)

function simulate()
    position_list = Vector{Float64}[]
    current_time = 1.
    h = 0.1 # s, set as initial step size right now but will store current step size
    h_new = h # s, will store the adaptive step size of the next step


    while current_time <= TOTAL_DURATION
        update(world, current_time, h)

        positionE = world.command_module.positionE
        positionH = world.command_module.positionH
        velocityE = world.command_module.velocityE
        velocityH = world.command_module.velocityH

        error_amt = norm(positionE - positionH) + TOTAL_DURATION * norm(velocityE - velocityH)
        h_new = min(0.5 * MARKER_TIME, max(0.1, h * sqrt(TOLERANCE / error_amt))) # Restrict step size to reasonable range

        current_time += h
        h = h_new

        push!(position_list, copy(world.command_module.position))
    end

    return position_list
end

@time pos = simulate()
writecsv("output.csv", pos)

Save Us!

So now, we can finally answer the question about how to save the astronauts. If we run the program and plot the output, we should see that the astronauts will return to the exact position they started at, right on earth:

A plot of Apollo 13's simulated trajectory when no corrective burns are fired.

In actuality, the astronauts had already fired one burn, called the mcc2 burn, which reduced their speed by about 7.04 meters/second. This is just enough to encourage a slow descent to the Lunar surface.

We can program that into the model by checking the time. If the time is correct, reduce velocity and set a flag to prevent the burn from firing again. In main.jl, in the simulate() method right after the update call, add this:

# in main.jl, in simulate() function

if !mcc2_burn_done && current_time >= 101104
    println("mcc2_burn fired")
    world.command_module.velocity -= 7.04 / norm(world.command_module.velocity) * world.command_module.velocity

    mcc2_burn_done = true
end

Make sure to initialize this variable with a value of false.

Here is what the trajectory would look like if the astronauts had fired the mcc2 burn and done nothing else (assuming they did not crash into the moon!):

A plot of Apollo 13's simulated trajectory when just the mcc2 burn is fired.

Now, our challenge is to fire a burn just large enough to save the astronauts from crashing into the moon, but small enough to prevent the astronauts from careening off into space forever. Add some code for a second burn (remember to initialize it as false). Also initialize a value for boost. Set it to 0 initially:

# in main.jl, in simulate() function

if !dps1_burn_done && current_time >= 212100
    println("corrective_burn fired")
    world.command_module.velocity += boost / norm(world.command_module.velocity) * world.command_module.velocity

    dps1_burn_done = true
end

Now experiment with different values of boost, seeing how it affects the trajectory of the spacecraft.

Conclusion

Here are some comparisons of different values of boost

boost = 10
boost = 100 (ouch)
boost = -10 (yikes)

As you can see here, a boost of about 10 m/s seems about right.