# -*- coding: utf-8 -*-
"""
Created on Wed Sep  3 12:08:33 2025

@author: TSI2
"""

import numpy as np
import matplotlib.pyplot as plt

#Paramètre à modifier : donner une valeur négative à beta pour obtenir un système instable
#Coefficient de frottement 
beta = 1000 #kg.s-1 

#------------------------------------------------------------------------------------------
#Paramètres
k = 100 #kg.s-2 constante de raideur
m = 100 #kg masse
omega_0 = np.sqrt(k/m) #rad.s-1 pulsation de coupure amortisseur
f_0 = omega_0/(2*np.pi) #s-1 fréquence de coupure amortisseur
facteur_qualité = omega_0*beta/k

T = np.sqrt(k/m)  #s période propre système

T_e = T/3  #s période exictation
f_e = 1/T_e #Hz fréquence excitation
pas_temps = T_e/100 #s pas de temps
instants = np.arange(0,T_e*40,pas_temps)
e = 1*np.sin(2*np.pi/T_e*instants) # excitation 
#e = np.zeros(len(instants)) # excitation 

s = np.zeros(len(instants)) #initialisation du tableau des positions s
s[0] = 0 #m position initiale

s_prime = np.zeros(len(instants)) #initialisation du tableau des vitesse ds/dt
s_prime[0] = 0 #m vitesse initiale

#Méthode d'Euler
for k in range(len(s)-1):
    s[k+1] = s[k] + pas_temps*s_prime[k]
    s_prime[k+1] = s_prime[k] + pas_temps/m*( e[k] - k*s[k] - beta*s_prime[k]  )

#------------------------------------------------------------------------------------------
#Tracé de la figure
figure_1, axes = plt.subplots(2)
#figure_1.suptitle('f =  %.2E ; f_0 = %.2E'.format(f_e, f_0))
figure_1.suptitle("$f$ = {:.2e} Hz ; $f_0$ = {:.2e} Hz ; $Q$ = {:.2e}".format(f_e, f_0, facteur_qualité))
axes[0].plot(instants, e, label='forçage')
axes[0].set(ylabel='Forçage (N)')
axes[0].legend()
axes[0].grid()

axes[1].plot(instants, s, label='position')
axes[1].set(xlabel='Instant (s)', ylabel='Sortie (m)')
axes[1].legend()
axes[1].grid()

plt.show()