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