## Importation des modules principaux
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as si 

##Fonction de Wien
def f(X,t):
    #on extrait les valeurs de s et sa dérivée
    s=X[0]
    derivee_de_s=X[1]
    #on définit la dérivée de s et sa dérivée 
    dX=np.zeros(2)  # on crée un tableau numpy de taille 2 que l'on remplit ensuite...
    
    dX[0]=X[1]
    dX[1]=-w0**2*s-1/Q*w0*derivee_de_s+w0*derivee_de_e(s,derivee_de_s)
    return dX


## Modélisation de l'ALI

def derivee_de_e(s,ds):
    if np.abs(A*s)>Vsat:
        return 0
    else:
        return ds*A
    
def sortie_ALI(s):
    if np.abs(A*s)>Vsat:
        return np.sign(s)*Vsat
    else:
        return s*A
sortie_ALI=np.vectorize(sortie_ALI)

## Définition des paramètres importants du système

T=1 #période de l'oscillateur
w0=2*np.pi/T #pulsation propre associée
Q=1/3 #facteur de qualité
A=3.2 #gain de l'amplificateur non inverseur
Vsat=15 #tension de saturation
Npoints=1000 #nombre de pas de temps 
tmax=10*T #temps max

## Oscillateur à pont de Wien 

cini=[.1,0] #Conditions initales
t=np.linspace(0,tmax,Npoints) #définition du vecteur temps
solution=si.odeint(f,cini,t) #utilisation de odeint pour résoudre 
s_sol=solution[:,0] #on extrait s 
ds_sol=solution[:,1] #et sa dérivée

e=sortie_ALI(s_sol) #on peut aussi calculer la sortie de l'ALI


## On passe aux figures
plt.close('all')

plt.plot(t,s_sol)
plt.title("Graphe de s(t)")
plt.xlabel("t")
plt.ylabel("s(t)")
#plt.xlim(3,6)#permet de zoomer sur une partie du temps si on le souhaite
plt.ylim(-5.5,5.5)# permet de choisir l'échelle verticale
# 
# plt.figure()
# plt.plot(t,e)
# plt.title("Graphe de e(t)")
# plt.xlabel("t")
# plt.ylabel("e(t)")
# #plt.xlim(3,6)#permet de zoomer sur une partie du temps si on le souhaite
# plt.ylim(-15.5,15.5)# permet de choisir l'échelle verticale

plt.figure()
plt.plot(s_sol,ds_sol)
plt.title("Portrait de phase s'(s)")
plt.xlabel("s(t)")
plt.ylabel("s'(t)")

plt.show()