import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as scin

# Parametres de l'atmosphère et condition initiale en USI
M = 29e-3  # masse molaire moyenne de l'air
g = 9.81 # champ de pesanteur
R = 8.314 # constante des GP
T0 = 288 # température moyenne en z = 0
alpha = 6e-3 # gradient thermique
P0 = 1 # pression en z = 0

# Nombre d’iterations, altitude maximale et pas d'espace en USI
n = 11001
zmax = 11e3
dz = zmax/(n-1)

# Listes des altitudes
z = [i*dz for i in range(n)]
# Création d’une liste pour les pressions obtenues via la méthode d'Euler
P1 = [0 for i in range(n)]
# Définition de la fonction F(P,z)
def F(P,z):
    return(-M*g*P/(R*(T0-alpha*z)))
# Initialisation
P1[0] = P0
# Resolution par la methode d’Euler
for i in range(n-1):
    #P1[i+1] = P1[i]-M*g*dz*P1[i]/(R*(T0-alpha*z[i]))
    P1[i+1] = P1[i]+F(P1[i],z[i])*dz



# Resolution par la fonction odeint
P2=scin.odeint(F,P0,z)

# Representation graphique
plt.plot(z,P1,z,P2)
plt.xlabel("z (m)")
plt.ylabel("P (bar)")
plt.legend(["Euler","odeint"],loc="best")
plt.grid()
plt.show()


# Affichage des differentes distances caracteristiques du probleme
print("distance caract associée au gradient thermique en m:",1/(alpha*T0))
print("dz en m:",dz)