import numpy as np
import matplotlib.pyplot as plt

#Vidange du reservoir

alpha=1e-3
g=9.81
Ldt=[0.1,1,5,10,30,60,250]

def euler(dt):
    h=[1.5]
    t=[0]
    while h[-1]>0:
        h.append(h[-1]-alpha*dt*np.sqrt(2*g*h[-1]))
        t.append(t[-1]+dt)
    return t,h

plt.figure()
for i in Ldt:
    T,H=[],[]
    T,H=euler(i)
    plt.plot(T,H,label="hauteur d'eau en fonction du temps, dt= "+str(i))
#solution analytique
T=np.arange(0,500,0.1)
plt.plot(T,[(-alpha*np.sqrt(g/2)*T[i]+np.sqrt(H[0]))**2 for i in range(len(T))])
plt.legend()

# exercice 2

from scipy.integrate import odeint
t=np.linspace(0,5,100)
def systE1(Y,t,g,l):
    theta,omega=Y
    dtheta=omega
    domega=-g/l*np.sin(theta)
    return [dtheta,domega]

def systE2(Y,t,g,l):
    theta,omega=Y
    dtheta=omega
    domega=-g/l*theta
    return [dtheta,domega]

CI=[[np.pi/4,0],[np.pi/20,0]]
g,l=9.81,0.3
for i in CI:
    plt.figure()
    solE1=odeint(systE1,i,t,args=(g,l))
    solE2=odeint(systE2,i,t,args=(g,l))
    plt.plot(t,solE1[:,0],label=('Solution E1, CI= '+str(i)))
    plt.plot(t,solE2[:,0],label=('Solution E2, CI= '+str(i)))
    plt.legend()
    plt.show()
    
#exercice 3

N=int(1e6)

#valeurs de l enonce:
ce=4180
m1=0.0542
m2=0.0573
T1=20.2
T2=47.3
Tf=31.2
um=0.0001
uT=0.1

#on calcule les valeurs de C avec les tirages aleatoires de valeurs

Lm1=np.random.normal(m1,um,N)
Lm2=np.random.normal(m2,um,N)
LT1=np.random.normal(T1,uT,N)
LT2=np.random.normal(T2,uT,N)
LTf=np.random.normal(Tf,uT,N)

LC=ce*(Lm2*(LT2-LTf)+Lm1*(LT1-LTf))/(LTf-LT1)

C=LC.mean()
uC=LC.std(ddof=1)

# on utilise l incertitude elargie pour l intervalle de confiance a 95%

print('la capacite thermique du calorimetre est de C = ('+str(C)+' ± '+str(2*uC)+') J/K')

#affichage de l'histogramme
plt.figure()
plt.hist(LC,bins='rice',color='blue',alpha=0.5)
plt.axvline(C,color='red')
plt.axvline(C-2*uC,linestyle='-',color='red')
plt.axvline(C+2*uC,linestyle='-',color='red')
plt.title('C = ('+str(C)+' ± '+str(2*uC)+')')
plt.show()

#exercice 4
N=int(1e6)
LL=np.random.normal(1,0.0025,N)
LT=np.random.normal(2.006,0.005,N)

Lg=4*np.pi**2*LL/LT**2

g=Lg.mean()
ug=Lg.std(ddof=1)

plt.figure()
plt.hist(Lg,bins='rice',color='blue',alpha=0.5)
plt.axvline(g,color='red')
plt.axvline(g-2*ug,linestyle='-',color='red')
plt.axvline(g+2*ug,linestyle='-',color='red')
plt.title('g = ('+str(g)+' ± '+str(2*ug)+')')
plt.show()

print('la valeur du champ de pesanteur est g = ('+str(g)+' ± '+str(2*ug)+') m/s/s')

g=9.81
Lm=np.random.normal(2.45e-3,0.005e-3,N)
LD=np.random.normal(8e-3,0.05e-3,N)
Lv=np.random.normal(0.625,0.03,N)
Lrho=np.random.normal(998,3,N)

LCd=8*Lm*g/np.pi/Lrho/LD**2/Lv**2

Cd=LCd.mean()
uCd=LCd.std(ddof=1)

plt.figure()
plt.hist(LCd,bins='rice',color='blue',alpha=0.5)
plt.axvline(Cd,color='red')
plt.axvline(Cd-2*uCd,linestyle='-',color='red')
plt.axvline(Cd+2*uCd,linestyle='-',color='red')
plt.title('Cd = ('+str(Cd)+' ± '+str(2*uCd)+')')
plt.show()

print('la valeur du coefficient de trainee est Cd = ('+str(Cd)+' ± '+str(2*uCd)+')')