using DifferentialEquations
using LinearAlgebra
using Plots
1plotlyjs()
- 1
- To get pretty, interactive plots
Can you shoot Earth from the Moon? Simulate lunar ballistic launches with Julia. Explore orbital mechanics, gravity, and drag in a simplified Moon-to-Earth cannon model.
September 25, 2024
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.
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.
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()
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…
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.
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)
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]
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: