from math import *
import numpy as np
import matplotlib.pyplot as plt

## Snell-Descartes
# q1 On a théoriquement a=n1/n2, et b=0.

# q2 et q3
I1 = [0, 0.17, 0.35, 0.52, 0.70, 0.87, 1.05, 1.22]
I2 = [0, 0.12, 0.22, 0.35, 0.42, 0.51, 0.67, 0.71]

x = [sin(i1) for i1 in I1]
y = [sin(i2) for i2 in I2]
a,b = np.polyfit(x,y,1)
print(a,b)

plt.figure()
plt.plot(x,y,'o', label='données')
plt.plot(x, [a*x+b for x in x], label='modèle ax+b')
plt.legend()
plt.grid()
plt.show()

# Les points semblent bien alignés, donc on peut valider le modèle et la loi de Descartes.
# On a a=0.687, or a=n1/n2 et n1=1 (air), donc n2 = 1/a = 1.45.


## Vitesse du son
# q4 On pose x = sqrt(T) et y=c pour se ramener à une loi linéaire.
c = [310, 325, 331, 342, 358, 370, 378]
T = [240, 260, 280, 300, 320, 340, 360]

x = [T**.5 for T in T]
y = c
uy = [c*0.01 for c in c] # pour q5 : calcul de l'incertitude sur y
a,b = np.polyfit(x,y,1)
print(a,b)

plt.figure()
#plt.plot(x,y,'o') # pour q4 
plt.errorbar(x,y,yerr=uy,fmt='o', label='données') # pour q6 : tracé avec barres d'erreurs
plt.plot(x, [a*x+b for x in x], label='modèle ax+b')
plt.legend()
plt.grid()
plt.show()

# Les points semblent bien alignés, et avec la q6 et le tracé des barres d'incertitude
# on confirme que le modèle est validé (la droite passe par les barres ou presque).
# On obtient a = 19.75 et b = 3.4. La valeur de a est plus proche de celle du modèle adiabatique.


## Vitesse du son avec incertitudes sur a et b
# q7

liste_a = []
liste_b = []
N = 10000
for i in range(N):
    ytirage = []
    for k in range(len(y)):
        ytirage.append(np.random.uniform(y[k]-uy[k]*sqrt(3), y[k]+uy[k]*sqrt(3)))
    
    a,b = np.polyfit(x,ytirage,1)
    liste_a.append(a)
    liste_b.append(b)

a = np.mean(liste_a)
b = np.mean(liste_b)
ua = np.std(liste_a)
ub = np.std(liste_b)
print(a,ua)
print(b,ub)

plt.figure()
plt.errorbar(x,y,yerr=uy,fmt='o', label='données')
plt.plot(x, [a*x+b for x in x], label='modèle ax+b')
plt.legend()
plt.grid()
plt.show()

# q8
# on obtient a = 19.75 et u(a) = 3.4,
# et b = 3.4 et u(b) = 19.
# On remarque que la valeur 0 est bien dans l'intervalle b +/- u(b).
# Pour a, on calcule l'écart normalisé avec la valeur du modèle adiabatique :
print(abs(a-20.1)/ua)  # 0.28 < 2, modèle validé.
# et l'écart normalisé avec la valeur du modèle isotherme :
print(abs(a-16.9)/ua)  # 2.6 > 2, modèle rejeté.