# -*- coding: utf-8 -*-
"""
Created on Tue Jun  4 16:07:37 2019

@author: tmsth
"""

import numpy as np
import matplotlib.pyplot as plt

import os
donnees=open('traj.txt',"r")
lignes=donnees.readlines()
donnees.close()
taille=len(lignes)


z1=[]
z2=[]
pas=1./30.
for i in range (len(lignes)-1):
    l=lignes[i]
    l=l.strip()


for i in range (len(lignes)-1):
    l=lignes[i]
    l=l.strip()
    l=l.split('\t')
    for i in range (len(l)):
        l[i]=l[i].replace(',','.')
    z1.append(float(l[0]))
    z2.append(float(l[1]))


#L'appareil prend 30 photos par seconde, je crée une liste des temps
#de la même taille que la liste des mesures de position
temps=[i*pas for i in range(len(z1))]


#je trace la courbe position en fonction du temps
plt.plot(temps,z1,'*')
plt.grid()
plt.xlabel('t(s)')
plt.ylabel('y(m)')
plt.show()    
    
#Le régime permanent est rapidement atteint, je crée une liste position qui
#ne prend que les positions à partir du 10 ième point de mesure: soit les 
#positions une fois le régime permanent atteint

position=[z1[i] for i in range(10,len(z1))]


#je calcule la vitesse en chaque point de la liste position en utilisant
# v=Delta z/ Delta t =(position[i+1]-position[i])/pas


vl=[]
for i in range(len(position)-1):
    vl.append((position[i+1]-position[i])/pas)
    
#Je calcule la vitesse limite moyenne et l'inccertitude sur cette vitesse moyenne    
vlmoy=np.mean(vl)
uvl=np.std(vl)

#j'en déduis h: la vitesse limite vl vérifie
#g*(1-(meau/m)))-(h*vl**2/m)=0 on en déduit h=m*g*(1-(meau/m))/vl**2

m=0.200
um=1e-4
R=2.67e-2
uR=0.05e-2
rho=1e3
meau=rho*4*np.pi*R**3/3
g=9.81
hmoy=m*g*(1-(meau/m))/vlmoy**2

print('h=',hmoy)

#J'applique Monte Carlo pour calculer l'incertitude sur h:



# j'applique la méthode d'Euler avec la valeur trouvée pour h pour comparer
# la trajectoire théorique et la trajectoire expérimentale


z0=0.
v0=0.
vth=[v0]
zth=[z0]
for i in range(0,len(z1)-1):
    a=g*(1-(meau/m))-(hmoy*vth[i]**2/m)
    vth.append(vth[i]+a*pas)
    zth.append(zth[i]+vth[i]*pas)

plt.plot(temps,z1,'*')
plt.plot(temps,zth)
plt.grid()
plt.xlabel('t(s)')
plt.ylabel('y(m)')
plt.show()    
        

#Même chose pour z2
plt.plot(temps,z2,'*')
plt.grid()
plt.xlabel('t(s)')
plt.ylabel('y(m)')
plt.show()    
    
#Le régime permanent est rapidement atteint, je crée une liste position qui
#ne prend que les positions à partir du 10 ième point de mesure: soit les 
#positions une fois le régime permanent atteint

position=[z2[i] for i in range(10,len(z2))]


#je calcule la vitesse en chaque point de la liste position en utilisant
# v=Delta z/ Delta t =(position[i+1]-position[i])/pas


vl=[]
for i in range(len(position)-1):
    vl.append((position[i+1]-position[i])/pas)
    
#Je calcule la vitesse limite moyenne et l'inccertitude sur cette vitesse moyenne    
vlmoy=np.mean(vl)
uvl=np.std(vl)

#j'en déduis h: la vitesse limite vl vérifie
#g*(1-(meau/m)))-(h*vl**2/m)=0 on en déduit h=m*g*(1-(meau/m))/vl**2

m=0.200
um=1e-4
R=2.67e-2
uR=0.05e-2
rho=1e3
meau=rho*4*np.pi*R**3/3
g=9.81
hmoy=m*g*(1-(meau/m))/vlmoy**2

print('h=',hmoy)