A Deep Dive into the Hohmann Transfer with Python
Every kilogram of propellant a spacecraft carries is a kilogram it can’t use for payload, instruments, or crew. That’s why orbital mechanics engineers obsess over a single quantity: Δv (delta-v), the change in velocity required to move a spacecraft from one orbit to another. Less Δv means less fuel, which means cheaper, lighter, and more capable missions.
In this article, we’ll explore the Hohmann transfer orbit — the classic, fuel-optimal way to move between two circular orbits using just two engine burns — and then push further into the bi-elliptic transfer, a three-burn strategy that can actually beat Hohmann when the orbits are far enough apart. We’ll build everything in Python, verify the math numerically, visualize the trajectories in 3D, and benchmark a slow implementation against a fast, vectorized one.
1. The Physics of a Hohmann Transfer
A Hohmann transfer connects two coplanar circular orbits (radii $r_1$ and $r_2$) using an elliptical “bridge” orbit whose perigee touches the inner orbit and whose apogee touches the outer orbit.
The semi-major axis of this transfer ellipse is:
$$
a_t = \frac{r_1 + r_2}{2}
$$
The circular velocities at each orbit come from the vis-viva equation applied to a circular path:
$$
v_1 = \sqrt{\frac{\mu}{r_1}}, \qquad v_2 = \sqrt{\frac{\mu}{r_2}}
$$
where $\mu = GM$ is the gravitational parameter of the central body (for Earth, $\mu \approx 398600.4418\ \text{km}^3/\text{s}^2$).
The velocities on the transfer ellipse itself, at perigee ($r_1$) and apogee ($r_2$), come from the full vis-viva equation:
$$
v_p = \sqrt{\mu \left( \frac{2}{r_1} - \frac{1}{a_t} \right)}, \qquad
v_a = \sqrt{\mu \left( \frac{2}{r_2} - \frac{1}{a_t} \right)}
$$
The two burns are:
$$
\Delta v_1 = v_p - v_1 \quad \text{(speed up to leave the inner circular orbit)}
$$
$$
\Delta v_2 = v_2 - v_a \quad \text{(speed up to circularize at the outer orbit)}
$$
$$
\Delta v_{\text{total}} = |\Delta v_1| + |\Delta v_2|
$$
The time of flight is simply half the period of the transfer ellipse:
$$
t_{\text{transfer}} = \pi \sqrt{\frac{a_t^3}{\mu}}
$$
This two-burn strategy is provably the minimum-Δv way to move between two circular orbits — as long as the ratio $r_2/r_1$ stays below about 11.94. Beyond that ratio, something surprising happens.
2. Bi-Elliptic Transfer: Beating Hohmann with a Detour
If the target orbit is far enough away, it can actually be cheaper to overshoot to a very large intermediate radius $r_b$, coast out on one ellipse, then coast back in on a second ellipse to reach $r_2$ — using three burns instead of two.
$$
a_1 = \frac{r_1 + r_b}{2}, \qquad a_2 = \frac{r_b + r_2}{2}
$$
$$
\Delta v_1 = \sqrt{\mu\left(\frac{2}{r_1}-\frac{1}{a_1}\right)} - \sqrt{\frac{\mu}{r_1}}
$$
$$
\Delta v_2 = \sqrt{\mu\left(\frac{2}{r_b}-\frac{1}{a_2}\right)} - \sqrt{\mu\left(\frac{2}{r_b}-\frac{1}{a_1}\right)}
$$
$$
\Delta v_3 = \sqrt{\frac{\mu}{r_2}} - \sqrt{\mu\left(\frac{2}{r_2}-\frac{1}{a_2}\right)}
$$
The intermediate radius $r_b$ is a free design variable — and finding the $r_b$ that minimizes $\Delta v_1 + \Delta v_2 + \Delta v_3$ is a genuine, textbook optimization problem. That’s exactly what we’ll hand to scipy.optimize below.
3. Worked Example
We’ll compute a real mission profile: launching from a 300 km altitude LEO and transferring to geostationary orbit (GEO, r = 42,164 km). We’ll also run a second scenario — LEO to a 100,000 km orbit — where the ratio $r_2/r_1 \approx 15$ exceeds the 11.94 threshold, letting us demonstrate the bi-elliptic advantage.
4. Full Python Source Code
This script is self-contained and runs top-to-bottom without modification.
1 | # ============================================================= |
5. Code Walkthrough
Section 1 — Constants. MU_EARTH and R_EARTH are the two numbers every geocentric orbital calculation depends on. Keeping them as named constants (rather than magic numbers scattered through the code) makes the script easy to adapt to other bodies — swap in the Moon’s or Mars’ $\mu$ and everything else still works.
Section 2 — hohmann_transfer(). This function is a direct, literal translation of the equations from Part 1. It returns a dictionary rather than a tuple so that every value is self-documenting when you inspect the result later (result['dv_total'] is much clearer than result[5]).
Section 3 — The LEO→GEO example. This is the “hello world” of orbital transfers. Notice that $\Delta v_1 \approx 2.4\ \text{km/s}$ and $\Delta v_2 \approx 1.5\ \text{km/s}$ — the first burn (leaving LEO) is always the expensive one because the spacecraft is moving fastest there, and a given velocity change buys the least “orbit-shape change per km/s” close to the planet.
Section 4 — plot_hohmann_3d(). The trajectory is computed with the standard polar conic-section formula $r(\theta) = \dfrac{a(1-e^2)}{1+e\cos\theta}$, using the eccentricity of the transfer ellipse derived from $r_1$ and $r_2$. Even though the transfer is physically planar, plotting it on 3D axes (with z=0) alongside a rendered Earth sphere makes the geometry — and the 180° half-ellipse “bridge” between the two circles — far more intuitive than a flat 2D plot.
Section 5 — Bi-elliptic transfer and scipy.optimize.minimize_scalar. This is the heart of the “minimize fuel consumption” theme. bi_elliptic_transfer() takes the intermediate radius $r_b$ as an explicit parameter, and minimize_scalar searches over it (bounded between just above $r_2$ and 25× $r_2$) to find the $r_b$ that truly minimizes total Δv. Because we chose $r_2/r_1 \approx 15$, which is above the well-known 11.94 crossover ratio, the optimizer should find a bi-elliptic solution that beats plain Hohmann — a nice numerical confirmation of the analytical theory.
Section 6 — Naive loop vs. vectorized NumPy. Scanning a $60 \times 60$ grid of orbit pairs with a Python-level double for loop (3,600 function calls, each doing several sqrt calls) is the classic “slow” pattern. The vectorized version replaces the entire loop with array broadcasting: np.meshgrid builds the 2D grid of $r_1, r_2$ pairs at once, and every formula (V1, V2, VP, VA, DV_TOTAL) is evaluated on the whole array simultaneously using NumPy’s compiled C backend. The final sanity check (max numerical diff) confirms both methods agree to floating-point precision — the speedup is “free,” with zero loss of accuracy.
Section 7 — 3D Delta-v surface. This turns the entire vectorized grid into a single 3D surface, letting you visually scan “which combinations of $r_1, r_2$ are cheap or expensive to connect.” The steep rise near small $r_1$ visually confirms why leaving a low, tight orbit is always the dominant fuel cost.
6. Where to See the Results
Run all seven sections in order in a single Colab cell (or split across cells — the functions and variables carry over naturally). Here’s where each output belongs when you paste your results back into this post:
Console output — Hohmann LEO→GEO and Bi-elliptic minimization:
===== Hohmann Transfer: LEO -> GEO ===== r1 (LEO radius) : 6678.14 km r2 (GEO radius) : 42164.00 km Transfer semi-major a : 24421.07 km v1 (circular @ LEO) : 7.7258 km/s v_p (transfer perigee) : 10.1515 km/s Delta-v1 (1st burn) : 2.4257 km/s v_a (transfer apogee) : 1.6078 km/s v2 (circular @ GEO) : 3.0747 km/s Delta-v2 (2nd burn) : 1.4668 km/s Total Delta-v : 3.8926 km/s Transfer time : 5.2750 hours
Image — 3D Hohmann transfer trajectory (Section 4):
Console output — Speed comparison (Section 6):
===== Delta-v Minimization: Hohmann vs Bi-Elliptic ===== r1 = 6678.14 km, r2 = 100000.00 km, ratio r2/r1 = 14.97 Hohmann total Delta-v : 4.1427 km/s Optimal intermediate radius : 2499999.94 km Bi-elliptic total Delta-v : 4.0393 km/s Delta-v saving : 0.1034 km/s ===== Speed Comparison: Naive Loop vs Vectorized NumPy ===== Grid size : 60 x 60 = 3600 points Naive loop time : 87.67 ms Vectorized NumPy time : 0.92 ms Speedup : 95.7x Max numerical diff : 0.00e+00 km/s
Image — 3D Delta-v surface (Section 7):
7. Takeaways
The Hohmann transfer remains the workhorse of orbital mechanics because it’s provably optimal for the vast majority of real missions — but as this example shows, “optimal” always comes with an asterisk: the 11.94 radius-ratio threshold where a three-burn bi-elliptic detour quietly becomes the cheaper option. Turning that insight into code is straightforward once the vis-viva equation is in hand, and letting scipy.optimize search the design space numerically is a good habit to build even when — as here — you could also derive the answer by hand. The vectorized NumPy pattern in Section 6 is worth keeping in your toolbox too: any time you’re evaluating the same closed-form physics formula across a grid of parameters, reach for broadcasting before you reach for a loop.
