R2A Landing Velocity Estimation

main.py

""" 
`pip install numpy scipy matplotlib`
`pip install -e git+https://github.com/AeroPython/scikit-aero#egg=scikit-aero`

Author: @rubendibattista
"""

# import matplotlib
# matplotlib.use('QT5Agg')  # Needed for OSX with virtualenv
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from skaero.atmosphere import coesa

mass = 50
CDx = 0.5
CDy = 0.4
g = 9.81
D = 175E-3

# === Cross-flow surfaces ===
S = np.pi*D**2/4
H = 5  # Length of the rocket
S_lat = D*H
surf_ratio = 0.9  # To take into account inclination

Sx = surf_ratio*S_lat + (1-surf_ratio)*S
Sy = surf_ratio*S + (1-surf_ratio)*S_lat

print(Sx)
print(Sy)


def dX(X, t):
    _ = X[0]
    vx = X[1]
    y = X[2]
    vy = X[3]

    _, _, _, rho = coesa.table(y)


    Dx = 0.5*rho*vx**2*Sx*CDx
    Dy = 0.5*rho*vy**2*Sy*CDy

    dx = vx
    dvx = - Dx/mass

    dy = vy
    dvy = -g + Dy/mass

    return [dx, dvx, dy, dvy]


apogee = 10E3  # Apogee

# Initial Conditions
x0 = 0
vx0 = 100
y0 = apogee
vy0 = 0

X0 = [
    x0,
    vx0,
    y0,
    vy0
]

t = np.linspace(0, 130, 100)

X = odeint(dX, X0, t)

# plt.plot(t, X[:, 0], label='x')
plt.plot(t, X[:, 2], label='y')
plt.xlabel('t[s]')
plt.ylabel('y[m]')
plt.legend(loc='best')
plt.tight_layout()

plt.figure()
plt.plot(t, X[:, 3], label=r'$v_y$')
plt.plot(t, X[:, 1], label=r'$v_x$')
plt.plot(t, np.sqrt(X[:, 1]**2 + X[:, 3]**2), label=r'$|v|$')
plt.xlabel('t[s]')
plt.ylabel(r'Velocity $[m/s]$')
plt.legend(loc='best')
plt.tight_layout()


plt.show()