import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#Parámetros problema

gamma = 1.4

p_0 = 101325 #Pa
T_0 = 288.15 #K
rho_0 = 1.225 #Kg/m^3
a_0 = np.sqrt(gamma * p_0 / rho_0)

#ISA
# h siempre en metros h<11000 m

def p_ISA(h):
    return p_0 * (1 - 22.558e-6 * h) ** 5.2559

def T_ISA(h):
    return T_0 * (1 - 22.558e-6 * h)

def rho_ISA(h):
    return rho_0 * (1 - 22.558e-6 * h) ** 4.2559

def a_ISA(h):
    return a_0 * (1 - 22.558e-6 * h) ** 0.5

def dif_presion(V_TAS, h):
    ''' Devuelve la diferencia de presión en Pa que
    mediría el instrumento si fuese perfecto'''
    rho = rho_ISA(h)
    p = p_ISA(h)
    V = V_TAS
    
    y = (gamma - 1) / (2 * gamma) * rho / p * V ** 2 + 1
    y = (y ** ( gamma / (gamma -1)) - 1) * p
    return y

def V_TAS(dif_presion, h):
    p = p_ISA(h)
    rho = rho_ISA(h)
    dp = dif_presion
    
    y =  2 * gamma / (gamma - 1.) * p / rho
    y = y *( (dp/p + 1) ** ((gamma - 1) / gamma) - 1)
    return np.sqrt(y)
    
def V_CAS(h, Vel_TAS):
    dp = dif_presion(Vel_TAS, h)
    return V_TAS(dp, 0)

def V_EAS(h, V_TAS):
    return rho_ISA(h) / rho_0 * V_TAS

h = np.linspace(0, 11000, 1000)
V = np.arange(50, 350, 50)

fig, ax = plt.subplots()

for ii in V:
    ax.plot(h, dif_presion(ii, h)/p_0, label=u'$V_{TAS} = %i m/s$' %ii )

ax.set_ylabel(u'$\Delta p /p_0$')
ax.set_xlabel(u'$h(m)$')
ax.legend()
ax.grid()

h = np.linspace(0, 11000, 1000)
V = np.linspace(0, 340, 1000)

VV, hh = np.meshgrid(V,h)

plt.contourf(VV, hh, dif_presion(VV,hh)/p_0,200)
plt.colorbar()
plt.contour(VV, hh, dif_presion(VV,hh)/p_0, 10, colors = 'black')

plt.title(u'$\Delta p /p_0$')
plt.ylabel(u'$h(m)$')
plt.xlabel(u'$V_{TAS}(m/s)$')

V = 150
h = np.linspace(0, 11000, 1000)

plt.plot(h, V*np.ones_like(h))
plt.plot(h, V_EAS(h, V))
plt.plot(h, V_CAS(h, V))

V = np.arange(50, 400, 50)
h = np.linspace(0, 11000, 1000)

fig, ax = plt.subplots(nrows= np.size(V))
fig.set_size_inches(10,15)

counter = 0

for ii in V:
    axn = ax[counter]
    
    axn.plot(h, ii*np.ones_like(h), label = u'$V_{TAS}$', lw = 2)
    axn.plot(h, V_EAS(h, ii), label = u'$V_{EAS}$')
    axn.plot(h, V_CAS(h, ii), label = u'$V_{CAS}$')
    
    axn.set_ylim(0,ii*1.10)
    axn.set_xlim(0, 11000)
    axn.set_xlabel(u'$h(m)$')
    axn.set_ylabel(u'$V(m/s)$')
    
    axn.legend(loc=3)
    axn.grid()
    
    counter += 1