import numpy as np
import matplotlib.pyplot as plt
plt.close('all')


### ÉTUDE DYNAMIQUE
### Mesure de la durée de 20 oscillations pour différentes masses suspendues
m = np.array([0, 5, 10, 20, 30, 40, 50, 60, 80, 100]) # g

T = np.array([7.95, 8.81, 10.27, 11.78, 13.00, 14.33, 14.97, 15.78, 17.16, 18.36])/20 # s
u_T = 1/20

# plt.figure()
# plt.errorbar(m,T,yerr=u_T,marker='o',ls='None')
# plt.xlabel('Masse suspendue (g)')
# plt.ylabel("Période d'oscillation (s)")

T2 = T**2
u_T2 = 2*T*u_T

a,b = np.polyfit(m[0:5],T2[0:5],1)

plt.figure()
plt.errorbar(m,T2,yerr=u_T2,marker='o',ls='None')
plt.plot(m, a*m+b, 'r-')
plt.xlabel('Masse suspendue (g)')
plt.ylabel("T^2 (s^2)")
plt.show()

k_exp = 4*np.pi**2 / a
m0_exp = b/a

print('k =', k_exp*1e-3, 'N.m-1')
print('m0 =', m0_exp, 'g')


### Estimation des incertitudes sur k et m0 par simulation Monte-Carlo :
N_simu = 10000
k_sim = []
m0_sim = []

for n in range(N_simu):
    T2_sim = np.random.normal(T2,u_T2)
    a,b = np.polyfit(m[0:5],T2_sim[0:5],1)
    k_sim.append(4*np.pi**2 / a)
    m0_sim.append(b/a)
    
u_k_exp = np.std(k_sim)
u_m0_exp = np.std(m0_sim)

print('u_k =', u_k_exp*1e-3, 'N.m-1')
print('u_m0 =', u_m0_exp, 'g')
    

# Valeurs attendues
m_vib = 79.95 * 48/53 # réglet de 53 cm fixé à la graduation 5 cm
m_eq = 33/140 * m_vib
print('m_dyn_th =', m_eq, 'g')


### ÉTUDE STATIQUE
### Mesure de lé déflexion en fonction de la masse suspendue

m = np.array([0, 5, 10, 20, 30, 40, 50]) # g
delta = np.array([8.7, 10.0, 11.3, 13.1, 15.0, 17.1, 18.6]) # cm

u_delta = 0.3/np.sqrt(3) # mesure précise à 3 mm près

g = 9.81    # m.s-2

a,b = np.polyfit(m,delta,1)
k_stat = g/a
m0_stat = b/a

plt.figure()
plt.errorbar(m,delta,yerr=u_delta,marker='o',ls='None')
plt.plot(m, a*m+b, 'r-')
plt.xlabel('Masse suspendue (g)')
plt.ylabel("Déflexion (cm)")
plt.show()

print('k_stat =', k_stat/10, 'N.m-1')
print('m0_stat =', m0_stat, 'g')

# Valeur attendue
print('m_stat_th =', 3/8 * m_vib, 'g')