------------
import numpy as np
import matplotlib.pyplot as plt
# --- FIZIKAI ÉS PLANETÁRIS KONSTANSOK ---
G = 6.67430e-11 # Gravitációs állandó (m^3 kg^-1 s^-2)
M_EARTH = 5.972e24 # Föld tömege (kg)
R_EARTH = 6371000 # Föld sugara (m)
# --- LÉGKÖRI MODELL PARAMÉTEREK ---
RHO_0 = 1.225 # Tengerszinti levegősűrűség (kg/m^3)
H_SCALE = 8500 # Skálamagasság (m) - kb. 8.5 km-enként feleződik a sűrűség
# --- ASZTEROIDA PARAMÉTEREK (2 km átmérő) ---
DIAMETER = 2000 # Átmérő (m)
RADIUS_ASTEROID = DIAMETER / 2
CROSS_SECTION = np.pi * (RADIUS_ASTEROID**2) # Homlokfelület (m^2)
DENSITY_ASTEROID = 3000 # Átlagos kőzet-aszteroida sűrűség (kg/m^3)
MASS_ASTEROID = (4/3) * np.pi * (RADIUS_ASTEROID**3) * DENSITY_ASTEROID # Tömeg (kg)
C_D = 1.2 # Légellenállási alaktényező (gömb/szabálytalan szikla közelítése)
# --- KEZDETI FELTÉTELEK ---
h_perigee = 300000 # 300 km magasság a perigeumban (m)
r_perigee = R_EARTH + h_perigee
# Sebesség meghatározása (hiperbolikus pálya, ~14 km/s a perigeumban)
v_escape = np.sqrt(2 * G * M_EARTH / r_perigee)
v_perigee = v_escape + 3000
# Mivel a légellenállás nem konzervatív erő (veszteséges), a pályát eltoljuk,
# és egy távoli pontból indítjuk a légkör felett (pl. 2500 km magasan).
h_start = 2500000 # 2500 km kezdőmagasság
r_start = R_EARTH + h_start
# Hiperbolikus pályaenergia (E) és impulzusmomentum (L) kiszámítása a kívánt perigeumhoz
E = 0.5 * v_perigee**2 - (G * M_EARTH / r_perigee)
L = r_perigee * v_perigee
# Kezdőpozíció szöge a hiperbola egyenletéből
p = L**2 / (G * M_EARTH)
e = np.sqrt(1 + (2 * E * L**2) / (G * M_EARTH)**2)
cos_theta = (p / r_start - 1) / e
theta_start = -np.arccos(cos_theta) # negatív szög, mert közeledik
# Kezdő koordináták és sebességek polár koordinátákból forgatva
x0 = r_start * np.cos(theta_start)
y0 = r_start * np.sin(theta_start)
v_radial = (G * M_EARTH * e * np.sin(theta_start)) / L
v_tangential = L / r_start
vx0 = v_radial * np.cos(theta_start) - v_tangential * np.sin(theta_start)
vy0 = v_radial * np.sin(theta_start) + v_tangential * np.cos(theta_start)
state_0 = np.array([x0, y0, vx0, vy0])
# --- DIFFERENCIÁLEGYENLET LÉGELLENÁLLÁSSAL ---
def gravity_and_drag_derivs(t, state):
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
h = r - R_EARTH
# 1. Gravitációs gyorsulás
ax_g = -G * M_EARTH * x / r**3
ay_g = -G * M_EARTH * y / r**3
# 2. Légköri sűrűség kiszámítása (ha h < 0, a szikla becsapódott)
if h > 0:
rho = RHO_0 * np.exp(-h / H_SCALE)
else:
rho = 0
# 3. Légellenállási erő (Drag) és gyorsulás
v_mag = np.sqrt(vx**2 + vy**2)
drag_force = 0.5 * rho * v_mag**2 * C_D * CROSS_SECTION
a_drag = drag_force / MASS_ASTEROID
# A légellenállás a sebességvektorral ellentétes irányú
ax_d = -a_drag * (vx / v_mag) if v_mag > 0 else 0
ay_d = -a_drag * (vy / v_mag) if v_mag > 0 else 0
return np.array([vx, vy, ax_g + ax_d, ay_g + ay_d])
# RK4 Integrátor lépés
def rk4_step(f, t, state, dt):
k1 = f(t, state)
k2 = f(t + dt/2, state + k1 * dt/2)
k3 = f(t + dt/2, state + k2 * dt/2)
k4 = f(t + dt, state + k3 * dt)
return state + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4)
# --- SZIMULÁCIÓ FUTTATÁSA ---
dt = 0.5 # Időlépés másodpercben
current_state = state_0.copy()
trajectory = [current_state]
max_steps = 16000
for _ in range(max_steps):
current_state = rk4_step(gravity_and_drag_derivs, 0, current_state, dt)
trajectory.append(current_state)
# Ellenőrzések az aktuális lépésben
r_curr = np.sqrt(current_state[0]**2 + current_state[1]**2)
if r_curr < R_EARTH:
print("A szikla becsapódott a Földbe!")
break
if r_curr > r_start and len(trajectory) > 1000:
break
trajectory = np.array(trajectory)
x_coords = trajectory[:, 0] / 1000
y_coords = trajectory[:, 1] / 1000
# Energiaveszteség kiszámítása
v_start_mag = np.sqrt(vx0**2 + vy0**2)
v_end_mag = np.sqrt(trajectory[-1, 2]**2 + trajectory[-1, 3]**2)
energy_loss_pct = (1 - (v_end_mag**2 / v_start_mag**2)) * 100
# --- VIZUALIZÁCIÓ ---
fig, ax = plt.subplots(figsize=(9, 9))
# Föld és légkör rajzolása
earth_circle = plt.Circle((0, 0), R_EARTH / 1000, color='royalblue', label='Föld felszín')
ax.add_patch(earth_circle)
# Légköri rétegek szemléltetése
ax.add_patch(plt.Circle((0, 0), (R_EARTH + 100000) / 1000, color='deepskyblue', alpha=0.15, label='Sűrűbb légkör (<100 km)'))
ax.add_patch(plt.Circle((0, 0), (R_EARTH + 300000) / 1000, color='skyblue', alpha=0.08, label='Ritka exoszféra (<300 km)'))
# Pálya kirajzolása
ax.plot(x_coords, y_coords, color='orangered', linewidth=2, label='Aszteroida pálya (RK4 + Drag)')
# Formázás
ax.set_aspect('equal', 'box')
ax.set_xlabel('X pozíció (km)')
ax.set_ylabel('Y pozíció (km)')
ax.set_title('Earth-grazing szimuláció légellenállással (Drag Model)')
ax.grid(True, linestyle=':', alpha=0.5)
ax.legend(loc='upper right')
# Fókuszálás a Földközelségre
ax.set_xlim(-8000, 8000)
ax.set_ylim(-4000, 8000)
# Információs panel
info_text = (
f"Szikla átmérő: {DIAMETER/1000:.1f} km\n"
f"Szikla tömeg: {MASS_ASTEROID:.2e} kg\n"
f"Minimális magasság: {h_perigee/1000:.0f} km\n"
f"Sebességveszteség: {v_start_mag - v_end_mag:.6f} m/s\n"
f"Kinetikus energia csökkenés: {energy_loss_pct:.8f}%"
)
plt.gca().text(0.03, 0.03, info_text, transform=plt.gca().transAxes,
bbox=dict(facecolor='white', alpha=0.85, edgecolor='gray'))
plt.show()
-----------
-------------
---------
Nincsenek megjegyzések:
Megjegyzés küldése