# -*- coding: utf-8 -*-
"""
Created on Fri Sep 12 07:22:00 2025

@author: erwan
"""

import numpy as np
import matplotlib.pyplot as plt


#Paramètres à changer (positif, négatif, nul)
b_1 = -0.1


b_0 = 1
b_2 = 1

nombre_périodes = 30
s_lim = 15 #V valeur limites

omega_0 = np.sqrt(b_0/b_2) #rad.s-1 pulsation propre
f_0 = omega_0/(2*np.pi) #s fréquence propre

duree = nombre_périodes/f_0 #s durée de l'étude
dt = duree/1000 #s pas de temps

instants = np.arange(0, duree+dt, dt)

s = np.zeros(len(instants))
s[0] = 0.001 #V Initialisation petite valeur
s_prim = np.zeros(len(instants))



#Equation différentielle
#b_0 s + b_1 s_prim + b_2 s_second = 0
#s(k+1) = s(k) + dt*s_prim(k)
#s_prim(k+1) = s_prim(k) - dt*( b_0/b_2 s(k) - b_1/b_2 s_prim(k) )

for k in range(len(instants)-1):
    if b_1 == 0:
        s[k+1] = s[0]*np.cos(omega_0*(k+1)*dt )
    else:
        s[k+1] = s[k] + dt*s_prim[k]
        if np.abs(s[k+1]) >= s_lim:
            s[k+1] = s_lim*np.sign(s[k+1])
        s_prim[k+1] = s_prim[k] - dt/b_2*( b_0*s[k] + b_1*s_prim[k] )
    
#------------------------------------------------------------------------------------------
#Tracé de la figure
figure_1, axes = plt.subplots(1)
axes.plot(instants, s, label='sortie')
if b_1 < 0 :
    axes.plot(instants, s_lim*np.ones(len(instants)),'r--', label='saturation')
    axes.plot(instants, -s_lim*np.ones(len(instants)),'r--')
    axes.set_ylim(-s_lim*1.5,s_lim*1.5)
axes.set(ylabel='Signal $s$ (V)')
axes.set(xlabel='Instants $t$ (s)')
axes.set_title('$b_0 =$ 1, $b_1 = $'+str(b_1)+' , $b_2 =$ 1' )
axes.legend()
axes.grid()
plt.show()