## Calcul d'intégrale et lecture de données
import numpy as np
import matplotlib.pyplot as plt

##########################
## Importation des données
##########################
## Si le dossier courant de travail (cwd pour current working directory) est bien défini
# (ex: si on ouvre l'éditeur Python par un clic droit - Ouvrir avec sur le fichier .py du programme) 
file = 'Radiosondage_Trappes_20211102_12h_tab.csv'
# données obtenues sur https://donneespubliques.meteofrance.fr

radiosondage = np.loadtxt(file,delimiter='\t',skiprows=3)
# Remarques:
# _ici on n'importe pas les 3 premières lignes (étiquettes)
# _dans ce fichier .csv, le séparateur de colonnes est ici la tabulation ( \t ),
# mais peut être aussi la virgule, le point virgule ou l'espace: ',' ou ';' ou ' '

nb_niv = np.shape(radiosondage)[0] # nombre de niveaux verticaux

## Si le chemin d'accès au fichier est à définir:
#import os # module d'interaction avec le système de fichier de l'OS utilisé
#os.getcwd() # pour connaître le dossier courant de travail
#os.chdir("C:\\Users\\Home\\TP_info") # pour définir le cwd sous windows (les \ doivent être doublés)
#os.chdir("/home/armel/TP_info") # pour définir le cwd sous UNIX (Linux, MacOS)

## Remarque: éviter les caractères spéciaux (accents, cédilles, espaces, apostrophes...)
## dans les noms de fichiers s'ils sont manipulés par des programmes, cela simplifie...

## On sépare les différentes variables physiques
P = radiosondage[:,0] # pression mesurée (Pa)
z = radiosondage[:,1] # altitude mesurée (m)
T = radiosondage[:,2] # température mesurée (K)

plt.plot(T,z/1000,'r-')
plt.title('Température')
plt.xlabel('Température (K)')
plt.ylabel('Altitude (km)')
plt.grid()
plt.show()

## CALCUL DE LA PRESSION HYDROSTATIQUE RÉELLE
M = 29e-3 # kg.mol^-1 masse molaire de l'air
g = 9.81 # m.s^-2 champ de pesanteur à Paris
R = 8.314 # J.K^-1.mol^-1 constante des gaz parfaits
Phyd = np.ones(nb_niv) # création du tableau des pressions (rempli de 1)
Phyd[0] = P[0] # pression de surface
# Intégration par la méthode des rectangles à gauche, à pas variable...
S = 0 # initialisation de la somme
for k in range(nb_niv-1):
    S = S + (z[k+1] - z[k]) / T[k]
    Phyd[k+1] = Phyd[0] * np.exp( -M*g/R * S )

## CALCUL DE LA PRESSION HYDROSTATIQUE ISOTHERME
Tmoy = T[0] # Température supposée uniforme = température de surface
Piso = P[0] * np.exp( -M*g/R/Tmoy * (z-z[0]) )

plt.plot(Phyd,z/1000,'r-',label='Calculée')
plt.plot(P,z/1000,'k-',label='Donnée (radiosondage)')
plt.plot(Piso,z/1000,'g-',label='Isotherme')
plt.title('Pression hydrostatique')
plt.xlabel('Pression (Pa)')
plt.ylabel('Altitude (km)')
plt.legend(loc='upper right')
plt.grid()
plt.show()


## À faire en plus:
## tracer la température de rosée, la hampe de vent en fonction de l'altitude (avec longueur du vecteur ?)