simulate launch of a rocket with two stages in python

main.py
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
G = 6.674e-11    # gravitational constant
M_e = 5.972e24   # mass of Earth
R_e = 6.371e6    # radius of Earth
M_s1 = 2e5       # mass of stage 1
M_s2 = 1e5       # mass of stage 2
v_e = 7.9e3      # exhaust velocity
A = 25           # cross-sectional area
Cd = 0.6         # drag coefficient
t_burn1 = 100    # burn time of stage 1
t_burn2 = 50     # burn time of stage 2
h0 = 0           # initial height
v0 = 0           # initial velocity

# Equations of motion
def rocket(t, y):
    """
    Differential equations for rocket motion.
    y = [h, v, m]
    """
    h, v, m = y
    if t <= t_burn1:
        dm_dt = -M_s1 / t_burn1
        F_thrust = v_e * dm_dt
    elif t <= t_burn1 + t_burn2:
        dm_dt = -M_s2 / t_burn2
        F_thrust = v_e * dm_dt
    else:
        dm_dt = 0
        F_thrust = 0
    F_grav = -G * M_e * m / (R_e + h)**2
    F_drag = -0.5 * Cd * A * v**2
    dh_dt = v
    dv_dt = (F_thrust + F_grav + F_drag) / m
    return [dh_dt, dv_dt, dm_dt]

# Initial conditions
y0 = [h0, v0, M_s1 + M_s2]

# Solve differential equations
sol = solve_ivp(rocket, [0, 500], y0, dense_output=True)

# Plot results
t = np.linspace(0, 500, 1000)
y = sol.sol(t)
plt.plot(t, y[0])
plt.xlabel('Time (s)')
plt.ylabel('Height (m)')
plt.show()
1353 chars
54 lines