#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
PRESSION ATMOSPHERIQUE DANS LE MODELE ISA
Created on Sun Apr  4 11:45:19 2021

@author: Emilie Fremont
"""

## Importation des bibliotheques
## ------------------------------------

import numpy as np                  # pour la manipulation des tableaux
import matplotlib.pyplot as plt     # pour les representations graphiques
from scipy.integrate import odeint  # pour la resolution des equations differentielles


## Definition des constantes du probleme
## -------------------------------------

g = 9.81                            # acceleration de la pesanteur (en m/s^2)
Mair = 29e-3                        # masse molaire de l'air (en kg/mol)
R = 8.314                           # constante du gaz parfait (en J/K/mol)
Tsol = 288                          # temperature de l'atmosphere au niveau du sol (en K)
Psol = 1.013e5                      # pression de l'atmosphere au niveau du sol (en Pa)


## Definition du gradient thermique vertical selon le modèle ISA
## -------------------------------------------------------------

def kISA(z):
    """ z est l'altitude en metres. La fonction renvoie la valeur du gradient thermique
    vertical à l'altitude z (en K/m). """
    if 0 <= z < 11e3: return -6.5e-3
    elif z < 20e3: return 0
    elif z < 32e3: return 1.0e-3
    elif z < 47e3: return 2.8e-3
    elif z < 51e3: return 0
    elif z < 71e3: return -2.8e-3
    elif z < 85e3: return -2.0e-3
    else: return None
    
    
## -----------------------------------------------------------------------------------------------
## A - Resolution numerique
## -----------------------------------------------------------------------------------------------

## Definition du systeme differentiel a resoudre
    
def systDiff(TP,z):
    """ 
    TP désigne le vecteur inconnu de dimension 2 (TP[0] : temperature ; TP[1] : pression) ;
    z désigne l'altitude. 
    La fonction renvoie respectivement la derivee de la temperature et la dérivee de
    la pression a l'altitude z.
    """
    # Lois prevues par le modele theorique
    dT = kISA(z)                    # derivee verticale de la temperature
    dP = - Mair*g/R*TP[1]/TP[0]     # derivee verticale de la pression
    
    return [dT, dP]

## Definition des conditions aux limites

CAL = [Tsol, Psol]

## Definition de l'ensemble des valeurs de z pour lesquelles on cherche la solution
## numerique approchee du systeme differentiel 

z = np.linspace(0, 85e3, 10000)     # on choisit 10000 points régulierement espaces entre 0
                                    # et 85 km d'altitude

## La resolution en elle-meme
                                    
TP = odeint(systDiff, CAL, z)
T = TP[:,0]                         # extraction des valeurs de la temperature
P = TP[:,1]                         # extraction des valeurs de la pression

## Representation graphique des resultats

def PisoT(z, T0 = 288):             # pour comparaison
    """ Calcule la valeur de P a l'altitude z selon le modele isotherme. 
    Par défaut, la température est fixée à 288 K. """
    return Psol*np.exp(-Mair*g*z/(R*T0))

plt.figure(figsize=(13,6))
plt.subplot(1,2,1)                  # graphique de gauche
plt.title("Evolution de la température avec l'altitude")
plt.plot(T, z*1e-3, 'b-')
plt.xlim(180,300), plt.xlabel(r"$T$ (en K)")
plt.ylim(0,85), plt.ylabel(r"$z$ (en km)")
plt.grid()
plt.subplot(1,2,2)                  # graphique de droite
plt.title("Evolution de la pression avec l'altitude")
plt.semilogx(P*1e-5, z*1e-3, 'b-', label = "Modèle ISA")
plt.semilogx(PisoT(z)*1e-5, z*1e-3,'m:', label = r"Modèle isotherme ($T=T_{\rm sol}$)")
plt.xlim(1e-6,1), plt.xlabel(r"$P$ (en bar)")
plt.ylim(0,85), plt.ylabel(r"$z$ (en km)")
plt.legend(loc=0)
plt.grid(which = 'both')
plt.show()

plt.figure()                        # comparaison des modèles ISA et isotherme
plt.plot(z*1e-3, abs(P-PisoT(z))/P, 'y-')
plt.xlim(0,20), plt.xlabel(r"$z$ (en km)")
plt.ylim(0,1), plt.ylabel(r"écart relatif : $\dfrac{\vert P-P_{\rm iso}\vert}{P}$")
plt.grid()
plt.show()