Skip to content

Commit df91ee6

Browse files
Add Hohmann transfer simulation and animation
Implement Hohmann transfer orbit simulation with visualization.
1 parent 49ec67a commit df91ee6

1 file changed

Lines changed: 128 additions & 0 deletions

File tree

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
# Der Code ist in Englisch, da ich eher Englisch spreche.
2+
3+
import numpy as np
4+
import matplotlib.pyplot as plt
5+
from matplotlib.animation import FuncAnimation
6+
7+
# ================== CONSTANTS ==================
8+
mu = 3.986e14
9+
r_earth = 6.371e6
10+
g0 = 9.80665
11+
12+
Isp = 320
13+
m0 = 1000
14+
15+
r1 = 7e6
16+
r2 = 4.2e7
17+
18+
# ================== HOHMANN TRANSFER ==================
19+
a_t = (r1 + r2) / 2
20+
e_t = (r2 - r1) / (r1 + r2)
21+
22+
T1 = 2*np.pi*np.sqrt(r1**3 / mu)
23+
Tt = np.pi*np.sqrt(a_t**3 / mu)
24+
T2 = 2*np.pi*np.sqrt(r2**3 / mu)
25+
26+
N1, Nt, N2 = 200, 300, 200
27+
28+
# ================== ORBIT 1 ==================
29+
t1 = np.linspace(0, T1, N1)
30+
theta1 = 2*np.pi*t1/T1
31+
x1 = r1*np.cos(theta1)
32+
y1 = r1*np.sin(theta1)
33+
34+
# ================== TRANSFER ORBIT ==================
35+
def kepler(M, e):
36+
E = M.copy()
37+
for _ in range(10):
38+
E -= (E - e*np.sin(E) - M) / (1 - e*np.cos(E))
39+
return E
40+
41+
t_t = np.linspace(0, Tt, Nt)
42+
M = np.sqrt(mu/a_t**3) * t_t
43+
E = kepler(M, e_t)
44+
45+
theta_t = 2*np.arctan2(
46+
np.sqrt(1+e_t)*np.sin(E/2),
47+
np.sqrt(1-e_t)*np.cos(E/2)
48+
)
49+
50+
r_t = a_t * (1 - e_t*np.cos(E))
51+
x_t = r_t*np.cos(theta_t)
52+
y_t = r_t*np.sin(theta_t)
53+
54+
# ================== ORBIT 2 ==================
55+
T2 = 2*T2
56+
t2 = np.linspace(0, T2, N2)
57+
theta2 = theta_t[-1] + 2*np.pi*t2/T2
58+
x2 = r2*np.cos(theta2)
59+
y2 = r2*np.sin(theta2)
60+
61+
# ================== COMBINE ==================
62+
x = np.concatenate([x1, x_t, x2])
63+
y = np.concatenate([y1, y_t, y2])
64+
r = np.sqrt(x**2 + y**2)
65+
66+
# ================== SPEED & ENERGY ==================
67+
v = np.zeros_like(r)
68+
69+
v[:N1] = np.sqrt(mu / r1)
70+
v[N1:N1+Nt] = np.sqrt(mu * (2/r[N1:N1+Nt] - 1/a_t))
71+
v[N1+Nt:] = np.sqrt(mu / r2)
72+
73+
energy = 0.5*v**2 - mu/r
74+
75+
# ================== FIGURE ==================
76+
fig = plt.figure(figsize=(11, 6))
77+
gs = fig.add_gridspec(2, 2)
78+
79+
ax_orbit = fig.add_subplot(gs[:, 0])
80+
ax_speed = fig.add_subplot(gs[0, 1])
81+
ax_energy = fig.add_subplot(gs[1, 1])
82+
83+
# ================== ORBIT VIEW ==================
84+
ax_orbit.set_aspect('equal')
85+
lim = r2 * 1.15
86+
ax_orbit.set_xlim(-lim, lim)
87+
ax_orbit.set_ylim(-lim, lim)
88+
ax_orbit.set_title("Hohmann-Transfer")
89+
90+
ax_orbit.add_patch(plt.Circle((0,0), r_earth, color='skyblue', zorder=0))
91+
ax_orbit.plot(x1, y1, '--', color='gray')
92+
ax_orbit.plot(x_t, y_t, '--', color='orange')
93+
ax_orbit.plot(x2, y2, '--', color='gray')
94+
95+
sat, = ax_orbit.plot([], [], 'ro', markersize=6)
96+
97+
# ================== SPEED GRAPH ==================
98+
ax_speed.set_title("Geschwindigkeit vs Zeit")
99+
ax_speed.set_ylabel("m/s")
100+
ax_speed.set_xlim(0, len(x))
101+
ax_speed.set_ylim(v.min()*0.95, v.max()*1.05)
102+
line_v, = ax_speed.plot([], [], 'r')
103+
104+
# ================== ENERGY GRAPH ==================
105+
ax_energy.set_title("Spezifische Energie")
106+
ax_energy.set_xlabel("Zeitschritt")
107+
ax_energy.set_ylabel("J/kg")
108+
ax_energy.set_xlim(0, len(x))
109+
ax_energy.set_ylim(energy.min()*1.05, energy.max()*0.95)
110+
line_e, = ax_energy.plot([], [], 'b')
111+
112+
# ================== ANIMATION ==================
113+
def update(i):
114+
sat.set_data([x[i]], [y[i]])
115+
line_v.set_data(np.arange(i), v[:i])
116+
line_e.set_data(np.arange(i), energy[:i])
117+
return sat, line_v, line_e
118+
119+
anim = FuncAnimation(
120+
fig,
121+
update,
122+
frames=len(x),
123+
interval=20,
124+
blit=False
125+
)
126+
127+
plt.tight_layout()
128+
plt.show()

0 commit comments

Comments
 (0)