# TP 11 - Mesure de viscosité - Programme 2
################################################################################
# Cellule 1 - Importation des bibliothèques

import numpy as np
import matplotlib.pyplot as plt
################################################################################
# Cellule 2 - Validation de la loi théorique donnant h(t) par régression linéaire

# Mesures
L =                                                                           # Longueur du tube en m             A COMPLETER
t = np.array([])                                                              # Tableau des temps t en s          A COMPLETER
hcm = np.array([])                                                            # Tableau des hauteurs h en cm      A COMPLETER
h = 0.01*hcm                                                                  # Tableau des hauteurs h en m

# Fonction à tracer pour vérifier la loi h(t)
Y =                                                                                 # A COMPLETER

# Précision des mesures : demi-intervalle de variation
deltaL =                            # L mesuré avec une règle graduée en mm         # A COMPLETER
deltat =                            # t mesuré grossièrement à la seconde près      # A COMPLETER
deltah =                            # h mesuré grossièrement à l'oeil au mm près    # A COMPLETER

# Incertitudes-type sur t en supposant une loi de probabilité uniforme
ut = deltat/np.sqrt(3)

# Incertitude-type sur Y calculée par simulation Monte-Carlo
N = 50000                                                                      # Nombre N d'expériences simulées
Ymcarlo = np.zeros([N,len(t)])                                                 # Initialisation du tableau des valeurs des Y[i,j] 
                                                                               #     N lignes : une pour chaque expérience simulée
                                                                               #     len(t)colonnes : une pour chaque instant tj
uY=np.zeros(len(t))                                                            # Initialisation du tableau des incertitudes-types sur Y à chaque instant tj

for j in range(0,len(t)):                                                      # Pour l'instant de mesure  tj

    for i in range (0,N):                                                      # Pour l'expérience numéro i
        hi = np.random.uniform(h[j]-deltah,h[j]+deltah)                        # Tirage aléatoire d'une valeur de h entre hj-deltah et hj+deltah
        Li = np.random.uniform(L-deltaL,L+deltaL)                              # Tirage aléatoire d'une valeur de L entre L-deltaL et L+deltaL
        Ymcarlo[i,j] = np.log(hi+Li)                                           # Calcul de Y à l'instant tj, pour l'expérience i
        
    uY[j] = np.std(Ymcarlo[:,j], ddof = 1)                                     # Ecart-type de tous les Yi à l'instant tj

# Tracé des points de mesure avec leurs barres d'incertitudes
plt.subplot(1,2,1)
plt.errorbar(t,Y,xerr=ut,yerr=uY,fmt = ',',color='k',label='Points expérimentaux') 

# Régression linéaire de Y en fonction de t
p = np.polyfit(t,Y,1)     
a = p[0]             # Pente
b = p[1]             # Ordonnée à l'origine

# Equation de la droite de régression Y = at + b
YRL = a*t+b

# Tracé de la droite de régression
plt.subplot(1,2,1)
plt.title('Régression linéaire de Y = Ln(h+L)')
plt.xlabel('t')
plt.ylabel('Y')
plt.plot(t,YRL,'r--',linewidth=1,label='Modèle affine')
plt.legend()  
plt.grid()
plt.show()

# Calcul des résidus et des écarts normalisés
résidu =                                         # Ecart entre la valeur de Y issue de la mesure et celle issue de la régression linéaire     # A COMPLETER
z =                                              # Ecart normalisé = résidu/incertitude-type sur Y                                            # A COMPLETER

# Tracé des écarts normalisés
plt.subplot(1,2,2)
plt.title('Ecarts normalisés')
plt.xlabel('t')
plt.ylabel('écarts normalisés')
plt.plot(t,z,'o',color='k')                      # Tracé des écarts normalisés
plt.grid()
plt.show()
################################################################################
# Cellule 3 - Estimation de l'incertitude-type sur a et b par simulation Monte-Carlo

# Précision sur Y : demi-intervalle de variation
deltaY = uY*np.sqrt(3)

N = 50000              # Nombre d'expériences simulées
a = np.zeros(N)        # Initialisation du tableau dans lequel on va stocker les N pentes a
b = np.zeros(N)        # Initialisation du tableau dans lequel on va stocker les N ordonnées à l'origine b

for i in range(0,N):                       # Pour chacune des N expériences
   
    ti = np.random.uniform(t-deltat,t+deltat)      # Tableau des temps t pour l'expérience i par tirage aléatoire entre t-deltat et t+deltat
    Yi = np.random.uniform(Y-deltaY,Y+deltaY)      # Tableau des Y pour l'expérience i par tirage aléatoire entre Y-deltaY et Y+deltaY

    pi = np.polyfit(ti,Yi,1)               # Régression linéaire pour l'expérience i
    
    a[i] = pi[0]                           # Pente pour l'expérience i
    b[i] = pi[1]                           # Ordonnée à l'origine pour l'expérience i

amoy = np.mean(a)                         # Pente moyenne          
ua = np.std(a, ddof = 1)                  # Ecart-type de a = Incertitude-type sur a
print('La pente moyenne vaut amoy =',amoy)
print("L'incertitude-type sur la pente vaut u(a) =",ua)

bmoy = np.mean(b)                          # Ordonnée à l'origine moyenne
ub = np.std(b, ddof = 1)                   # Ecart-type de b = Incertitude-type b
print("L'ordonnée à l'origine moyenne vaut bmoy =",bmoy)
print("L'incertitude-type sur l'ordonnée à l'origine vaut u(b) =",ub)
################################################################################
# Cellule 4 - Calcul de la viscosité du glycérol et son incertitude-type par simulation Monte-Carlo

# Valeurs numériques des paramètres de la formule donnant la viscosité
µ =         # Masse volumique du glycérol en kg/m3                              A COMPLETER   
g =         # Champ de pesanteur en m/s²                                        A COMPLETER
a = amoy    # Pente moyenne de la régression linéaire
L =         # Longueur du tube en m                                             A COMPLETER
d =         # Diamètre du tube en m                                             A COMPLETER
D =         # Diamètre du récipient en m                                        A COMPLETER

# Précision : demi-intervalle de variation
deltaa = ua*np.sqrt(3)              # Précision sur a = incertitude-type*racine de 3     
deltaL =                            # L mesuré avec une règle graduée en mm              A COMPLETER  
deltad =                            # d mesuré avec une règle graduée en mm              A COMPLETER
deltaD =                            # D mesuré avec une règle graduée en mm              A COMPLETER

N = 50000                           # Nombre d'expériences simulées
eta = np.zeros(N)                   # Initialisation du tableau dans lequel on va stocker les viscosités

for i in range(0,N):                # Pour chacune des N expériences
   
    ai =                            # Tirage aléatoire d'une valeur de a pour l'expérience i entre a-deltaa et a+deltaa       A COMPLETER
    Li =                            # Tirage aléatoire d'une valeur de L pour l'expérience i entre L-deltaL et L+deltaL       A COMPLETER
    di =                            # Tirage aléatoire d'une valeur de d pour l'expérience i entre d-deltad et d+deltad       A COMPLETER
    Di =                            # Tirage aléatoire d'une Valeur de D pour l'expérience i entre D-deltaD et D+deltaD       A COMPLETER

    eta[i] =                        # Calcul de la viscosité pour l'expérience i                                              A COMPLETER 

# Viscosité moyenne et incertitude-type
etamoy = np.mean(eta)                         # Viscosité moyenne          
ueta = np.std(eta, ddof = 1)                  # Incertitude-type sur la viscosité
print('La viscosité moyenne vaut : eta =',etamoy)
print("L'incertitude-type sur la viscosité vaut : u(eta) =",ueta) 

# Histogramme des N mesures de viscosité
plt.title('Histogramme de la viscosité')
plt.xlabel('viscosité')
plt.ylabel('Effectif')
plt.hist(eta,bins='rice')
plt.grid()
plt.show()
################################################################################