How to perform a soft landing on Mars

So, let’s say that you’re a spacecraft. Or in one. And you just arrived at the Martian Sphere of Influence, with a 15000km circular orbit. How do you land safely? You have no parachute, and a carton of eggs that must not break (a friendly reminder that there are probably no hens on Mars). A good idea would perhaps be to lower your periapsis enough to touch the Martian atmosphere, and utilize the resistance of the air to lower your apoapsis without spending any fuel.
How much fuel will I need? How long will it take? Are you sure that my eggs won’t break?
All these questions can be answered by formulating a model for the system and performing a simulation!
We start with the Martian atmosphere. A quick duckduckgoing (equivalent to googling but less invasive) points us to this NASA model, which we happily implement, only to plot it and realize that there is an error in the model, and the temperature goes below -273.15°C at an altitude of around 110km (which is probably not a very good idea). Disclaimer: this is an outdated model based on some quite old curve fits, and NASA provides a multitude of other more realistic and accurate options but since (thankfully) nobody assigned me the task of safely landing some eggs on Mars this model will be deemed sufficient. The side effect of the below absolute zero temperature is that the density asymptotically spikes on both sides around this point:
Atmospheric model error
Luckily, another duckduckgoing indicates that the folks at the University of Colorado have already devised a solution that does not include implementing a more advanced model, a polynomial curve fitting (the document disappeared as I was writing this post, please use the wayback machine). After patching up our atmospheric function, we move on to the spacecraft.
There is no upper limit for the complexity of the spacecraft dynamics model, and the lower limit is dictated by one’s conscience. So let’s just consider that the spacecraft is an infinitely small point in 2 dimensional space with fixed aerodynamic characteristics! Just like real life.
By conveniently placing Mars in the origin of the Cartesian coordinate system (0,0), the distance of the spacecraft from the center of Mars (yes, it is a perfect sphere with uniformly distributed mass in case you were wondering) becomes:

\(r = \sqrt{x^2+y^2}\\\)

where x,y are the coordinates of the spacecraft

In similar fashion we calculate the absolute velocity:

\(v = \sqrt{v_x^2+v_y^2}\\\)

where \(v_x,v_y\) are the velocity components of the spacecraft

The magnitude of the force of gravity is calculated using Newton’s law of universal gravitation:

\(F_{grav} = \mu \frac{m}{r^2}\\\)

where \(\mu\) is the standard gravitational parameter for Mars (pick the right units!) and m the mass of the spacecraft

Since we have the speed of the spacecraft and its altitude (remember to subtract Mars’ radius from r) to plug it in the atmosphere model, we can calculate the magnitudes of the aerodynamic forces:

F_D =\frac{1}{2} \rho v^2 C_D A \\
F_L =\frac{1}{2} \rho v^2 C_L A \\
C_L = \frac{L}{D} C_D \\

where \(F_D\) is the magnitude of the drag force, \(F_L\) the magnitude of the lift force, \(\rho\) the density derived by the atmospheric model, \(C_D\) the drag coefficient, \(C_L\) the lift coefficient, \(A\) the aerodynamically active area, \(v\) the magnitude of the velocity and \(\frac{L}{D}\) the lift to drag ratio

And yes, all the aerodynamic characteristics are constant since we’re aiming for maximum realism, we’ve been over that already!
I assume that you are familiar with drag and lift, \(\frac{L}{D}\) is just a design parameter specifying the relationship between the two in a capsule. Together with the capsule’s ballistic coefficient, they are the two major factors that dictate the trajectory of the spacecraft when aerobraking, or in general whenever it is flying (or falling) through the atmosphere.
In order to not break our precious eggs, we of course need to strap a huge ass rocket motor to the capsule to decelerate.
Every time the engine is producing thrust, it does so by spewing massive quantities of exhaust gasses with furious anger from the nozzle. To do so, it has to expend significant amounts of fuel and oxidizer (in the case of liquid rocket motors), which of course means that the spacecraft’s mass is dropping. The rate of change of the mass can be calculated:

\(\dot{m} = \frac{F_{thrust}}{g_0 I_{sp}}\\\)

where \(\dot{m}\) the rate of change of the spacecraft mass, \(F_{thrust}\) the magnitude of the thrust, \(g_0\) is the standard gravity and \(I_{sp}\) the specific impulse of the motor

Specific impulse is just a measure of how efficiently the rocket uses up its propellant, which in reality directly relates to the velocity of the exhaust gasses. Higher exhaust gas velocities produce more momentum with the same mass (amount of propellant used) which means more furious anger with less weight loss.

Now that we have the magnitudes of everything, we need to add the directions. Do do so we can just normalize the vector components to get the unit vector that’s pointing in the right direction. For example with velocity:

v_{1x} = \frac{v_x} {v} \\
v_{1y} = \frac{v_y} {v} \\\)

where the resulting \(v_{1x}, v_{1y}\) are the components of the normalized velocity vector

This vector will be used to set the direction of most of the forces, namely the drag, lift and motor thrust. The drag and motor thrust act in the opposite direction of this vector. Drag does that by definition, and we pick that for the thrust because we want to decelerate (the term for that is that we are burning in the retrograde direction). For the lift, there are two cases. Since we considered the aerodynamic characteristics constant, it can either point to the sky (lift-up) or to the ground (lift-down). A configuration with downward lift results in more abrupt deceleration and higher forces, as the spacecraft dives faster in the thicker layers of the lower atmosphere. On the other hand, when the lift vector is pointing upwards the vertical descent of the spacecraft is slower and it is subjected to less stress, with the disadvantage of a longer aerobraking duration. To minimize the stress on our zygotes, we pick the lift-up configuration and calculate the direction as a right angle to the velocity vector:

l_{1x} = v_{1x}cos(\theta_{L})-v_{1y}sin(\theta_{L}) \\
l_{1y} = v_{1x}sin(\theta_{L})+v_{1y}cos(\theta_{L}) \\\)

where \(\theta_{L}\) is the angle of the lift vector relative to the velocity

Finally, Newton’s second law of motion can be used to calculate the acceleration of the spacecraft from the acting forces, leading to the complete differential equation system. As mentioned, the thrust and drag are always applied in the retrograde direction (hence the use of the velocity unit vector), and the gravitational attraction to the direction of the origin (position unit vector).

a_x = \frac{\frac{x}{r} F_{grav} – v_{1x} F_D – v_{1x} F_{thrust} + l_{1x} F_L}{m} \\
a_y = \frac{\frac{y}{r} F_{grav} – v_{1y} F_D – v_{1y} F_{thrust} + l_{1y} F_L}{m} \\

And the complete system:

\dot{y} = \left[ \begin{array}{c} v_x \\ v_y \\ a_x \\ a_y \\ \dot{m} \end{array} \right], y_0 = \left[ \begin{array}{c} x_0 \\ y_0 \\ v_{x0} \\ v_{y0} \\ m_0 \end{array} \right], \quad
\dot{y}(t) = f(t,y(t),p), \quad y(t_0) = y_0 \\

where \(y_0\) are the initial states for each variable and p is the parameters vector (coefficients and the such)

To calculate the initial states, we place the spacecraft at \(x_0 = 15000, y_0 = 0\) (just a point in the desired 15000km circular orbit, on the x axis for simplicity) and calculate the initial velocity using

v_{apo0} = \sqrt{\mu_{Mars} \frac{2}{R_{A0}} – \frac{1}{a}} \\

where \(R_{A0}\) is the apoapsis radius (15000km), and a the semi-major axis of the orbit (again 15000km because it is circular)

To control when the motor is active or not, a simple if checks if the current time is within the pre-specified burn intervals.

After implementing everything and putting them together, the system is integrated using MATLAB’s (I’m sorry open source ecosystem, I truly am) built in explicit Dormand-Prince 5-4 (ode45). While dynamic orbital systems are already stiff and the addition of thrust intervals with abrupt changes make the system even stiffer, the simulation duration is small enough to allow a very small tolerance (hence time step) to ensure sufficiently correct integration.

After educated guesses on the spacecraft’s characteristics based on a scaled SpaceX crew dragon capsule and sufficient fiddling around with the manipulated variables, we get the following output (of course some prettifying stuff are added such as highlight the powered parts of the trajectory and a gradient of the logarithm of the atmospheric density):

Mars landing simulation
Mars landing simulation

I played enough with the numbers to achieve a touchdown velocity of 0.295 m/s (THE EGGS DID NOT BREAK) with a fuel mass fraction of 0.718. This means that 72% of the initial mass is fuel that got burnt in the process.

Next step? Strap an optimizer (most likely fmincon) to find the optimal solution with the smallest mass fraction. Because there are discrete phases in the simulation (thus discontinuities) this is not as simple as I originally anticipated.

That’s it! Thank you for reading!