Skip to content

LEO to GEO Hohmann Transfer

In this example we'll compute and visualize a Hohmann transfer from a low Earth orbit (LEO) at 400 km altitude to geostationary orbit (GEO) at 35,786 km altitude. The Hohmann transfer is the most fuel-efficient two-impulse maneuver for transferring between coplanar circular orbits, making it fundamental to spacecraft mission design.

We'll use the NumericalOrbitPropagator with event callbacks to apply impulsive velocity changes at perigee (departure) and apogee (arrival). The example demonstrates the complete workflow from calculating the required delta-v values to visualizing the orbit geometry and trajectory profiles.

Problem Setup

Our mission is to transfer a spacecraft from a circular LEO parking orbit at 400 km altitude to geostationary orbit at 35,786 km altitude. At GEO altitude, a satellite with zero inclination orbits the Earth once per sidereal day, appearing stationary relative to the ground - ideal for communications and weather satellites.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import os
import pathlib
import sys

import numpy as np
import plotly.graph_objects as go

import brahe as bh

bh.initialize_eop()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# LEO starting orbit: 400 km circular altitude
r_leo = bh.R_EARTH + 400e3  # meters

# GEO target orbit: 35,786 km altitude (geostationary)
r_geo = bh.R_EARTH + 35786e3  # meters

print(
    f"LEO radius: {r_leo / 1e3:.1f} km (altitude: {(r_leo - bh.R_EARTH) / 1e3:.0f} km)"
)
print(
    f"GEO radius: {r_geo / 1e3:.1f} km (altitude: {(r_geo - bh.R_EARTH) / 1e3:.0f} km)"
)

Hohmann Transfer Theory

Mathematical Background

The Hohmann transfer uses an elliptical transfer orbit that is tangent to both the initial and final circular orbits. The spacecraft performs two impulsive burns:

  1. First burn at perigee (LEO altitude): Increases velocity to enter the transfer ellipse
  2. Second burn at apogee (GEO altitude): Circularizes the orbit at the target altitude

The transfer orbit has its perigee at the initial orbit radius \(r_1\) and apogee at the final orbit radius \(r_2\).

Transfer orbit semi-major axis:

\[a_{transfer} = \frac{r_1 + r_2}{2}\]

Circular orbit velocity (from the vis-viva equation for \(e = 0\)):

\[v_{circ} = \sqrt{\frac{\mu}{r}}\]

Velocity at perigee and apogee of transfer ellipse (vis-viva equation):

\[v = \sqrt{\mu \left(\frac{2}{r} - \frac{1}{a}\right)}\]

Delta-v for each burn:

\[\Delta v_1 = v_{perigee,transfer} - v_{LEO} = \sqrt{\mu \left(\frac{2}{r_1} - \frac{1}{a_{transfer}}\right)} - \sqrt{\frac{\mu}{r_1}}\]
\[\Delta v_2 = v_{GEO} - v_{apogee,transfer} = \sqrt{\frac{\mu}{r_2}} - \sqrt{\mu \left(\frac{2}{r_2} - \frac{1}{a_{transfer}}\right)}\]

Transfer time (half the orbital period of the transfer ellipse):

\[T_{transfer} = \pi \sqrt{\frac{a_{transfer}^3}{\mu}}\]

Delta-V Calculations

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# Circular orbit velocities (vis-viva equation for circular orbit: v = sqrt(mu/r))
v_leo = np.sqrt(bh.GM_EARTH / r_leo)
v_geo = np.sqrt(bh.GM_EARTH / r_geo)

# Transfer orbit parameters
# Semi-major axis is the average of the two radii
a_transfer = (r_leo + r_geo) / 2

# Eccentricity of transfer ellipse
e_transfer = (r_geo - r_leo) / (r_geo + r_leo)

# Velocities on transfer orbit using vis-viva equation: v^2 = mu(2/r - 1/a)
v_perigee_transfer = np.sqrt(bh.GM_EARTH * (2 / r_leo - 1 / a_transfer))
v_apogee_transfer = np.sqrt(bh.GM_EARTH * (2 / r_geo - 1 / a_transfer))

# Delta-v magnitudes
dv1 = v_perigee_transfer - v_leo  # First burn: prograde at perigee (LEO)
dv2 = v_geo - v_apogee_transfer  # Second burn: prograde at apogee (GEO)

# Transfer time: half the period of the transfer ellipse
transfer_time = np.pi * np.sqrt(a_transfer**3 / bh.GM_EARTH)

print("\nHohmann Transfer Parameters:")
print(f"  Transfer semi-major axis: {a_transfer / 1e3:.1f} km")
print(f"  Transfer eccentricity:    {e_transfer:.4f}")
print(f"  LEO circular velocity:    {v_leo / 1e3:.3f} km/s")
print(f"  GEO circular velocity:    {v_geo / 1e3:.3f} km/s")
print(f"  First burn (perigee):     {dv1 / 1e3:.3f} km/s")
print(f"  Second burn (apogee):     {dv2 / 1e3:.3f} km/s")
print(f"  Total delta-v:            {(dv1 + dv2) / 1e3:.3f} km/s")
print(f"  Transfer time:            {transfer_time / 3600:.2f} hours")

The total delta-v of approximately 3.85 km/s is substantial - this is why geostationary satellites require large launch vehicles or are launched into geostationary transfer orbits (GTO) where the rocket provides the first burn and the satellite provides the circularization burn.

Implementation

Initial State

We create an initial circular orbit at LEO altitude. The spacecraft starts at the point that will become the perigee of the transfer orbit:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# Create initial epoch
epoch = bh.Epoch.from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)

# Initial state: circular LEO orbit (at perigee of transfer)
# Keplerian elements: [a, e, i, raan, argp, M]
oe_initial = np.array([r_leo, 0.0001, 0.0, 0.0, 0.0, 0.0])
state_initial = bh.state_koe_to_eci(oe_initial, bh.AngleFormat.DEGREES)

print("\nInitial State (ECI):")
print(
    f"  Position: [{state_initial[0] / 1e3:.1f}, {state_initial[1] / 1e3:.1f}, {state_initial[2] / 1e3:.1f}] km"
)
print(
    f"  Velocity: [{state_initial[3] / 1e3:.3f}, {state_initial[4] / 1e3:.3f}, {state_initial[5] / 1e3:.3f}] km/s"
)

Event Callbacks for Impulsive Maneuvers

One way to implement the impulsive burns is by propagating to the burn times, extracting the state, applying the delta-v, and restarting propagation with a new state and propagator. However, brahe provides a cleaner approach using event callbacks. An event callback is a function that is invoked when a specified event occurs during propagation. In this case, we use TimeEvent detectors to trigger the burns at the calculated transfer times. Each callback receives the event epoch and current state, and returns the modified state along with an EventAction indicating whether to continue propagation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Define event callbacks for impulsive maneuvers
# Each callback receives the event epoch and current state,
# and returns (new_state, EventAction)


def first_burn_callback(event_epoch, event_state):
    """Apply first delta-v at departure (prograde burn at perigee)."""
    new_state = event_state.copy()
    # Add delta-v in velocity direction (prograde)
    v = event_state[3:6]
    v_hat = v / np.linalg.norm(v)
    new_state[3:6] += dv1 * v_hat
    print(f"First burn applied: dv = {dv1 / 1e3:.3f} km/s (prograde)")
    return (new_state, bh.EventAction.CONTINUE)


def second_burn_callback(event_epoch, event_state):
    """Apply second delta-v at arrival (prograde burn at apogee)."""
    new_state = event_state.copy()
    v = event_state[3:6]
    v_hat = v / np.linalg.norm(v)
    new_state[3:6] += dv2 * v_hat
    print(f"Second burn applied: dv = {dv2 / 1e3:.3f} km/s (prograde)")
    return (new_state, bh.EventAction.CONTINUE)

The callbacks apply delta-v in the prograde direction (along the velocity vector) to raise the orbit. Returning EventAction.CONTINUE tells the propagator to proceed with the modified state.

Single Propagator with Events

We use a single propagator with TimeEvent detectors to trigger the burns at the appropriate times. This is cleaner than multi-stage propagation for simple maneuver sequences:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
# Create a single propagator with event-based maneuvers
# This is cleaner than multi-stage propagation for simple burns

# Timing
burn1_time_s = 1.0  # First burn shortly after start
burn2_time_s = burn1_time_s + transfer_time  # Second burn at apogee
geo_period = bh.orbital_period(r_geo)
total_time = burn2_time_s + geo_period  # Continue for one GEO orbit

print("\nPropagation Timeline:")
print(f"  Burn 1 at t = {burn1_time_s:.1f} s")
print(f"  Burn 2 at t = {burn2_time_s / 3600:.2f} hours")
print(f"  Total simulation: {total_time / 3600:.2f} hours")

# Create the propagator with two-body dynamics
prop = bh.NumericalOrbitPropagator(
    epoch,
    state_initial,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.two_body(),
    None,
)

# Add time-based events for the burns
event1 = bh.TimeEvent(epoch + burn1_time_s, "First Burn").with_callback(
    first_burn_callback
)
event2 = bh.TimeEvent(epoch + burn2_time_s, "Second Burn").with_callback(
    second_burn_callback
)

prop.add_event_detector(event1)
prop.add_event_detector(event2)

# Propagate through both burns plus one GEO orbit
print("\nPropagating...")
prop.propagate_to(epoch + total_time)
print("  Complete!")

# Verify final orbit
final_koe = prop.state_koe_osc(prop.current_epoch, bh.AngleFormat.DEGREES)
final_altitude = final_koe[0] - bh.R_EARTH
print("\nFinal GEO Orbit:")
print(f"  Semi-major axis: {final_koe[0] / 1e3:.1f} km")
print(f"  Altitude:        {final_altitude / 1e3:.1f} km (target: 35786 km)")
print(f"  Eccentricity:    {final_koe[1]:.6f}")

Why Two-Body Dynamics?

We use the two-body force model (ForceModelConfig.two_body()) for this example because it provides the idealized Keplerian motion that matches the analytical Hohmann transfer theory. In practice, perturbations from Earth's oblateness, atmospheric drag (at LEO), and third-body effects would cause deviations that require trajectory correction maneuvers.

Sampling the Trajectory

We sample the trajectory at regular intervals to create the visualization data. The single propagator stores the complete trajectory including the state discontinuities from the burns:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Sample trajectory data for plotting
# The single propagator stores the complete trajectory with maneuvers
times_hours = []
altitudes_km = []
velocities_km_s = []

dt = 60.0  # 1-minute sampling
t = 0.0

while t <= total_time:
    current_epoch = epoch + t
    state = prop.state_eci(current_epoch)

    r_mag = np.linalg.norm(state[:3])
    v_mag = np.linalg.norm(state[3:6])

    times_hours.append(t / 3600.0)
    altitudes_km.append((r_mag - bh.R_EARTH) / 1e3)
    velocities_km_s.append(v_mag / 1e3)

    t += dt

print(f"\nSampled {len(times_hours)} trajectory points")

Visualizations

Orbit Geometry

The following plot shows a top-down view of the transfer geometry. The spacecraft departs from the LEO parking orbit (dashed blue), follows the transfer ellipse (solid red) for half an orbit, and arrives at GEO (dashed green):

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
# Create 2D top-down view of orbit geometry
theta = np.linspace(0, 2 * np.pi, 200)

# Earth circle
earth_x = bh.R_EARTH * np.cos(theta) / 1e3
earth_y = bh.R_EARTH * np.sin(theta) / 1e3

# LEO orbit circle
leo_x = r_leo * np.cos(theta) / 1e3
leo_y = r_leo * np.sin(theta) / 1e3

# GEO orbit circle
geo_x = r_geo * np.cos(theta) / 1e3
geo_y = r_geo * np.sin(theta) / 1e3

# Transfer ellipse (only upper half: theta from 0 to pi)
theta_transfer = np.linspace(0, np.pi, 100)
p_transfer = a_transfer * (1 - e_transfer**2)
r_transfer = p_transfer / (1 + e_transfer * np.cos(theta_transfer))
transfer_x = r_transfer * np.cos(theta_transfer) / 1e3
transfer_y = r_transfer * np.sin(theta_transfer) / 1e3

fig_geometry = go.Figure()

# Earth (filled)
fig_geometry.add_trace(
    go.Scatter(
        x=earth_x.tolist(),
        y=earth_y.tolist(),
        fill="toself",
        fillcolor="lightblue",
        line=dict(color="steelblue", width=1),
        name="Earth",
        hoverinfo="name",
    )
)

# LEO orbit
fig_geometry.add_trace(
    go.Scatter(
        x=leo_x.tolist(),
        y=leo_y.tolist(),
        mode="lines",
        line=dict(color="blue", width=2, dash="dash"),
        name=f"LEO ({(r_leo - bh.R_EARTH) / 1e3:.0f} km)",
    )
)

# GEO orbit
fig_geometry.add_trace(
    go.Scatter(
        x=geo_x.tolist(),
        y=geo_y.tolist(),
        mode="lines",
        line=dict(color="green", width=2, dash="dash"),
        name=f"GEO ({(r_geo - bh.R_EARTH) / 1e3:.0f} km)",
    )
)

# Transfer ellipse
fig_geometry.add_trace(
    go.Scatter(
        x=transfer_x.tolist(),
        y=transfer_y.tolist(),
        mode="lines",
        line=dict(color="red", width=3),
        name="Transfer Orbit",
    )
)

# Burn 1 marker (at LEO, right side)
fig_geometry.add_trace(
    go.Scatter(
        x=[r_leo / 1e3],
        y=[0],
        mode="markers+text",
        marker=dict(size=15, color="red", symbol="star"),
        text=[f"Burn 1<br>{dv1 / 1e3:.2f} km/s"],
        textposition="bottom right",
        textfont=dict(size=10),
        name="Burn 1",
        showlegend=False,
    )
)

# Burn 2 marker (at GEO, left side - apogee)
fig_geometry.add_trace(
    go.Scatter(
        x=[-r_geo / 1e3],
        y=[0],
        mode="markers+text",
        marker=dict(size=15, color="red", symbol="star"),
        text=[f"Burn 2<br>{dv2 / 1e3:.2f} km/s"],
        textposition="top left",
        textfont=dict(size=10),
        name="Burn 2",
        showlegend=False,
    )
)

fig_geometry.update_layout(
    title="LEO to GEO Hohmann Transfer Geometry",
    xaxis_title="X (km)",
    yaxis_title="Y (km)",
    yaxis_scaleanchor="x",
    showlegend=True,
    legend=dict(x=0.02, y=0.98),
    height=600,
    margin=dict(l=60, r=40, t=60, b=60),
)

Altitude Profile

The altitude profile shows the spacecraft climbing from 400 km to nearly 36,000 km over the ~5.3 hour transfer:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# Create altitude vs time plot
fig_altitude = go.Figure()

fig_altitude.add_trace(
    go.Scatter(
        x=times_hours,
        y=altitudes_km,
        mode="lines",
        line=dict(color="blue", width=2),
        name="Altitude",
    )
)

# Reference lines for initial and target altitudes
fig_altitude.add_hline(
    y=400,
    line_dash="dash",
    line_color="gray",
    annotation_text="LEO: 400 km",
    annotation_position="top right",
)

fig_altitude.add_hline(
    y=35786,
    line_dash="dash",
    line_color="gray",
    annotation_text="GEO: 35,786 km",
    annotation_position="top right",
)

# Burn markers
burn1_time_hr = burn1_time_s / 3600.0
burn2_time_hr = burn2_time_s / 3600.0

fig_altitude.add_vline(
    x=burn1_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 1: {dv1 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_altitude.add_vline(
    x=burn2_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 2: {dv2 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_altitude.update_layout(
    title="Altitude During Hohmann Transfer",
    xaxis_title="Time (hours)",
    yaxis_title="Altitude (km)",
    showlegend=False,
    height=500,
    margin=dict(l=60, r=40, t=60, b=60),
)

Velocity Profile

The velocity profile reveals the characteristic behavior of elliptical orbits - the spacecraft is fastest at perigee and slowest at apogee. The impulsive burns appear as discontinuities:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# Create velocity vs time plot
fig_velocity = go.Figure()

fig_velocity.add_trace(
    go.Scatter(
        x=times_hours,
        y=velocities_km_s,
        mode="lines",
        line=dict(color="blue", width=2),
        name="Velocity",
    )
)

# Reference lines for circular orbit velocities
fig_velocity.add_hline(
    y=v_leo / 1e3,
    line_dash="dash",
    line_color="gray",
    annotation_text=f"LEO: {v_leo / 1e3:.2f} km/s",
    annotation_position="top right",
)

fig_velocity.add_hline(
    y=v_geo / 1e3,
    line_dash="dash",
    line_color="gray",
    annotation_text=f"GEO: {v_geo / 1e3:.2f} km/s",
    annotation_position="bottom right",
)

# Burn markers
fig_velocity.add_vline(
    x=burn1_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 1: +{dv1 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_velocity.add_vline(
    x=burn2_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 2: +{dv2 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_velocity.update_layout(
    title="Velocity During Hohmann Transfer",
    xaxis_title="Time (hours)",
    yaxis_title="Velocity (km/s)",
    showlegend=False,
    height=500,
    margin=dict(l=60, r=40, t=60, b=60),
)

Transfer Summary

Parameter Value
Initial altitude 400 km (LEO)
Final altitude 35,786 km (GEO)
First burn \(\Delta v_1\) 2.40 km/s
Second burn \(\Delta v_2\) 1.46 km/s
Total \(\Delta v\) 3.85 km/s
Transfer time 5.29 hours

Full Code Example

geo_hohmann_transfer.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
import os
import pathlib
import sys

import numpy as np
import plotly.graph_objects as go

import brahe as bh

bh.initialize_eop()

# Configuration for output files
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# LEO starting orbit: 400 km circular altitude
r_leo = bh.R_EARTH + 400e3  # meters

# GEO target orbit: 35,786 km altitude (geostationary)
r_geo = bh.R_EARTH + 35786e3  # meters

print(
    f"LEO radius: {r_leo / 1e3:.1f} km (altitude: {(r_leo - bh.R_EARTH) / 1e3:.0f} km)"
)
print(
    f"GEO radius: {r_geo / 1e3:.1f} km (altitude: {(r_geo - bh.R_EARTH) / 1e3:.0f} km)"
)

# Circular orbit velocities (vis-viva equation for circular orbit: v = sqrt(mu/r))
v_leo = np.sqrt(bh.GM_EARTH / r_leo)
v_geo = np.sqrt(bh.GM_EARTH / r_geo)

# Transfer orbit parameters
# Semi-major axis is the average of the two radii
a_transfer = (r_leo + r_geo) / 2

# Eccentricity of transfer ellipse
e_transfer = (r_geo - r_leo) / (r_geo + r_leo)

# Velocities on transfer orbit using vis-viva equation: v^2 = mu(2/r - 1/a)
v_perigee_transfer = np.sqrt(bh.GM_EARTH * (2 / r_leo - 1 / a_transfer))
v_apogee_transfer = np.sqrt(bh.GM_EARTH * (2 / r_geo - 1 / a_transfer))

# Delta-v magnitudes
dv1 = v_perigee_transfer - v_leo  # First burn: prograde at perigee (LEO)
dv2 = v_geo - v_apogee_transfer  # Second burn: prograde at apogee (GEO)

# Transfer time: half the period of the transfer ellipse
transfer_time = np.pi * np.sqrt(a_transfer**3 / bh.GM_EARTH)

print("\nHohmann Transfer Parameters:")
print(f"  Transfer semi-major axis: {a_transfer / 1e3:.1f} km")
print(f"  Transfer eccentricity:    {e_transfer:.4f}")
print(f"  LEO circular velocity:    {v_leo / 1e3:.3f} km/s")
print(f"  GEO circular velocity:    {v_geo / 1e3:.3f} km/s")
print(f"  First burn (perigee):     {dv1 / 1e3:.3f} km/s")
print(f"  Second burn (apogee):     {dv2 / 1e3:.3f} km/s")
print(f"  Total delta-v:            {(dv1 + dv2) / 1e3:.3f} km/s")
print(f"  Transfer time:            {transfer_time / 3600:.2f} hours")

# Create initial epoch
epoch = bh.Epoch.from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)

# Initial state: circular LEO orbit (at perigee of transfer)
# Keplerian elements: [a, e, i, raan, argp, M]
oe_initial = np.array([r_leo, 0.0001, 0.0, 0.0, 0.0, 0.0])
state_initial = bh.state_koe_to_eci(oe_initial, bh.AngleFormat.DEGREES)

print("\nInitial State (ECI):")
print(
    f"  Position: [{state_initial[0] / 1e3:.1f}, {state_initial[1] / 1e3:.1f}, {state_initial[2] / 1e3:.1f}] km"
)
print(
    f"  Velocity: [{state_initial[3] / 1e3:.3f}, {state_initial[4] / 1e3:.3f}, {state_initial[5] / 1e3:.3f}] km/s"
)

# Define event callbacks for impulsive maneuvers
# Each callback receives the event epoch and current state,
# and returns (new_state, EventAction)


def first_burn_callback(event_epoch, event_state):
    """Apply first delta-v at departure (prograde burn at perigee)."""
    new_state = event_state.copy()
    # Add delta-v in velocity direction (prograde)
    v = event_state[3:6]
    v_hat = v / np.linalg.norm(v)
    new_state[3:6] += dv1 * v_hat
    print(f"First burn applied: dv = {dv1 / 1e3:.3f} km/s (prograde)")
    return (new_state, bh.EventAction.CONTINUE)


def second_burn_callback(event_epoch, event_state):
    """Apply second delta-v at arrival (prograde burn at apogee)."""
    new_state = event_state.copy()
    v = event_state[3:6]
    v_hat = v / np.linalg.norm(v)
    new_state[3:6] += dv2 * v_hat
    print(f"Second burn applied: dv = {dv2 / 1e3:.3f} km/s (prograde)")
    return (new_state, bh.EventAction.CONTINUE)



# Create a single propagator with event-based maneuvers
# This is cleaner than multi-stage propagation for simple burns

# Timing
burn1_time_s = 1.0  # First burn shortly after start
burn2_time_s = burn1_time_s + transfer_time  # Second burn at apogee
geo_period = bh.orbital_period(r_geo)
total_time = burn2_time_s + geo_period  # Continue for one GEO orbit

print("\nPropagation Timeline:")
print(f"  Burn 1 at t = {burn1_time_s:.1f} s")
print(f"  Burn 2 at t = {burn2_time_s / 3600:.2f} hours")
print(f"  Total simulation: {total_time / 3600:.2f} hours")

# Create the propagator with two-body dynamics
prop = bh.NumericalOrbitPropagator(
    epoch,
    state_initial,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.two_body(),
    None,
)

# Add time-based events for the burns
event1 = bh.TimeEvent(epoch + burn1_time_s, "First Burn").with_callback(
    first_burn_callback
)
event2 = bh.TimeEvent(epoch + burn2_time_s, "Second Burn").with_callback(
    second_burn_callback
)

prop.add_event_detector(event1)
prop.add_event_detector(event2)

# Propagate through both burns plus one GEO orbit
print("\nPropagating...")
prop.propagate_to(epoch + total_time)
print("  Complete!")

# Verify final orbit
final_koe = prop.state_koe_osc(prop.current_epoch, bh.AngleFormat.DEGREES)
final_altitude = final_koe[0] - bh.R_EARTH
print("\nFinal GEO Orbit:")
print(f"  Semi-major axis: {final_koe[0] / 1e3:.1f} km")
print(f"  Altitude:        {final_altitude / 1e3:.1f} km (target: 35786 km)")
print(f"  Eccentricity:    {final_koe[1]:.6f}")

# Sample trajectory data for plotting
# The single propagator stores the complete trajectory with maneuvers
times_hours = []
altitudes_km = []
velocities_km_s = []

dt = 60.0  # 1-minute sampling
t = 0.0

while t <= total_time:
    current_epoch = epoch + t
    state = prop.state_eci(current_epoch)

    r_mag = np.linalg.norm(state[:3])
    v_mag = np.linalg.norm(state[3:6])

    times_hours.append(t / 3600.0)
    altitudes_km.append((r_mag - bh.R_EARTH) / 1e3)
    velocities_km_s.append(v_mag / 1e3)

    t += dt

print(f"\nSampled {len(times_hours)} trajectory points")

# Create 2D top-down view of orbit geometry
theta = np.linspace(0, 2 * np.pi, 200)

# Earth circle
earth_x = bh.R_EARTH * np.cos(theta) / 1e3
earth_y = bh.R_EARTH * np.sin(theta) / 1e3

# LEO orbit circle
leo_x = r_leo * np.cos(theta) / 1e3
leo_y = r_leo * np.sin(theta) / 1e3

# GEO orbit circle
geo_x = r_geo * np.cos(theta) / 1e3
geo_y = r_geo * np.sin(theta) / 1e3

# Transfer ellipse (only upper half: theta from 0 to pi)
theta_transfer = np.linspace(0, np.pi, 100)
p_transfer = a_transfer * (1 - e_transfer**2)
r_transfer = p_transfer / (1 + e_transfer * np.cos(theta_transfer))
transfer_x = r_transfer * np.cos(theta_transfer) / 1e3
transfer_y = r_transfer * np.sin(theta_transfer) / 1e3

fig_geometry = go.Figure()

# Earth (filled)
fig_geometry.add_trace(
    go.Scatter(
        x=earth_x.tolist(),
        y=earth_y.tolist(),
        fill="toself",
        fillcolor="lightblue",
        line=dict(color="steelblue", width=1),
        name="Earth",
        hoverinfo="name",
    )
)

# LEO orbit
fig_geometry.add_trace(
    go.Scatter(
        x=leo_x.tolist(),
        y=leo_y.tolist(),
        mode="lines",
        line=dict(color="blue", width=2, dash="dash"),
        name=f"LEO ({(r_leo - bh.R_EARTH) / 1e3:.0f} km)",
    )
)

# GEO orbit
fig_geometry.add_trace(
    go.Scatter(
        x=geo_x.tolist(),
        y=geo_y.tolist(),
        mode="lines",
        line=dict(color="green", width=2, dash="dash"),
        name=f"GEO ({(r_geo - bh.R_EARTH) / 1e3:.0f} km)",
    )
)

# Transfer ellipse
fig_geometry.add_trace(
    go.Scatter(
        x=transfer_x.tolist(),
        y=transfer_y.tolist(),
        mode="lines",
        line=dict(color="red", width=3),
        name="Transfer Orbit",
    )
)

# Burn 1 marker (at LEO, right side)
fig_geometry.add_trace(
    go.Scatter(
        x=[r_leo / 1e3],
        y=[0],
        mode="markers+text",
        marker=dict(size=15, color="red", symbol="star"),
        text=[f"Burn 1<br>{dv1 / 1e3:.2f} km/s"],
        textposition="bottom right",
        textfont=dict(size=10),
        name="Burn 1",
        showlegend=False,
    )
)

# Burn 2 marker (at GEO, left side - apogee)
fig_geometry.add_trace(
    go.Scatter(
        x=[-r_geo / 1e3],
        y=[0],
        mode="markers+text",
        marker=dict(size=15, color="red", symbol="star"),
        text=[f"Burn 2<br>{dv2 / 1e3:.2f} km/s"],
        textposition="top left",
        textfont=dict(size=10),
        name="Burn 2",
        showlegend=False,
    )
)

fig_geometry.update_layout(
    title="LEO to GEO Hohmann Transfer Geometry",
    xaxis_title="X (km)",
    yaxis_title="Y (km)",
    yaxis_scaleanchor="x",
    showlegend=True,
    legend=dict(x=0.02, y=0.98),
    height=600,
    margin=dict(l=60, r=40, t=60, b=60),
)

# Create altitude vs time plot
fig_altitude = go.Figure()

fig_altitude.add_trace(
    go.Scatter(
        x=times_hours,
        y=altitudes_km,
        mode="lines",
        line=dict(color="blue", width=2),
        name="Altitude",
    )
)

# Reference lines for initial and target altitudes
fig_altitude.add_hline(
    y=400,
    line_dash="dash",
    line_color="gray",
    annotation_text="LEO: 400 km",
    annotation_position="top right",
)

fig_altitude.add_hline(
    y=35786,
    line_dash="dash",
    line_color="gray",
    annotation_text="GEO: 35,786 km",
    annotation_position="top right",
)

# Burn markers
burn1_time_hr = burn1_time_s / 3600.0
burn2_time_hr = burn2_time_s / 3600.0

fig_altitude.add_vline(
    x=burn1_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 1: {dv1 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_altitude.add_vline(
    x=burn2_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 2: {dv2 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_altitude.update_layout(
    title="Altitude During Hohmann Transfer",
    xaxis_title="Time (hours)",
    yaxis_title="Altitude (km)",
    showlegend=False,
    height=500,
    margin=dict(l=60, r=40, t=60, b=60),
)

# Create velocity vs time plot
fig_velocity = go.Figure()

fig_velocity.add_trace(
    go.Scatter(
        x=times_hours,
        y=velocities_km_s,
        mode="lines",
        line=dict(color="blue", width=2),
        name="Velocity",
    )
)

# Reference lines for circular orbit velocities
fig_velocity.add_hline(
    y=v_leo / 1e3,
    line_dash="dash",
    line_color="gray",
    annotation_text=f"LEO: {v_leo / 1e3:.2f} km/s",
    annotation_position="top right",
)

fig_velocity.add_hline(
    y=v_geo / 1e3,
    line_dash="dash",
    line_color="gray",
    annotation_text=f"GEO: {v_geo / 1e3:.2f} km/s",
    annotation_position="bottom right",
)

# Burn markers
fig_velocity.add_vline(
    x=burn1_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 1: +{dv1 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_velocity.add_vline(
    x=burn2_time_hr,
    line_dash="dot",
    line_color="red",
    annotation_text=f"Burn 2: +{dv2 / 1e3:.2f} km/s",
    annotation_position="top left",
)

fig_velocity.update_layout(
    title="Velocity During Hohmann Transfer",
    xaxis_title="Time (hours)",
    yaxis_title="Velocity (km/s)",
    showlegend=False,
    height=500,
    margin=dict(l=60, r=40, t=60, b=60),
)

# Validation
altitude_gain = final_altitude - (r_leo - bh.R_EARTH)
assert altitude_gain > 30000e3, f"Altitude gain too small: {altitude_gain / 1e3:.0f} km"
assert final_koe[1] < 0.01, f"Final orbit not circular enough: e = {final_koe[1]:.6f}"

print("\nExample validated successfully!")
print(f"  Altitude gain: {altitude_gain / 1e3:.0f} km")
print(f"  Final eccentricity: {final_koe[1]:.6f}")

See Also