#N'hésitez pas à jouer avec les valeurs, changer les fonctions, les conditions initiales... C'est tout l'intérêt
#de l'outil informatique ! 
#Pour utiliser ce fichier : il est plus simple de récupérer seulement la question/l'exercice qui vous intéresse, et de
#le mettre dans votre interface Python, sinon tout va compiler simultanément... 


#Exercice 1 : 
#------------

import numpy as np
import matplotlib.pyplot as plt
#Exercice 1 : méthode d'Euler pour l'équation différentielle d'ordre 1
R=500
C=100e-6
tau=R*C
K=5
n=1000 #Nombre de points pour le tracé des courbes
t=np.linspace(0.00001,10*tau,n) #La liste des temps. On ne la fait pas commencer à exactement
#0 pour éviter les divergences.
E=K*(t/tau) #On peut changer cette fonction pour essayer d'autres second membres !
dt=10*tau/n
uc=np.zeros(n)
uc[0]=0
for i in range(n-1):
    uc[i+1]=uc[i]+dt*(-uc[i]/tau+E[i]/tau) #Relation de récurrence de la méthode d'Euler
    
plt.close()
plt.plot(t,uc) #Tracé de uc(t) après résolution
plt.show()


#Exercice 2 : 
#------------
#Résolution de l'équation du pendule 

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as sp

g=9.81
l=1
w0=np.sqrt(g/l)
n=1000
tmax=3*np.pi*2/w0
t=np.linspace(0.001,tmax,n)

def pendule(X,t): #X est le vecteur (theta,theta'), car Odeint ne peut résoudre que des eq d'ordre 1
    return (X[1],-w0*w0*np.sin(X[0]))

theta0deg=85
theta0=theta0deg*np.pi/180
derivtheta0=0
X_ini=(theta0,derivtheta0) #Vecteur conditions initiales

X_sol=sp.odeint(pendule,X_ini,t) #Odeint a besoin de la fonction de l'équadiff, des conditions 
#initiales, et de la liste des temps. Il itère tout seul, et renvoie le tableau X_sol.
#La première colonne de X_sol contient theta(t), et la deuxième contient dtheta/dt(t).

plt.close()
plt.plot(t,X_sol[:,0]) #Xsol[:,0] est la fonction theta(t), c'est la "première colonne" de X_sol
plt.show()

def annulation(liste,epsilon):  #Pour chercher le premier zéro, on parcourt theta(t) jusqu'à trouver
    for i in range(n):          #la première annulation, à epsilon près.
        if liste[i]<epsilon:
            return i*tmax/n
epsilon=0.001
T=annulation(X_sol[:,0],epsilon)*4 #La première annulation se produit au quart de la période ! (voir le graphe)
print(T) 

#--------------
# Tracé de T=f(theta_0)

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as sp

theta_deg=np.array([10,20,30,40,50,60,70,80,85,88])
theta=theta_deg*np.pi/180
T=np.array([2.022,2.022,2.046,2.070,2.118,2.167,2.215,2.287,2.335,2.359]) 
#Valeurs obtenues avec le programme précédent, mais on pourrait automatiser...
g=9.81
l=1
w0=np.sqrt(g/l)
T0=2*np.pi/w0
plt.close()
plt.plot(theta,T)
plt.show()
regression=np.polyfit(theta,T,2) #La fonction polyfit retourne une liste contenant les
#coefficients du polynôme : si T=ax^2+bx+c, polyfit renvoie [a,b,c].
print(regression[0]/T0) #On ne veut que le coefficient devant theta_0^2, qui est 
#le premier de la liste.

#Exercice 3 :
#------------

#Filtre d'ordre 1 appliqué à cos^4(w0*t)

import numpy as np
import matplotlib.pyplot as plt

n=10000 #Nombre de points pour les tracés
w0=10000
t=np.linspace(0,2*2*np.pi/w0,n)

wc=2*w0 #Pulsation de coupure du filtre : on veut garder 4*w0.
K=1
U0=6
ue=U0*(3/8+1/2*np.cos(2*w0*t)+1/8*np.cos(4*w0*t)) #Fonction ue(t) après linéarisation

def H(w,wc,K): #Définition de la fonction de transfert
    return (K*1j*w/wc)/(1+1j*w/wc)

us_complexe=U0*(3/8*H(0,wc,K)+1/2*np.exp(1j*2*w0*t)*H(2*w0,wc,K)+1/8*np.exp(1j*4*w0*t)*H(4*w0,wc,K))
us_reel=np.real(us_complexe) #On applique la fonction de transfert à chaque morceau de Ue et on 
# prend la partie réelle du résultat pour avoir Us

plt.close()
plt.plot(t,us_reel)
plt.plot(t,ue) #Ue sera tracée en orange sur le graphique.
plt.show()
#On constate que le filtre n'est pas assez sélectif : le signal de sortie n'est pas
#sinusoïdal... On peut essayer avec un filtre passe-bande, beaucoup plus sélectif

#-----------
#Bonus : Application d'un filtre passe-bande au même signal d'entrée 

import numpy as np
import matplotlib.pyplot as plt

n=10000
w0=10000
t=np.linspace(0,2*2*np.pi/w0,n)

K=1 #On peut changer la valeur de K pour accentuer 
U0=6
ue=U0*(3/8+1/2*np.cos(2*w0*t)+1/8*np.cos(4*w0*t))

def Hbande(w,wc,K,Q): #Filtre passe-bande d'ordre 2
    return K/(1+1j*Q*(w/wc-wc/w))

Q=100 #On rappelle qu'un grand facteur de qualité rend le filtre plus sélectif !
wc=4*w0 #Il faut cette fois-ci choisir comme fréquence du filtre celle qu'on veut garder
#Attention ici au terme constant : H diverge en w=0... on remplace 0 par 0.00001
us_passebande=U0*(3/8*Hbande(0.00001,wc,K,Q)+1/2*np.exp(1j*2*w0*t)*Hbande(2*w0,wc,K,Q)+1/8*np.exp(1j*4*w0*t)*Hbande(4*w0,wc,K,Q))
us_passebande_reel=np.real(us_passebande)

plt.close()
plt.plot(t,us_passebande_reel)
plt.plot(t,ue)
plt.show() 
#C'est beaucoup mieux : le signal de sortie est sinusoïdal

#----------
# Application au signal triangle

import numpy as np
import matplotlib.pyplot as plt

n=10000
w0=10000
t=np.linspace(0,2*2*np.pi/w0,n)

wc=400*w0
K=100 #On peut amplifier le signal si wc>>w0, si on veut le dériver (cf E6)
U0=6
def H(w,wc,K): #Définition de la fonction de transfert
    return (K*1j*w/wc)/(1+1j*w/wc)

ue_complexe=np.zeros(n)
us_complexe=np.zeros(n)
p=1000 #Nombre de termes dans la somme de Fourier
for i in range(p): #On itère sur p pour calculer chaque composante, et on en profite pour calculer la sortie
    ue_complexe=ue_complexe-U0*8/np.pi*(1/(2*i+1)**2)*np.exp(1j*(2*i+1)*w0*t)
    us_complexe=us_complexe - U0*8/np.pi*(1/(2*i+1)**2)*np.exp(1j*(2*i+1)*w0*t)*H((2*i+1)*w0,wc,K)

us_reel=np.real(us_complexe)
ue_reel=np.real(ue_complexe)
plt.close()
plt.plot(t,us_reel)
plt.plot(t,ue_reel)
plt.show() #On voit qu'on dérive le filtre si wc>>w0 : c'est ce qu'on avait vu dans E6.