Important Note

This post ended up being more elementary than I initially thought it would be, but the last post on this blog was from 2 years ago so I’m going to send it. Getting to a satisfying MVP for this project was a larger undertaking than I thought, so I’m shipping from a less satisfying place.

Treat this blog as a soft introduction to using the Julia DifferentialEquations package, and hopefully in a follow on I’ll get deeper into simulating space mechanics like I originally planned.

Moon Cannon: Can You Shoot Earth from the Moon?

With technological advancements and decreasing launch costs, establishing permanent manufacturing and production facilities in Low Earth Orbit and on the lunar surface is becoming inevitable. As we move closer to this reality, innovative methods for transporting cargo between the Moon and Earth are worth exploring. One such method of transportation, though it may sound far-fetched, is launching a payload from the Moon to Earth using a cannon. The low gravity and lack of atmosphere on the Moon present a unique environment that enables this possiblity. The cannon is really a mathematical simplification, there are real solutions like a Railgun or SpinLaunch are realistic way to impart a massive, instantaneous, initial velocity to something without any sort of combustive propulsion.

Setting Up the Problem

We’ll be using Julia to simulate the trajectory of a payload shot from the Moon towards Earth. I’ll explain the math and code we are setting up after the environment is setup. Here are the packages we need:

using DifferentialEquations
using LinearAlgebra
using Plots
1plotlyjs()

1: To get pretty, interactive plots

Here are some fundamental constants we’ll use throughout the simulation:

const G = 6.67430e-11  # Gravitational constant (m³/kg·s²)
const M_Earth = 5.97e24  # Mass of Earth (kg)
const M_Moon = 7.34767309e22  # Mass of Moon (kg)
const R_Earth = 6.371e6  # Radius of Earth (m)
const R_Moon = 1.737e6  # Radius of Moon (m)
const d_Earth_Moon = 384400e3  # Average distance between Earth and Moon (m)
const payload_mass = 1000.0  # Mass of payload (kg)

I tried using Unitful for tracking the units, but ran into issues and upon searching the Julia forums for an answer I found an open thread explaining why DifferentialEquations.jl doesn’t play nice with Unitful that was opened by ME! It’s crazy how I never learn from my mistakes…

Differential Equations Setup

The core of this problem lies in setting up and solving the equations of motion. We’ll start with a basic restricted three-body problem, meaining only accounting for the Moon’s and Earth’s gravity. This is a gross oversimplification for a few reasons. The first being that in this simulation the Moon and the Earth will remain static, introducing their movement is more complicated than I want to cover in this introduction to the problem. The other oversimplification is ignoring the gravitational pull from the Sun and Jupiter, and evenutally maybe the rest of the planets in the solar system. Over the time period of this simulation we could get away with combining these into one static force vector that gets applied every step, but I would rather wait until I can properly build out the simulation.

The equations of motion are:

\[ \begin{aligned} \frac{d\vec{r}}{dt} &= \vec{v}, \\[10pt] \frac{d\vec{v}}{dt} &= \vec{a}_{\text{Moon}} + \vec{a}_{\text{Earth}} + \vec{a}_{\text{drag}}. \end{aligned} \]

Ok so it is a little crazy to simulate such a simple environment, and then throw atmospheric drag into the equation. But hear me out, defining success criteria for this kind of problem can be difficult, since a lot of solutions get so close that they should be counted as wins, in this case drag creates a larger target that also accounts for velocity. If our spacecraft gets close enough that drag pulls it to the surface of Earth then we can call the simulation a success.

const payload_area = 10.0  # Cross-sectional area of payload (m²)
const rho_0 = 1.225  # Sea-level atmospheric density (kg/m³)
const H = 7000  # Scale height of atmosphere (m)
const Cd = 0.47  # Drag coefficient (dimensionless)

function atmospheric_drag(position, velocity)

    altitude = norm(position) - R_Earth
    if altitude > 100000
        return [0.0, 0.0]
    end

    rho = rho_0 * exp(-altitude / H)
    v_mag = norm(velocity)
    drag_magnitude = 0.5 * rho * v_mag^2 * Cd * payload_area
    drag_force = -drag_magnitude * velocity / v_mag

    return drag_force
end

Next, we define the differential equation function in Julia:

function f(du, u, p, t)
    # du: derivative of state vector
    # u: state vector (position and velocity)
    # p: parameters (unused in this function)
    # t: time (unused in this function)

    # Extract position and velocity from state vector
    r = u[1:2]  # Position vector (x, y)
    v = u[3:4]  # Velocity vector (vx, vy)

    # Define positions of Moon and Earth
    r_Moon = [0.0, 0.0]  # Moon at origin
    r_Earth = [d_Earth_Moon, 0.0]  # Earth along x-axis

    # Calculate relative positions
    r_rel_Moon = r - r_Moon  # Spacecraft position relative to Moon
    r_rel_Earth = r - r_Earth  # Spacecraft position relative to Earth

    # Calculate gravitational accelerations
    # Using Newton's law of gravitation: F = G * M1 * M2 / r^2
    # Acceleration = Force / Mass = G * M / r^2
    a_Moon = -G * M_Moon * r_rel_Moon / norm(r_rel_Moon)^3
    a_Earth = -G * M_Earth * r_rel_Earth / norm(r_rel_Earth)^3

    # Calculate atmospheric drag
    drag_force = atmospheric_drag(r_rel_Earth, v)
    a_drag = drag_force / payload_mass  # Convert force to acceleration

    # Set up differential equations
    du[1:2] = v  # dr/dt = v
    du[3:4] = a_Moon + a_Earth + a_drag  # dv/dt = sum of accelerations
end

I tried to comment the code as much as possible to help this make sense to someone who has never seen DifferentialEquations.jl. The state vector can be weird at first sight but should be similar to what you see in Matlab.

Detecting Collisions

Julia is pretty fast, but this is at minimum a multiple day simulation so it is important to define end conditions so that we don’t go on forever for really bad solutions (we crash back into the Moon) or really great solutions (we nail Earth head on). To handle this we just set callbacks that check our position and see if we are inside of either bodies radius.

function Earth_collision(u, t, integrator)
    r = u[1:2]
    r_Earth = [d_Earth_Moon, 0.0]
    return norm(r - r_Earth) - R_Earth
end

function Moon_collision(u, t, integrator)
    r = u[1:2]
    r_Moon = [0.0, 0.0]
    return norm(r - r_Moon) - R_Moon
end

Earth_cb = ContinuousCallback(Earth_collision, terminate!)
Moon_cb = ContinuousCallback(Moon_collision, terminate!)
cb = CallbackSet(Earth_cb, Moon_cb)

Run the Simulation

Now like any good optimistic software engineer, I’m going to assume that I’ll get things right first time and not invest any time into a helper function to allow me to quickly iterate on the initial conditions.

1v0 = [1500, 1500]
2r0 = v0 / norm(v0) * (R_Moon + 10)
3tspan = (0.0, 10.0 * 24 * 3600)
4prob = ODEProblem(f, [r0; v0], tspan)
sol = solve(prob, Tsit5(), callback=Earth_cb, reltol=1e-8, abstol=1e-8)

5x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]

1: Guesstimate initial velocity
2: Make the initial position the same direction as the initial velocity, and 10 meters above the Lunar surface
3: Guesstimate timespan. 10 days should be more than plenty
4: Setup the Differential Equation
5: Grab the x and y positions from the solution.

Now lets plot the results:

plot(x, y, aspect_ratio=:equal, label="cannonball", xlabel="x (m)", ylabel="y (m)")
scatter!([0], [0], label="Moon", markersize=5, color="grey")
scatter!([d_Earth_Moon], [0], label="Earth", markersize=10, color="cornflowerblue")

title!("First Simulation!")

Well we clearly didn’t make it far, lets try again without the Earth included in the plot so we can see what happened. (Or you could use plotly to zoom in)

plot(x, y, aspect_ratio=:equal, label="cannonball", xlabel="x (m)", ylabel="y (m)")
scatter!([0], [0], label="Moon", markersize=5, color="grey")

title!("First Simulation! (But Zoomed in!)")

Spirographs are cool, but not really what we are looking for. Lets turn all of this code into a function so that we can iterate on the initial conditions faster.

function plot_trajectory(v0, title)
    r0 = v0 / norm(v0) * (R_Moon + 10)
    tspan = (0.0, 20.0 * 24 * 3600)
    prob = ODEProblem(f, [r0; v0], tspan)
    sol = solve(prob, Tsit5(), callback=Earth_cb, reltol=1e-8, abstol=1e-8)

    x = [u[1] for u in sol.u]
    y = [u[2] for u in sol.u]

    flight_time_days = sol.t[end] / (24 * 3600)

    final_pos = sol.u[end][1:2]
    r_Earth = [d_Earth_Moon, 0.0]

    distance_to_Earth = norm(final_pos - r_Earth) - R_Earth

    plot(x, y, aspect_ratio=:equal, label="cannonball", xlabel="x (m)", ylabel="y (m)")
    scatter!([0], [0], label="Moon", markersize=5, color="grey")
    scatter!([d_Earth_Moon], [0], label="Earth", markersize=10, color="cornflowerblue")
    title!(title)
end

Try the same initial condition to make sure we wrapped everything up correctly:

plot_trajectory([1500, 1500], "First Simulation, plotted from a function")

Looks exactly the same. Now what if we go faster:

plot_trajectory([2000, 2000], "First Simulation, plotted from a function")

Ok, medium speed?

plot_trajectory([1700, 1700], "First Simulation, plotted from a function")

On the bright side the simulation passes the smell test for what we expect an orbital simulation to look like. The last attempt seems to be about the right magnitude, but what if we pointed the cannon right at the Earth?

plot_trajectory([2400, 0], "First Simulation, plotted from a function")

Well this problem is certainly easier than I had originally envisioned. What if I wanted to shoot from the dark side of the Moon.

plot_trajectory([-2300, 800], "First Simulation, plotted from a function")

Also super easy…

Conclusion

So I think we learned two things here, the Julia Differential Equations library makes getting problems like this going really really easy, and when you make massive oversimplifications for a problem it makes it look really easy. I think from this I have the confidence that I needed to invest time into a far more realistic simulation of the problem. For now this was a great way to dust off my Julia and diff eq part of my brain.

Next Steps

From here I think the next logical step is fully simulating all of the bodies in the solar system. Whether I do that with differential equations like this example, or use something like Horizons is an exercise for the next post. I think another cool element of this problem would be giving the cannonball/payload some small amount of thrust it can use to see if it can circularize itself in LEO, but I just really want an excuse to turn this into an optimization problem. There is also the question of sending payloads from the Moon to Mars which opens up a ton of options for different things to optimize for like transfer time or minimum initial velocity by doing something crazy like a Venus flyby. Another cool end condition for this series would be trying to impact a landing site on another planet, but I really don’t want to deal with that much aerodynamics so I doubt I’ll get that crazy with it.

Reuse

CC BY 4.0