##TP8 noté partie 2
##euler
import numpy as np
import matplotlib.pyplot as plt # ça peut servir ...
#intiailisation des variables
Kt,Kn=1e-6,1e-6
m=10e-6
g=9.81
nu=0.5
gamma=5e-5

def F(X):
    return np.array([X[1],-Kt/m*X[0]-nu*g])
    
def G(Y):
    return np.array([Y[1],-Kn/m*Y[0]-gamma/m*Y[1]-g])
    
X1=np.array([0,0]) 
Y1= np.array([1e-5,0])
    
print(F(X1))
print(G(Y1))

#[  1.00000000e-03  -4.90500000e+00]
#[ 0.       -9.810001]
##
def euler(F,G, a, b, n):
    h = (b - a) / n
    T=np.linspace(a,b,n)
    X = np.zeros((n,2))
    Y = np.zeros((n,2))
    X[0]=X1
    Y[0]=Y1
    for k in range(n-1):
        X[k+1] = X[k]+h*F(X[k])
        Y[k+1] = Y[k]+h*G(Y[k]) 
    return T, X, Y
    
T,X,Y=euler(F,G, 0, 200, 200000)
print("temps",T)
print("x", X)
print("y",Y)


#mettre des titres aux axes et au graphique
plt.subplot(1,2,1)
plt.xlabel ('temps')
plt.ylabel ('deplacements')
#legende et correspondances des couleurs
plt.plot (T,X[:,0],'b',label=u'x(t)') 
plt.legend(loc=1)  

plt.subplot(1,2,2)
plt.xlabel ('temps')
#legende et correspondances des couleurs
plt.plot (T,Y[:,0],'r',label=u'y(t)')  
plt.legend(loc=1) 

##
plt.figure(2)
plt.subplot(1,2,1)
plt.xlabel ('x')
plt.ylabel ('y')
#legende et correspondances des couleurs
plt.plot (X[:,0],Y[:,0],'b',label=u'x en fonction de y ')
plt.legend(loc=1)

plt.subplot(1,2,2)
plt.xlabel ('y')
#portrait de phase
plt.plot (Y[:,0],Y[:,1],'r',label=u"y' en fonction de y") 
plt.legend(loc=1) 

## formule théorqiue à finir
# plt.figure(5)
# n=2000
# T,X,Y=euler(F,G, 0, 20, n)
# y_reel=[]
# ecart=0
# for i in range (len(T)):
#     b=np.sqrt((gamma/m/2)**2-(Kn/m)**2)
#     y_reel.append(np.cosh(np.sqrt(Kn/m)*T[i])*(np.exp(-gamma/m/2*T[i])))
#     #np.exp(-gamma/m/2*T[i])*(np.cos(b*T[i])+1/(b/np.sqrt((Kn/m)**2))*np.sin(b*T[i])))
#     #
#     
#     ecart+=abs(np.real(y_reel[i])-Y[i,0])
# plt.plot(T,y_reel)
# plt.plot(T,Y[:,0],'r')
# 
# def ecart(n):
#     T,X,Y=euler(F,G, 0, 0.2, n)
#     y_reel=[]
#     ecart=0
#     for i in range (len(T)):
#         y_reel.append(np.cos(np.sqrt(Kn/m)*T[i])*np.exp(-gamma/m/2*T[i])-1)
#         ecart+=abs(np.real(y_reel[i])-Y[i,0])
#     return ecart
#     
# n=np.linspace(2000,20000,10) 
# E=[]   
# for k in range(10):
#     E.append(ecart(int(n[k])))
#     
# plt.figure(3)
# plt.plot(n,E,'r+')
#