# -*- 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

#Equations différentielles ordinaires (EDO) (bibliothèque directement intégrée)
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=20000 # Nombre de points pour discrétiser le temps

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

def dx(x,t):
    #fonction qui calcule la derivée et la derivée seconde (dans la
    #ligne i de la matrice x on stocke le déplacement et sa dérivée)
    x1= x[0]
    dx1= x[1]
    #stockage dans un plus petit vecteur pour plus de clarté (optionnel)
    X1=np.array([x1,dx1])
    #écriture de la matrice
    K=np.array([[0,1],[-k1/m1,0]])  
    #calcul du terme suivant (produit matriciel np.dot)
    dX1=np.dot(K,X1)
    #recupération de la derivée et de la dérivée seconde
    dx1 = dX1[0]
    ddx1 = dX1[1]

    return [dx1, ddx1]
    

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

# On créer une matrice avec des zeros
# le nombre de lignes correspond au nombre de points
# le nombre de colones (2) corespond au déplacement et à sa dérivée.
x=np.zeros((nbpt,2))
x[0,0]=x0 # On donne la position pour t=0

# On parcour tous les points de calcul
# pour chaque point le : permet de travailler a la fois sur le déplacement 
# et sa dérivée en prenant les deux termes correspondants de la matrice.
for i in range(len(t)-1):
    x[i+1,:]=x[i,:]+np.dot(dx(x[i,:],t[i]),t[i+1]-t[i])
    
#  <----- Affichage des résultats ----->
pl.plot(t, x[:, 0], 'b') # Affichage des positions en fonction du temps
pl.xlabel('t') # Titre des abscisses
pl.ylabel(u'Déplacement') # Titre des ordonnées
#pl.title(u'Déplacement de la tour en fonction du temps') # Titre
pl.show() # Affichage