# -*- coding: utf-8 -*-
# Created on Thu Sep 12 08:15:00 2024
# @author: mhebding
# TP4 - Methodes numeriques



## 1 - Preparation

# 0.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

m = 5.3*10**-3 # masse du volant
g = 9.81 # acceleration de la pesanteur terrestre
Te = 45*10**-3 # periode d'echantillonnage

# 1.
pointage = np.loadtxt("pointage.txt",skiprows=1) # lis le fichier en sautant la 1ere ligne des titres
Xexp = pointage[:,5] # affecte toutes les lignes de la colonne 5 a Xexp
Yexp = pointage[:,6] # affecte toutes les lignes de la colonne 6 a Yexp
plt.figure('Donnees experimentales')
plt.plot(Xexp,Yexp, 'o', label='Donnees brutes')
plt.legend()
x = Xexp - Xexp [0] # replace l'origine au point de depart
y = -(Yexp - Yexp[0]) # idem et mets l'axe dans le bon sens

# 2.
def exp():
    plt.figure('Données expérimentales')
    plt.plot(x,y, 'o', label='Trajectoire du volant')
    plt.legend()
    plt.show()

exp()



## 2 - Modelisation
"""
Voici le minimum theorique, voir cours de physique pour le reste des questions
dv/dt + a/m v = g (1)
dv/dt + b/m v**2 = g (2)

En regime permanent dv/dt = 0, il vient donc :
(1) a = m*g/Vlim
(2) b = m*g/Vlim**2
"""

# 3.
Vxlim = (x[-1]-x[-2])/Te # v = xn+1-xn/dt
Vylim = (y[-1]-y[-2])/Te # idem
Vlim = np.sqrt(Vxlim**2+Vylim**2) # valeur de l'auteur 6.7 m/s
print('Vlim= '+str(round(Vlim,3))) # 6.799 VS 6.7

a = m*g/Vlim
b = m*g/Vlim**2

# 4.
Vx0 = (x[1]-x[0])/Te # meme strategie mais a t = 0
Vy0 = (y[1]-y[0])/Te # idem
V0 = np.sqrt(Vx0**2+Vy0**2) # valeur de l'auteur 40 m/s
print('V0= '+str(round(V0,3))) # 39.096 VS 40



## 3 - Methode d'Euler

# 5.
def A_lin(Vx,Vy):
    Ax = - a*Vx/m
    Ay = - g - a*Vy/m
    return Ax, Ay

def A_quad(Vx,Vy):
    V = np.sqrt(Vx**2 + Vy**2)
    Ax = - b*V*Vx/m
    Ay = - g - b*V*Vy/m
    return Ax, Ay

# 6.
dt = Te/50 # pas
N  = int((Te/dt)*len(x)) # nombre de points

def euler_V(acceleration):
    t = 0 # initialisations
    Vx = Vx0
    Vy = Vy0
    liste_t = [t]
    liste_Vx = [Vx]
    liste_Vy = [Vy]
    for n in range(N):
        t += dt
        liste_t.append(t)
        Ax, Ay = acceleration(Vx, Vy)
        Vx += Ax*dt
        Vy += Ay*dt
        liste_Vx.append(Vx)
        liste_Vy.append(Vy)
    return liste_t, liste_Vx, liste_Vy

# 7.
def euler_OM(Vx,Vy):
    x = 0
    y = 0
    liste_X = [x]
    liste_Y = [y]
    for n in range(len(Vx)):
        x += Vx[n]*dt
        liste_X.append(x)
        y += Vy[n]*dt
        liste_Y.append(y)
    return liste_X, liste_Y

# 8.
t, Vx_lin, Vy_lin = euler_V(A_lin) # resolution de v pour le modele lineaire
X_lin, Y_lin = euler_OM(Vx_lin, Vy_lin) # positions pour le modele lineaire

t, Vx_quad, Vy_quad = euler_V(A_quad) # resolution de v pour le modele quadratique
X_quad, Y_quad = euler_OM(Vx_quad, Vy_quad) # positions pour le modele quadratique

def comparaison():
    plt.figure('Comparaison des modèles aux données')
    plt.plot(x,y, 'o', label='Trajectoire du volant') # données experimentales
    plt.plot(X_quad, Y_quad, 'r', label='Modèle quadratique') # modele v**2
    plt.plot(X_lin, Y_lin, 'g', label='Modèle linéaire') # modele v
    plt.legend()
    plt.show() # affichage

comparaison()

# 9.
def auteur():
    V0 = 40
    Vx0 = 40*np.cos(52)
    Vy0 = 40*np.sin(52)
    comparaison()

auteur()



## 4 - Resolution avec une bibliotheque

# 10.
"""
Cinem est un tableau numpy (= une liste de listes) qui contient les vecteurs cinematiques du volant a l'instant t
Cinem[0] = liste des x(t)
Cinem[1] = liste des y(t)
Cinem[2] = liste des Vx(t)
Cinem[3] = liste des Vy(t)

DerivC est une fonction qui, connaissant Cinem a l'instant t renvoie sa dérivee (idem methode d'Euler)
"""

def DerivC(Cinem,t): # prends X, Y, Vx, Vy
    Ax, Ay = A_quad(Cinem[2],Cinem[3])
    return [Cinem[2],Cinem[3], Ax, Ay] # retourne Vx, Vy, Ax, Ay

t = np.linspace(0,Te*len(x),N)
Cinem = odeint(DerivC,[0,0,Vx0,Vy0],t)
X_odeint = Cinem[:,0]
Y_odeint = Cinem[:,1]

def exp_euler_odeint():
    # INSTRUCTIONS GRAPHES EXP, EULER VS ODEINT
    plt.figure('Comparaison Euler et odeint')
    plt.plot(x, y, marker='o', label='expérimental')
    plt.plot(X_quad, Y_quad, '-', color='b', label='euler')
    plt.plot(X_odeint, Y_odeint, color='red', label='odeint')
    plt.legend()
    plt.show()

exp_euler_odeint()