Simulating orbits, for a trans-lunar injection

July 16, 2026

Cyprien Hachem

no generative AI written text (read more ->)

Link to the trans-lunar injection simulator ->

Intro

This post aims to be a short introduction on how one might simulate an orbit.

In classical mechanics, Newton's law of universal gravitation governs orbits. It describes the force that two bodies exert on each other:

F = G\frac{m_1 m_2}{r^2}

$G = 6.6743 \times 10^{-11} \, \mathrm{m^3 \, kg^{-1} \, s^{-2}}$ is the gravitational constant
$m_1$ and $m_2$ are the masses of two bodies
$r$ is the distance between the two bodies

The law itself is very simple. In the two bodies case, the equation for the trajectories of motion can be worked out analytically, which means that given initial conditions (i.e. the positions and velocities of the two bodies), the complete shape of the trajectory can be determined; the position or velocity of one of the body can be calculated accurately at any point in time.

The problem goes from very easy to incredibly complicated when you introduce a third body. This is the well known three-body problem. At this point, in the general case, no closed-form solution exists, meaning that it is not possible to find some kind of formula which says where and when the bodies will be at any given point in time.

At this point, numerical simulation of orbits becomes the most practical method. But there are many design choices to make when numerically computing orbits, the rabbit hole of complexity goes deep depending on the level of precision desired.

The n-body problem

Given the fact that the general n-body problem does not have any closed-form solutions, there are two main ways to approach the problem:

On one hand, in special cases, the problem simplifies. For example, the restricted n-body problem considers that the third body is massless.

For the general case, we must resort to numerical integration, i.e. simulate numerically the trajectories of the n bodies. This is the approach used in the simulator ->. It involves 2 parts: modelling the system, and integration. It is especially advantageous because we can represent a lot of effects in the simulation, simply by adding terms to our model. The integration part then "just takes care of it".

Modelling the n-body system Earth-Moon-Spacecraft

The simplest way to model the Earth-Moon-Spacecraft system is by considering all the perturbations which affect the spacecraft's motion. In other words, we sum all the forces by the different bodies which act on the spacecraft, as well as any additional force we wish to model.

More specifically, in our case, there are 3 celestial bodies which act upon the spacecraft: the Earth, the Moon, and the Sun. The forces of each celestial body on the spacecraft are easily expressed using Newton's law of gravitation.

a_\mathrm{g} = GM_\odot \frac{ r_\odot - r }{{\Vert r_\odot - r \Vert}^3} +GM_\mathrm{Earth} \frac{ r_\mathrm{Earth} - r }{{\Vert r_\mathrm{Earth} - r \Vert}^3} +GM_\mathrm{Moon} \frac{ r_\mathrm{Moon} - r }{{\Vert r_\mathrm{Moon} - r \Vert}^3}

$\odot$ is a symbol used to denote the Sun.

Other perturbations also exist our spacecraft in space. With our model, we can just add more perturbation terms.

The Earth's somewhat oblate shape changes its gravitational field a small amount. This is model using the spherical harmonic model of Earth's gravity. By adding the second harmonic of Earth's gravity, sometimes referred to as the J2, we take into account the equatorial bulge caused by Earth's rotation around its axis.

J_2 \approx 1.08263 \times 10^{-3}

I greatly simplified the gravitational potential by considering the spacecraft to be equatorial (which is definitely not true).

a_{J_2} = \frac{3}{2} J_2 \, G M_\mathrm{Earth} R_\mathrm{Earth}^2\frac{r_\mathrm{Earth} - r}{\Vert r_\mathrm{Earth} - r \Vert ^5}

The Sun's photons impact a (very small) force on the spacecraft, due to the fact that photons have momentum that is transferred to our object when they land. This is called radiation pressure (SRP), or more specifically for the Sun, the solar radiation pressure. Its value near the Earth (at 1 AU) is roughly:

P_\odot \approx 4.56 \, \mu \mathrm{Pa}

This is the pressure value at 1 AU, near the Earth. To find the pressure value at any distance from the Sun, we use the inverse-square law. And to obtain the acceleration of our spacecraft, we multiply the pressure by its cross-section area $A$ and its reflectivity coefficient $C_\mathrm{R}$ , which I estimated to be:

A \approx \pi R^2 + f\, n\, w h \approx 48 \, \mathrm{m}^2 \newline \begin{aligned} \text{where } & R &&\approx 2.5 \, \mathrm{m} && \text{(capsule radius)} \\ & n &&= 4 && \text{(solar array wings)} \\ & w &&\approx 7 \, \mathrm{m}, \quad h \approx 2 \, \mathrm{m} && \text{(panels dimensions)} \\ & f &&\approx 0.5 && \text{(mean projection factor)} \end{aligned} \\ C_\mathrm{R} \approx 1.3

Hence, the Sun radiation pressure (SRP):

a_\mathrm{SRP} = \frac{P_\odot C_\mathrm{R} A}{m} \frac{(1 \, \mathrm{AU})^2}{{\Vert r_\odot - r \Vert}^2}

We now have all the perturbation terms we want to model. In our simulation, we can just add them all together:

a = a_\mathrm{g} + a_{J_2} + a_\mathrm{SRP}

Integrating numerically

In order to actually perform the simulation, we need to integrate the perturbation model.

There are many methods. The one I used is the Dormand-Prince method, which is part of the family of the Runge-Kutta methods.

Using the method is very simple. At each step, we compute the next step by using the form:

y_{n+1} = y_n + h \sum_{i=1}^{s} b_i k_i \\ \begin{aligned} \text{with the state vector } & y_n && = \begin{pmatrix} r_{n} \\ v_{n}\end{pmatrix} \\ \text{the time step } & h \\ \text{and the functions } & k_1 && = f(t_n, y_n), \\ & k_2 && = f(t_n+c_2h, y_n+(a_{21}k_1)h), \\ & k_3 && = f(t_n+c_3h, y_n+(a_{31}k_1+a_{32}k_2)h), \\ & \ \vdots \\ & k_s && = f(t_n+c_sh, y_n+(a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1})h). \\ & f(t,y) &&= \frac{dy}{dt} = \begin{pmatrix} v \\ a_g + a_{J_2} + a_{SRP} \end{pmatrix} \\ \end{aligned}

And the coefficients $c_i, a_{ij}, b_i$ are given by a table. In the Dormand-Prince method, we use the coefficients

\begin{array}{c|ccccccc} 0 \\[0.5em] \frac{1}{5} & \frac{1}{5} \\[0.5em] \frac{3}{10} & \frac{3}{40} & \frac{9}{40} \\[0.5em] \frac{4}{5} & \frac{44}{45} & -\frac{56}{15} & \frac{32}{9} \\[0.5em] \frac{8}{9} & \frac{19372}{6561} & -\frac{25360}{2187} & \frac{64448}{6561} & -\frac{212}{729} \\[0.5em] 1 & \frac{9017}{3168} & -\frac{355}{33} & \frac{46732}{5247} & \frac{49}{176} & -\frac{5103}{18656} \\[0.5em] 1 & \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} \\[0.3em] \hline \rule{0pt}{1.2em} & \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} & 0 \\[0.5em] & \frac{5179}{57600} & 0 & \frac{7571}{16695} & \frac{393}{640} & -\frac{92097}{339200} & \frac{187}{2100} & \frac{1}{40} \end{array}