# -*- coding: utf-8 -*-
"""
Created on Fri Jul  8 14:38:50 2022

@author: starons
"""

"""
Effectue une régression linéaire pour voir si la loi de Cauchy décrit 
correctement le comportement du prisme.La méthode de Monte-Carlo permet de 
simuler la variabilité des mesures de l'indice
pour évaluer les incertitudes-type sur les paramètres du modèle.
"""

# On importe les bibliothèques utiles
# -----------------------------------
import numpy as np
from matplotlib import pyplot as plt
from numpy import random

# On rentre les données expérimentales
# ------------------------------------
Lambda = np.array([404.7,435.8,480.0,546.1,578.1,615])      # longueur d'onde en nm 
n = np.array([1.7761,1.7652,1.7475,1.7358,1.7294,1.7229])   # indice de réfraction
X = 1/Lambda**2                                             # création d'une variable adéquate

# On procède à une première régression linéaire
# ---------------------------------------------
A,B = np.polyfit(X,n,1)      

# On introduit les incertitudes-type sur les mesures de l'indice
# --------------------------------------------------------------
u_n = 0.0034 # la même pour tous les valeurs... douteux !

# On procède à N regressions linéaires
# ------------------------------------
N = 1000
a = np.zeros(N)                       # pour stocker les N pentes
b = np.zeros(N)                       # pour stocker les N ordonnées à l'origine 
for i in range(N):
    nmc = random.normal(n,u_n)        # On utilise ici une loi normale pour simuler la i ème série de 6 valeurs de n 
    a[i],b[i] = np.polyfit(X,nmc,1)   # i ème régression linéaire 
pente = np.mean(a)
oo = np.mean(b)

# Représentations graphiques
# --------------------------
plt.errorbar(X,n,yerr=u_n,fmt='.',label="points expérimentaux")
plt.plot(X,A*X+B,label="droite de régression simple")
plt.plot(X,pente*X+oo,label="droite de régression après simulation Monte-Carlo")
plt.xlabel("X"),plt.ylabel("n")
plt.title("Validation de la loi de Cauchy")
plt.grid(),plt.legend()

# Affichage des résultats
# -----------------------
print ("Après simulation de Monte Carlo pour tenir compte de la variabilité des mesures, on a :")
print("A = {:.0f} +/- {:.0f}".format(np.mean(a),np.std(a,ddof=1)))
print("B = {:.4f} +/- {:.4f}".format(np.mean(b),np.std(b,ddof=1)))
print("Avec la régression simple on avait :")
print("A = ",A) 
print("B = ",B)