# -*- coding: utf-8 -*-

# <----- Importation des bibliothèques ----->

import numpy as np #Bibliothèque utilisée pour le calcul
import pylab as pl #Bibliothèque utilisée pour les affichages

# Bibliothèque pour résoudre les
# équations différentielles ordinaires (EDO) 
from scipy.integrate import odeint 

# <----- Mise en données ----->

m1=264000000. # Masse équivalente du batiment
k1=225000000. # Raideur du batiment
x0=0.25 # Sollicitation initiale
temps1=240 # Temps maximum d'étude
nbpt=200000 # Nombre de points pour discrétiser le temps
m2= <-COMPLETER-> # (float)
k2= <-COMPLETER->
f= <-COMPLETER->

# <----- Ecriture de la dérivée ----->

def dx(x,t):
    # Fonction qui calcule la dérivée et la dérivée seconde
    # dans la ligne i de la matrice x on stocke le déplacement
    # de la tour, celui de l'amortisseur et leurs 2 dérivées. 
    x1= x[0]        # déplacement de la tour
    x2= x[1]        # déplacement de l'amortisseur
    dx1= x[2]       # vitesse de la tour
    dx2= x[3]       # vitesse de l'amortisseur
    
    #On les stocke dans un plus petit vecteur pour plus de clarté (optionnel)
    X1=<-COMPLETER->
    #Ecriture de la matrice
    K= <-COMPLETER->
    #Calcul du terme suivant (produit matriciel np.dot)
    dX1=<-COMPLETER->
    #Récupération des dérivées et des dérivées secondes
    <-COMPLETER->
    
    
    
    return [<-COMPLETER->]
    

# <----- Résolution ----->
 
# Découpage de l'intervalle avec un pas de temps régulier
t = np.linspace(0, temps1, nbpt)

# Résolution avec euler explicite
x= <-COMPLETER->
x[0,0]=x0
for i in range(len(t)-1):
    x[i+1,:]= <-COMPLETER->
    
# Résolution avec Scipy.integrate
x0_scipy = [x0, 0, 0, 0]
x_scipy = odeint(dx, x0_scipy, t)


#  <----- Affichage des résultats ----->

pl.plot(t, x[:, 0], 'b'label='Tour') # Affichage des positions de la tour en fonction du temps
pl.plot(t, x_scipy[:, 0], 'r') # Affichage de la solution avec Scipy.integrate
pl.xlabel('t') # Titre des abscisses
pl.ylabel('Déplacement') # Titre des ordonnées
pl.title('Déplacement de la tour en fonction du temps') # Titre
pl.show() # Affichage

<-COMPLETER-> # Affichage des positions de l'amortisseur en fonction du temps
<-COMPLETER-> # Affichage de la solution avec Scipy.integrate

pl.xlabel('t') # Titre des abscisses
pl.ylabel('Déplacement') # Titre des ordonnées
pl.title('Déplacement de l amortisseur en fonction du temps') # Titre
pl.legend()
pl.show() # Affichage