🔬 Amaliyot: Raketa Simulyatsiyasi
Laboratoriya maqsadi
Bu laboratoriyada raketa dinamikasini simulyatsiya qilamiz. Orbital mexanika, trajectory optimization va boshqarish tizimlarini o'rganamiz.
Kerakli kutubxonalar
pip install numpy matplotlib scipy
1-laboratoriya: Raketa Vertikal Ko'tarilishi
Oddiy raketa modeli
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
class SimpleRocket:
"""
Oddiy bir bosqichli raketa modeli
State: [h, v, m]
- h: balandlik (m)
- v: tezlik (m/s)
- m: massa (kg)
"""
def __init__(self):
# Raketa parametrlari
self.m0 = 50000 # Boshlang'ich massa (kg)
self.m_fuel = 40000 # Yoqilg'i massasi (kg)
self.m_dry = 10000 # Quruq massa (kg)
self.m_dot = 300 # Yoqilg'i sarfi (kg/s)
self.Isp = 300 # Solishtirma impuls (s)
self.g0 = 9.81 # Gravitatsiya
self.Cd = 0.3 # Drag koeffitsienti
self.A = 10 # Cross-section area (m²)
# Atmosfera modeli (soddalashtirilgan)
self.rho0 = 1.225 # Sea level havo zichligi
self.H = 8500 # Scale height (m)
# Dvigatel holati
self.engine_on = True
def thrust(self, m):
"""Dvigatel kuchi"""
if self.engine_on and m > self.m_dry:
return self.Isp * self.g0 * self.m_dot
return 0
def air_density(self, h):
"""Balandlikka qarab havo zichligi"""
if h < 0:
return self.rho0
return self.rho0 * np.exp(-h / self.H)
def drag(self, h, v):
"""Aerodinamik qarshilik"""
rho = self.air_density(h)
return 0.5 * rho * self.Cd * self.A * v * abs(v)
def gravity(self, h):
"""Gravitatsiya tezlanishi (Yer radiusini hisobga olgan holda)"""
R_earth = 6371000 # m
return self.g0 * (R_earth / (R_earth + h))**2
def dynamics(self, t, state):
"""Harakat tenglamalari"""
h, v, m = state
# Kuchlar
T = self.thrust(m)
D = self.drag(h, v)
g = self.gravity(h)
# Tezlanish
if m > 0:
a = (T - D) / m - g
else:
a = -g
# Massa o'zgarishi
if self.engine_on and m > self.m_dry:
dm = -self.m_dot
else:
dm = 0
return [v, a, dm]
def simulate_vertical_launch():
"""Vertikal uchirish simulyatsiyasi"""
rocket = SimpleRocket()
# Boshlang'ich holat
state0 = [0, 0, rocket.m0] # h=0, v=0, m=m0
# Event: yoqilg'i tugashi
def fuel_empty(t, state):
return state[2] - rocket.m_dry
fuel_empty.terminal = False
fuel_empty.direction = -1
# Event: yerga qaytish
def hit_ground(t, state):
return state[0]
hit_ground.terminal = True
hit_ground.direction = -1
# Simulyatsiya
sol = solve_ivp(
rocket.dynamics,
[0, 600],
state0,
method='RK45',
events=[fuel_empty, hit_ground],
dense_output=True,
max_step=0.5
)
t = sol.t
h = sol.y[0]
v = sol.y[1]
m = sol.y[2]
# Yoqilg'i tugash vaqti
if len(sol.t_events[0]) > 0:
t_burnout = sol.t_events[0][0]
idx_burnout = np.argmin(np.abs(t - t_burnout))
print(f"Yoqilg'i tugash vaqti: {t_burnout:.1f} s")
print(f"Burnout balandligi: {h[idx_burnout]/1000:.1f} km")
print(f"Burnout tezligi: {v[idx_burnout]:.1f} m/s")
# Max balandlik
h_max = np.max(h)
idx_max = np.argmax(h)
print(f"\nMax balandlik: {h_max/1000:.1f} km")
print(f"Max balandlik vaqti: {t[idx_max]:.1f} s")
print(f"Parvoz vaqti: {t[-1]:.1f} s")
# Vizualizatsiya
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Balandlik
axes[0, 0].plot(t, h/1000, 'b-', linewidth=2)
if len(sol.t_events[0]) > 0:
axes[0, 0].axvline(x=t_burnout, color='r', linestyle='--', label='Burnout')
axes[0, 0].axhline(y=h_max/1000, color='g', linestyle=':', alpha=0.5)
axes[0, 0].set_xlabel('Vaqt (s)')
axes[0, 0].set_ylabel('Balandlik (km)')
axes[0, 0].set_title('Balandlik vs Vaqt', fontsize=14)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# Tezlik
axes[0, 1].plot(t, v, 'r-', linewidth=2)
if len(sol.t_events[0]) > 0:
axes[0, 1].axvline(x=t_burnout, color='r', linestyle='--', label='Burnout')
axes[0, 1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axes[0, 1].set_xlabel('Vaqt (s)')
axes[0, 1].set_ylabel('Tezlik (m/s)')
axes[0, 1].set_title('Tezlik vs Vaqt', fontsize=14)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# Massa
axes[1, 0].plot(t, m/1000, 'g-', linewidth=2)
axes[1, 0].axhline(y=rocket.m_dry/1000, color='r', linestyle='--', label='Quruq massa')
axes[1, 0].set_xlabel('Vaqt (s)')
axes[1, 0].set_ylabel('Massa (ton)')
axes[1, 0].set_title('Massa vs Vaqt', fontsize=14)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# Tezlanish
a = np.gradient(v, t)
axes[1, 1].plot(t, a/rocket.g0, 'm-', linewidth=2)
axes[1, 1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
if len(sol.t_events[0]) > 0:
axes[1, 1].axvline(x=t_burnout, color='r', linestyle='--', label='Burnout')
axes[1, 1].set_xlabel('Vaqt (s)')
axes[1, 1].set_ylabel('Tezlanish (g)')
axes[1, 1].set_title('Tezlanish vs Vaqt', fontsize=14)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
print("="*50)
print("RAKETA VERTIKAL UCHIRISH SIMULYATSIYASI")
print("="*50)
simulate_vertical_launch()
2-laboratoriya: Gravity Turn Maneuvri
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
class GravityTurnRocket:
"""
Gravity turn maneuvrini bajaradigan raketa
State: [r, theta, vr, vtheta, m]
- r: Yer markazidan masofa (m)
- theta: burchak (rad)
- vr: radial tezlik (m/s)
- vtheta: tangensial tezlik (m/s)
- m: massa (kg)
"""
def __init__(self):
# Yer parametrlari
self.R_earth = 6371000 # Yer radiusi (m)
self.mu = 3.986e14 # Gravitatsion parametr (m³/s²)
# Raketa parametrlari
self.m0 = 500000 # Boshlang'ich massa
self.m_fuel = 400000 # Yoqilg'i
self.m_dry = 100000 # Quruq massa
self.m_dot = 2000 # Yoqilg'i sarfi
self.Isp = 350 # Solishtirma impuls
self.g0 = 9.81
# Aerodinamika
self.Cd = 0.2
self.A = 50
self.rho0 = 1.225
self.H = 8500
# Gravity turn parametrlari
self.pitch_rate = 0.002 # rad/s - sekin burilish
self.initial_pitch = np.radians(85) # boshlang'ich pitch
def thrust(self, m):
if m > self.m_dry:
return self.Isp * self.g0 * self.m_dot
return 0
def air_density(self, h):
if h < 0:
return self.rho0
return self.rho0 * np.exp(-h / self.H)
def pitch_angle(self, t, h):
"""Gravity turn pitch burchagi"""
# Vertikal boshlash, sekin burilish
if h < 1000: # Vertikal ko'tarilish
return np.pi/2
# Gravity turn
pitch = self.initial_pitch - self.pitch_rate * t
return max(pitch, 0) # Gorizontaldan pastga tushmasin
def dynamics(self, t, state):
r, theta, vr, vtheta, m = state
h = r - self.R_earth
# Gravitatsiya
g = self.mu / r**2
# Thrust
T = self.thrust(m)
pitch = self.pitch_angle(t, h)
# Thrust components (body frame -> inertial)
Tr = T * np.sin(pitch) / m if m > 0 else 0
Ttheta = T * np.cos(pitch) / m if m > 0 else 0
# Drag
v = np.sqrt(vr**2 + vtheta**2)
rho = self.air_density(h)
D = 0.5 * rho * self.Cd * self.A * v**2
if v > 0 and m > 0:
Dr = -D * vr / (m * v)
Dtheta = -D * vtheta / (m * v)
else:
Dr = Dtheta = 0
# Harakat tenglamalari (polar koordinatalar)
dr = vr
dtheta = vtheta / r
dvr = vtheta**2 / r - g + Tr + Dr
dvtheta = -vr * vtheta / r + Ttheta + Dtheta
dm = -self.m_dot if m > self.m_dry else 0
return [dr, dtheta, dvr, dvtheta, dm]
def simulate_gravity_turn():
"""Gravity turn simulyatsiyasi"""
rocket = GravityTurnRocket()
# Boshlang'ich holat
r0 = rocket.R_earth
theta0 = 0
vr0 = 0
vtheta0 = 0
m0 = rocket.m0
state0 = [r0, theta0, vr0, vtheta0, m0]
# Simulyatsiya
sol = solve_ivp(
rocket.dynamics,
[0, 500],
state0,
method='RK45',
dense_output=True,
max_step=1.0
)
t = sol.t
r = sol.y[0]
theta = sol.y[1]
vr = sol.y[2]
vtheta = sol.y[3]
m = sol.y[4]
h = r - rocket.R_earth
v = np.sqrt(vr**2 + vtheta**2)
# Orbital parametrlar
# Orbital energy
E = 0.5 * v**2 - rocket.mu / r
# Semi-major axis
a = -rocket.mu / (2 * E)
# Eccentricity vector
h_vec = r * vtheta # specific angular momentum
e = np.sqrt(1 + 2 * E * h_vec**2 / rocket.mu**2)
print(f"Final balandlik: {h[-1]/1000:.1f} km")
print(f"Final tezlik: {v[-1]:.1f} m/s")
print(f"Orbital tezlik (bu balandlikda): {np.sqrt(rocket.mu/r[-1]):.1f} m/s")
print(f"\nOrbital parametrlar:")
print(f"Semi-major axis: {a[-1]/1000:.1f} km")
print(f"Eccentricity: {e[-1]:.4f}")
# Vizualizatsiya
fig = plt.figure(figsize=(16, 10))
# Trajectory
ax1 = fig.add_subplot(2, 2, 1, projection='polar')
ax1.plot(theta, r/1000, 'b-', linewidth=2)
# Yer
earth_theta = np.linspace(0, 2*np.pi, 100)
ax1.fill(earth_theta, np.ones(100) * rocket.R_earth/1000, color='lightblue', alpha=0.7)
ax1.set_title('Trajectory (Polar)', fontsize=14)
ax1.set_rlabel_position(45)
# Balandlik va tezlik
ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(t, h/1000, 'b-', linewidth=2, label='Balandlik')
ax2.set_xlabel('Vaqt (s)')
ax2.set_ylabel('Balandlik (km)', color='b')
ax2.tick_params(axis='y', labelcolor='b')
ax2.grid(True, alpha=0.3)
ax2b = ax2.twinx()
ax2b.plot(t, v/1000, 'r-', linewidth=2, label='Tezlik')
ax2b.set_ylabel('Tezlik (km/s)', color='r')
ax2b.tick_params(axis='y', labelcolor='r')
ax2.set_title('Balandlik va Tezlik', fontsize=14)
# Downrange vs Altitude
ax3 = fig.add_subplot(2, 2, 3)
downrange = rocket.R_earth * theta
ax3.plot(downrange/1000, h/1000, 'g-', linewidth=2)
ax3.set_xlabel('Downrange (km)')
ax3.set_ylabel('Balandlik (km)')
ax3.set_title('Trajectory Profile', fontsize=14)
ax3.grid(True, alpha=0.3)
# Pitch angle
ax4 = fig.add_subplot(2, 2, 4)
pitch = [rocket.pitch_angle(ti, hi) for ti, hi in zip(t, h)]
ax4.plot(t, np.degrees(pitch), 'm-', linewidth=2)
ax4.set_xlabel('Vaqt (s)')
ax4.set_ylabel('Pitch burchagi (°)')
ax4.set_title('Gravity Turn Profili', fontsize=14)
ax4.grid(True, alpha=0.3)
ax4.axhline(y=0, color='k', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
print("="*50)
print("GRAVITY TURN SIMULYATSIYASI")
print("="*50)
simulate_gravity_turn()
3-laboratoriya: Orbital Mexanika
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
class OrbitalMechanics:
"""Orbital mexanika simulyatsiyasi"""
def __init__(self):
# Yer parametrlari
self.R_earth = 6371000 # m
self.mu = 3.986e14 # m³/s²
def orbital_elements_to_state(self, a, e, i, omega, Omega, nu):
"""
Orbital elementlardan pozitsiya va tezlikka
a: semi-major axis (m)
e: eccentricity
i: inclination (rad)
omega: argument of periapsis (rad)
Omega: RAAN (rad)
nu: true anomaly (rad)
"""
# Orbital plane da
p = a * (1 - e**2) # semi-latus rectum
r = p / (1 + e * np.cos(nu))
# Pozitsiya (orbital plane)
r_orbital = np.array([r * np.cos(nu), r * np.sin(nu), 0])
# Tezlik (orbital plane)
h = np.sqrt(self.mu * p)
v_orbital = np.array([
-self.mu / h * np.sin(nu),
self.mu / h * (e + np.cos(nu)),
0
])
# Rotation matrices
R3_Omega = self._rot_z(-Omega)
R1_i = self._rot_x(-i)
R3_omega = self._rot_z(-omega)
Q = R3_Omega @ R1_i @ R3_omega
r_eci = Q @ r_orbital
v_eci = Q @ v_orbital
return r_eci, v_eci
def _rot_x(self, angle):
c, s = np.cos(angle), np.sin(angle)
return np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
def _rot_z(self, angle):
c, s = np.cos(angle), np.sin(angle)
return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
def two_body_dynamics(self, t, state):
"""Ikki jism muammosi"""
r = state[:3]
v = state[3:]
r_mag = np.linalg.norm(r)
a = -self.mu / r_mag**3 * r
return np.concatenate([v, a])
def hohmann_transfer(self, r1, r2):
"""
Hohmann transfer hisoblash
r1: boshlang'ich orbita radiusi
r2: yakuniy orbita radiusi
"""
# Transfer orbita
a_t = (r1 + r2) / 2
# Tezlik o'zgarishlari
v1 = np.sqrt(self.mu / r1) # Circular velocity at r1
v2 = np.sqrt(self.mu / r2) # Circular velocity at r2
v_t1 = np.sqrt(self.mu * (2/r1 - 1/a_t)) # Transfer velocity at r1
v_t2 = np.sqrt(self.mu * (2/r2 - 1/a_t)) # Transfer velocity at r2
delta_v1 = v_t1 - v1
delta_v2 = v2 - v_t2
# Transfer vaqti
T_transfer = np.pi * np.sqrt(a_t**3 / self.mu)
return delta_v1, delta_v2, T_transfer
def simulate_orbit():
"""Orbit simulyatsiyasi"""
orb = OrbitalMechanics()
# ISS ga o'xshash orbita
a = orb.R_earth + 420000 # 420 km altitude
e = 0.001 # Nearly circular
i = np.radians(51.6) # ISS inclination
omega = 0
Omega = 0
nu = 0
r0, v0 = orb.orbital_elements_to_state(a, e, i, omega, Omega, nu)
state0 = np.concatenate([r0, v0])
# Simulyatsiya (bir nechta orbit)
T_orbit = 2 * np.pi * np.sqrt(a**3 / orb.mu)
t_span = [0, 3 * T_orbit]
sol = solve_ivp(
orb.two_body_dynamics,
t_span,
state0,
method='RK45',
dense_output=True,
max_step=60
)
t = sol.t
r = sol.y[:3]
v = sol.y[3:]
# 3D vizualizatsiya
fig = plt.figure(figsize=(14, 6))
ax1 = fig.add_subplot(121, projection='3d')
# Orbita
ax1.plot(r[0]/1e6, r[1]/1e6, r[2]/1e6, 'b-', linewidth=1, label='Orbit')
# Yer
u = np.linspace(0, 2*np.pi, 50)
v_sphere = np.linspace(0, np.pi, 50)
x_earth = orb.R_earth/1e6 * np.outer(np.cos(u), np.sin(v_sphere))
y_earth = orb.R_earth/1e6 * np.outer(np.sin(u), np.sin(v_sphere))
z_earth = orb.R_earth/1e6 * np.outer(np.ones(np.size(u)), np.cos(v_sphere))
ax1.plot_surface(x_earth, y_earth, z_earth, color='lightblue', alpha=0.6)
ax1.set_xlabel('X (1000 km)')
ax1.set_ylabel('Y (1000 km)')
ax1.set_zlabel('Z (1000 km)')
ax1.set_title('ISS Orbit', fontsize=14)
ax1.legend()
# Balandlik vs vaqt
ax2 = fig.add_subplot(122)
altitude = np.linalg.norm(r, axis=0) - orb.R_earth
ax2.plot(t/3600, altitude/1000, 'b-', linewidth=1)
ax2.set_xlabel('Vaqt (soat)')
ax2.set_ylabel('Balandlik (km)')
ax2.set_title('Orbital Altitude', fontsize=14)
ax2.grid(True, alpha=0.3)
# Orbital period
print(f"Orbital period: {T_orbit/60:.1f} minutes")
print(f"Altitude: {(a - orb.R_earth)/1000:.1f} km")
print(f"Velocity: {np.linalg.norm(v0):.1f} m/s")
plt.tight_layout()
plt.show()
def simulate_hohmann():
"""Hohmann transfer simulyatsiyasi"""
orb = OrbitalMechanics()
# LEO dan GEO ga transfer
r1 = orb.R_earth + 400000 # LEO: 400 km
r2 = orb.R_earth + 35786000 # GEO: 35786 km
dv1, dv2, T_transfer = orb.hohmann_transfer(r1, r2)
print("="*50)
print("HOHMANN TRANSFER: LEO -> GEO")
print("="*50)
print(f"LEO radius: {r1/1000:.0f} km")
print(f"GEO radius: {r2/1000:.0f} km")
print(f"\nDelta-V1 (LEO da): {dv1:.1f} m/s")
print(f"Delta-V2 (GEO da): {dv2:.1f} m/s")
print(f"Total Delta-V: {dv1 + dv2:.1f} m/s")
print(f"Transfer vaqti: {T_transfer/3600:.1f} soat")
# Simulyatsiya
# Boshlang'ich LEO
state0 = np.array([r1, 0, 0, 0, np.sqrt(orb.mu/r1), 0])
# Phase 1: LEO
T_leo = 2 * np.pi * np.sqrt(r1**3 / orb.mu)
sol1 = solve_ivp(orb.two_body_dynamics, [0, T_leo], state0,
method='RK45', max_step=60)
# Burn 1
state_burn1 = np.array([r1, 0, 0, 0, np.sqrt(orb.mu/r1) + dv1, 0])
# Phase 2: Transfer
sol2 = solve_ivp(orb.two_body_dynamics, [0, T_transfer], state_burn1,
method='RK45', max_step=300)
# Burn 2
state_burn2 = sol2.y[:, -1].copy()
state_burn2[4] += dv2 # vy ga qo'shish
# Phase 3: GEO
T_geo = 2 * np.pi * np.sqrt(r2**3 / orb.mu)
sol3 = solve_ivp(orb.two_body_dynamics, [0, T_geo], state_burn2,
method='RK45', max_step=300)
# Vizualizatsiya
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Orbitalar
ax1 = axes[0]
# LEO
theta = np.linspace(0, 2*np.pi, 100)
x_leo = r1 * np.cos(theta) / 1e6
y_leo = r1 * np.sin(theta) / 1e6
ax1.plot(x_leo, y_leo, 'b--', linewidth=1, label='LEO', alpha=0.5)
# GEO
x_geo = r2 * np.cos(theta) / 1e6
y_geo = r2 * np.sin(theta) / 1e6
ax1.plot(x_geo, y_geo, 'g--', linewidth=1, label='GEO', alpha=0.5)
# Transfer
ax1.plot(sol2.y[0]/1e6, sol2.y[1]/1e6, 'r-', linewidth=2, label='Transfer')
# Yer
earth = plt.Circle((0, 0), orb.R_earth/1e6, color='lightblue')
ax1.add_patch(earth)
# Burn points
ax1.scatter(r1/1e6, 0, c='red', s=100, marker='*', zorder=5, label='Burn 1')
ax1.scatter(sol2.y[0, -1]/1e6, sol2.y[1, -1]/1e6, c='green', s=100, marker='*', zorder=5, label='Burn 2')
ax1.set_xlabel('X (1000 km)')
ax1.set_ylabel('Y (1000 km)')
ax1.set_title('Hohmann Transfer: LEO → GEO', fontsize=14)
ax1.set_aspect('equal')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
# Delta-V bar chart
ax2 = axes[1]
labels = ['ΔV₁ (LEO)', 'ΔV₂ (GEO)', 'Total']
values = [dv1, dv2, dv1 + dv2]
colors = ['red', 'green', 'blue']
bars = ax2.bar(labels, values, color=colors, alpha=0.7)
ax2.set_ylabel('Delta-V (m/s)')
ax2.set_title('Hohmann Transfer Delta-V', fontsize=14)
ax2.grid(True, alpha=0.3, axis='y')
for bar, val in zip(bars, values):
ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 50,
f'{val:.0f}', ha='center', va='bottom', fontsize=12)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
simulate_orbit()
simulate_hohmann()
4-laboratoriya: Raketa Boshqarish Tizimi
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
class RocketGuidance:
"""Raketa yo'naltirish va boshqarish tizimi"""
def __init__(self):
# Raketa parametrlari
self.m0 = 100000
self.m_fuel = 80000
self.m_dry = 20000
self.Isp = 320
self.g0 = 9.81
self.max_thrust = 1500000 # N
self.throttle = 1.0
# Boshqarish parametrlari
self.Kp_pitch = 0.5
self.Kd_pitch = 0.3
def thrust_magnitude(self, m):
if m > self.m_dry:
return self.max_thrust * self.throttle
return 0
def mass_flow_rate(self, m):
if m > self.m_dry:
T = self.thrust_magnitude(m)
return T / (self.Isp * self.g0)
return 0
def guidance_law(self, t, state, target_pitch):
"""
Proportional Navigation guidance
state: [x, y, vx, vy, pitch, pitch_rate, m]
"""
x, y, vx, vy, pitch, pitch_rate, m = state
# Pitch error
pitch_error = target_pitch - pitch
# PD control
pitch_command = self.Kp_pitch * pitch_error - self.Kd_pitch * pitch_rate
# Rate limit
max_rate = np.radians(5) # 5 deg/s
pitch_command = np.clip(pitch_command, -max_rate, max_rate)
return pitch_command
def dynamics(self, t, state, target_pitch_func):
x, y, vx, vy, pitch, pitch_rate, m = state
# Boshqarish
target_pitch = target_pitch_func(t, x, y)
pitch_cmd = self.guidance_law(t, state, target_pitch)
# Kuchlar
T = self.thrust_magnitude(m)
g = 9.81 * (6371000 / (6371000 + y))**2
ax = T * np.cos(pitch) / m
ay = T * np.sin(pitch) / m - g
# Mass
dm = -self.mass_flow_rate(m)
return [vx, vy, ax, ay, pitch_rate, pitch_cmd * 10, dm]
def simulate_guided_launch():
"""Boshqariladigan uchirish simulyatsiyasi"""
rocket = RocketGuidance()
def target_pitch(t, x, y):
"""Gravity turn pitch profili"""
if y < 1000:
return np.radians(89)
elif y < 10000:
return np.radians(80)
elif y < 30000:
return np.radians(60)
elif y < 60000:
return np.radians(40)
else:
return np.radians(20)
# Boshlang'ich holat
state0 = [0, 0, 0, 0, np.radians(89), 0, rocket.m0]
def dynamics_wrapper(t, state):
return rocket.dynamics(t, state, target_pitch)
# Simulyatsiya
sol = solve_ivp(
dynamics_wrapper,
[0, 200],
state0,
method='RK45',
max_step=0.5
)
t = sol.t
x, y, vx, vy, pitch, pitch_rate, m = sol.y
v = np.sqrt(vx**2 + vy**2)
flight_path = np.arctan2(vy, vx)
# Vizualizatsiya
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Trajectory
axes[0, 0].plot(x/1000, y/1000, 'b-', linewidth=2)
axes[0, 0].scatter(0, 0, c='green', s=100, marker='o', label='Launch')
axes[0, 0].scatter(x[-1]/1000, y[-1]/1000, c='red', s=100, marker='*', label='Burnout')
axes[0, 0].set_xlabel('Downrange (km)')
axes[0, 0].set_ylabel('Altitude (km)')
axes[0, 0].set_title('Trajectory', fontsize=14)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# Pitch vs Flight path
axes[0, 1].plot(t, np.degrees(pitch), 'b-', linewidth=2, label='Pitch (thrust direction)')
axes[0, 1].plot(t, np.degrees(flight_path), 'r--', linewidth=2, label='Flight path angle')
axes[0, 1].set_xlabel('Vaqt (s)')
axes[0, 1].set_ylabel('Burchak (°)')
axes[0, 1].set_title('Pitch Control', fontsize=14)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# Velocity
axes[1, 0].plot(t, v, 'g-', linewidth=2)
axes[1, 0].set_xlabel('Vaqt (s)')
axes[1, 0].set_ylabel('Tezlik (m/s)')
axes[1, 0].set_title('Tezlik', fontsize=14)
axes[1, 0].grid(True, alpha=0.3)
# Mass
axes[1, 1].plot(t, m/1000, 'm-', linewidth=2)
axes[1, 1].axhline(y=rocket.m_dry/1000, color='r', linestyle='--', label='Dry mass')
axes[1, 1].set_xlabel('Vaqt (s)')
axes[1, 1].set_ylabel('Massa (ton)')
axes[1, 1].set_title('Massa', fontsize=14)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Natijalar
print("="*50)
print("BOSHQARILADIGAN UCHIRISH NATIJALARI")
print("="*50)
print(f"Final altitude: {y[-1]/1000:.1f} km")
print(f"Final downrange: {x[-1]/1000:.1f} km")
print(f"Final velocity: {v[-1]:.1f} m/s")
print(f"Flight path angle: {np.degrees(flight_path[-1]):.1f}°")
if __name__ == "__main__":
simulate_guided_launch()
Xulosa
Bu laboratoriyalarda biz:
- ✅ Raketa vertikal uchirish simulyatsiyasini yaratdik
- ✅ Gravity turn manevrini amalga oshirdik
- ✅ Orbital mexanikani o'rgandik
- ✅ Hohmann transfer ni simulyatsiya qildik
- ✅ Raketa boshqarish tizimini yaratdik
Keyingi qadamlar
- Staging (ko'p bosqichli raketalar)
- Reentry va thermal analysis
- Landing guidance (SpaceX Falcon 9)
- Interplanetary transfers