Model Rocket Simulator

import matplotlib.pyplot as plt
from numpy import sign

# Motors

EstesA8  = { "t1": 0.22 , "t2": 0.30 , "t3": 0.65 , "t4": 0.72 , "fmax": 10.0 , "fsus":  2.2 }
EstesB6  = { "t1": 0.18 , "t2": 0.30 , "t3": 0.80 , "t4": 0.86 , "fmax": 12.0 , "fsus":  4.5 }
EstesC6  = { "t1": 0.18 , "t2": 0.24 , "t3": 1.80 , "t4": 1.84 , "fmax": 14.6 , "fsus":  4.5 }
EstesD12 = { "t1": 0.35 , "t2": 0.40 , "t3": 1.60 , "t4": 1.70 , "fmax": 30.0 , "fsus": 10.0 }
KlimaA6  = { "t1": 0.35 , "t2": 0.70 , "t3": 0.75 , "t4": 0.75 , "fmax":  7.0 , "fsus":  0.0 }
KlimaB4  = { "t1": 0.55 , "t2": 0.70 , "t3": 1.00 , "t4": 1.25 , "fmax":  8.0 , "fsus":  4.0 }
KlimaC2  = { "t1": 0.55 , "t2": 0.70 , "t3": 4.80 , "t4": 5.05 , "fmax":  5.0 , "fsus":  2.0 }
KlimaD3  = { "t1": 0.65 , "t2": 1.00 , "t3": 5.50 , "t4": 6.00 , "fmax":  9.0 , "fsus":  3.0 }
KlimaD9  = { "t1": 0.25 , "t2": 0.60 , "t3": 1.80 , "t4": 2.15 , "fmax": 25.0 , "fsus":  9.0 }

# Environment
g   = 9.8 # m/s^2
rho = 1.2 # kg/m^3
dt  = 0.001 # s
rod = 1 # m - rod length

# Rocket
diam  = 0.042 # m
m     = 0.100  # kg
cd    = 0.65
chute = 0.35  # m

motor = EstesC6
delay = 5

def getThrust(t):
    if t<=motor["t1"]:
        thrust = motor["fmax"]*t/motor["t1"]
    if motor["t1"]<t<motor["t2"]:
        thrust = motor["fmax"] + ( (t-motor["t1"])/(motor["t2"]-motor["t1"]) * (motor["fsus"]-motor["fmax"]) )
    if motor["t2"]<=t<=motor["t3"]:
        thrust = motor["fsus"]
    if motor["t3"]<t<=motor["t4"]:
        thrust = motor["fsus"] - ( (t-motor["t3"])/(motor["t4"]-motor["t3"]) * motor["fsus"] )
    if t>motor["t4"]:
        thrust=0
    return thrust

h=0
v=0
i=0
t=0
hold=0
tmax=100

liftoff = False
rodclr  = False
burnout = False
apogee  = False
deploy  = False
landing = False

lt  = []
lh  = []
lv  = []
la  = []
lq  = []
lFr = []
lFd = []
lFg = []
lFtot = []

while t<tmax and h>=0:
    Fr = getThrust(t)

    Fg = -m*g

    q  = 0.5 * rho * v**2
    Fd = q * 3.14159 * (diam/2)**2 * cd * -sign(v)

    Ftot = Fr + Fg + Fd
    a = Ftot/m

    v = v + a*dt
    h = h + v*dt

    if h<0 and not liftoff:
        h=0
        v=0

    if h>0 and not liftoff:
        h=0
        print("Liftoff at T+ %.2f" % t + " s")
        liftoff=True

    if h>rod and not rodclr:
        print("Rod clr at T+ %.2f" % t + " s, v=%.2f" % abs(v) + " m/s")
        rodclr=True

    if t>motor["t4"] and not burnout:
        print("Burnout at T+ %.2f" % t + " s, h=%.2f" % h + " m, v=%.2f" % v + " m/s")
        burnout=True

    if h<hold and not apogee:
        print("Apogee  at T+ %.2f" % t + " s, h=%.2f" % h + " m")
        apogee=True
    hold = h

    if t>motor["t4"]+delay and not deploy:
        print("Deploy  at T+ %.2f" % t + " s, h=%.2f" % h + " m, v=%.2f" % abs(v) + " m/s")
        diam = chute
        deploy=True
        if not apogee:
            print("Deployed before apogee.")

    if h<0 and apogee:
        print("Landing at T+%.2f" % t + " s, v=%.2f" % abs(v) + " m/s")
        if not deploy:
            print("Landed before deployment.")

    lt.append(t)
    lh.append(h)
    lv.append(v)
    la.append(a/g)
    lq.append(q)
    lFr.append(Fr)
    lFd.append(Fd)
    lFg.append(Fg)
    lFtot.append(Ftot)

    i=i+1
    t=t+dt

print()

plt.subplot(411)
plt.plot(lt,lh)
plt.ylabel('Altitude / [m]')
plt.grid(True)

plt.subplot(412)
plt.plot(lt,lv)
plt.ylabel('Speed / [m/s]')
plt.grid(True)

plt.subplot(413)
plt.plot(lt,la)
plt.ylabel('Acceleration / [G]')
plt.grid(True)

plt.subplot(414)
plt.plot(lt,lFr)
plt.xlabel('Time / [s]')
plt.ylabel('Dynamic Pressure / [Pa]')
plt.grid(True)

plt.show()