# TP 6 - MICHELSON EN LAME D'AIR - Etude du rayon des anneaux sombres
##################################################################################
# Cellule 1 : Importation des bibliotèques

import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rd
###############################################################################################################################
# Cellule 2 : Entrées des valeurs expérimentales

k = np.array([])         # Numéros des anneaux étudiés                                                   A COMPLETER

Dmin = np.array([])      # Diamètre minimum des anneaux en cm                                            A COMPLETER
Dmax = np.array([])      # Diamètre maximum des anneaux en cm                                            A COMPLETER      

R = 0.01*(Dmin + Dmax)/4          # Rayons moyens en m
deltaR = 0.01*(Dmax-Dmin)/4       # Demi-intervalle de variation de chaque valeur du rayon en m
uR = deltaR/np.sqrt(3)            # Incertitude-type sur les rayons en m  

y =                               # Fonction à tracer pour vérifier la loi R = f(k)                      A COMPLETER
uy =                              # Incertitude-type sur y                                               A COMPLETER
deltay = uy*np.sqrt(3)            # Demi-intervalle de variation de y
###############################################################################################################################
# Cellule 3 : Tracé des points de mesure avec les barres d'incertitude et tracé de la droite de régression linéaire

plt.errorbar(k,y,yerr=uy,fmt = ',',color='k',label='Points expérimentaux')   

p = np.polyfit(k,y,1)                                             # Régression linéaire y = ak + b
print('Pente de la droite a =',p[0],'m')                          # Valeur de a
print("Ordonnée à l'origine b = ",p[1],'m')                       # Valeur de b

plt.plot(k,p[0]*k+p[1],'r--',linewidth=1,label='Modèle affine')   # Tracé de la droite de régression

plt.title('Régression linéaire du rayon des anneaux au carré')
plt.xlabel('k')
plt.ylabel('R² en m²')
plt.legend()
plt.grid()
plt.show()

# L'ordonnée à l'origine est en théorie nulle. Il faut vérifier que b est très inférieur à la plus petite valeur de y
###############################################################################################################################
# Cellule 4 : Procédure de validation : analyse des écarts normalisés

résidus = y - (p[0]*k+p[1])   # Calcul des résidus (écart entre la valeur expérimentale et la droite de régression linéaire)
z = résidus/uy                # Calcul des écarts normalisés          

plt.plot(k,z,'o',color='k')    # Tracé des écarts normalisés
plt.title('Ecarts normalisés')
plt.xlabel('k')
plt.ylabel('Ecarts normalisés')
plt.grid()
plt.show()

# Si les écarts normalisés sont inférieurs à 2 en valeur absolue, la loi est validée
################################################################################################################################
# Cellule 5 : Simulation Monte-Carlo pour déterminer les incertitudes sur a et b

N = 50000             # Nombre d'expériences simulées
a = np.zeros(N)       # Initialisation du tableau dans lequel on va stocker les pentes a
b = np.zeros(N)       # Initialisation du tableau dans lequel on va stocker les ordonnées à l'origine b

for i in range(N):                          # Pour chacune des N expériences

  yi = rd.uniform(y-deltay,y+deltay)        # On simule la mesure des y pour l'expérience i
                                            # On prend une loi uniforme de valeur moyenne les y mesurés et de demi-intervalle deltay
                                            
  p = np.polyfit(k, yi, 1)                  # Régression linéaire de l'expérience i
  a[i] = p[0]                               # Pente pour l'expérience i
  b[i] = p[1]                               # Ordonnée à l'origine pour l'expérience i 
  
plt.subplot(1, 2, 1)                        # Histogramme des valeurs a
plt.hist(a, bins = 'rice', color = 'r')
plt.xlabel('Valeurs de a obtenues')
plt.ylabel("Nombre d'apparitions")
plt.grid()

plt.subplot(1, 2, 2)                        # Histogramme des valeurs de b
plt.hist(b, bins = 'rice', color = 'b')
plt.xlabel('Valeurs de b obtenues')
plt.ylabel("Nombre d'apparitions")
plt.grid()

plt.show()

amoy = np.mean(a)                           # a moyen
ua = np.std(a,ddof = 1)                     # Incertitude-type sur a 
print('Pente moyenne a =',amoy,'m')
print('Incertitude-type sur la pente u(a) =',ua,'m')

bmoy = np.mean(b)                           # b moyen
ub = np.std(b,ddof = 1)                     # Incertitude-type sur b
print("Ordonnée à l'origine moyenne b =",bmoy,'m')
print("Incertitude-type sur l'ordonnée à l'origine u(b) =",ub,'m')

# Les valeurs de a et de u(a) permettent alors de calculer l'épaisseur e de la lame d'air et son incertitude-type u(e)
#####################################################################################################################