import numpy as np
import matplotlib.pyplot as plt

""" Données expérimentales (avec incertitudes-types) """
lambd=511e-9    # Longueur d'onde en m

az0=4.05    # Azimut d'ordre 0 en degré
uaz0=0.03  # Incertitude a0 en degré

az_list=np.array([326.17,346.22,21.95,41.88])   # Liste des azimuts en degrés pour les différentes ordres
uazi=0.03   # Incertitude-type sur un azimut ai en degré

k_list=np.array([-2,-1,1,2])   # Liste des ordres d'interférences

""" Algorithme de Monte-Carlo pour un ajustement affine """

N=1000                    #Nombre de tirages simulés
M=len(k_list)   # Nombre de points mesurés (ici : 4 mesures)

a_list=np.zeros(N)    # Initialisation des listes à remplir pour les pas du réseau a de chaque régression
az_list_k=np.zeros(M)    # Initialisation des listes à remplir pour les azimuts de chaque tirage k

for k in range(N): # Procédure de tirage
    # Pour le tirage k, calcul aléatoire de chaque abscisse et chaque ordonnée
    # on choisit une distribution uniforme (rectangulaire) donc la fonction de 
    # tirage aléatoire à choisir est np.random.uniform
    # Rq : le np.sqrt(3) permet de passer de l'incertitude-type à la loi de probabilité
    # Si on avait choisi une distribution normale, 
    # il aurait fallu choisir np.random.normal
    
    az0_k=np.random.uniform(az0-np.sqrt(3)*uaz0,az0+np.sqrt(3)*uaz0)
    
    for i in range(M): 
        az_list_k[i]=np.random.uniform(az_list[i]-np.sqrt(3)*uazi,az_list[i]+np.sqrt(3)*uazi)
    
    y_list_k=np.sin((az_list_k-az0_k)*np.pi/180)
    
    cd,b=np.polyfit(k_list,y_list_k,1)  # Régression linéaire retournant cd et b tels que y=cd*x+b
    a_list[k]=lambd/cd
    
a_moy = np.mean(a_list)          #Calcul de la valeur moyenne
ua = np.std(a_list,ddof=1)       #Calcul de l'écart-type (ddof=1 pour correspondre à la définition de physique)

print('Pas du réseau a : ',a_moy, ' m')                  #Affichage des résultats
print('Incertitude-type u(a) : ',ua,' m')

plt.figure('Tirage du pas du réseau')
plt.hist(a_list,range=(min(a_list),max(a_list)),bins = 'rice')  # rice permet d'optimiser la taille des bins
plt.title('Histogramme des tirages aléatoires des pas du réseau')
plt.xlabel("a (m)")
plt.ylabel('Occurences')
plt.show()