import numpy as np
import matplotlib.pyplot as plt
from math import pi, sqrt, cos, atan

"""La fonction Filtrage(f,k) permet de calculer la réponse de l'harmonique de rang k pour le signal triangulaire de fréquence f.
Il faut donner la valeur de la fréquence f du signal triangulaire vers la ligne 54 ainsi que le rang des plus petit et plus grand harmoniques que l'on veut prendre en compte ligne 56"""

R = 500
C = 30e-9
omega0 = 1 / R / C
f0 = omega0 / 2 / pi
E = 10

def Gain(f):
    return 3 / sqrt(1 + (f / f0 - f0 / f)**2)

def Phase(f):
    return -atan(f / f0 - f0 / f)

def entree_a(f,k):
    if k%2 == 0:
        return 0
    else:
        return 8 * E / k**2 / pi**2

def sortie_a(f,k):
    if k%2 == 0:
        return 0
    else:
        return 8 * E / k**2 / pi**2 * Gain(k * f)

def sortie_p(f,k):
    if k%2 == 0:
        return 0
    else:
        return Phase(k * f)

def Filtrage(f,k):
    print("a_",k, "=", entree_a(f,k), " G(",k,"*",f,")*a_",k, "=", sortie_a(f,k), " phi(",k,"*",f,")=", sortie_p(f,k), "rad =", sortie_p(f,k) * 180 / pi,"°",sep='')

def somme_sortie(f,kmin,kmax,t):
    res = 0
    for i in range(kmin,kmax + 1):
        res = res + sortie_a(f,i) * cos(2 * pi * i * f * t + sortie_p(f,i))
    return res

def somme_entree(f,k,t):
    res = 0
    for i in range(k + 1):
        res = res + entree_a(f,i) * cos(2 * pi * i * f * t)
    return res

# Choix de la fréquence
f = 10e3 #2-a 
#f = 4e3 #2-b
#f = 20e3 #2-c

T = 1 / f # période
kmin = 1 # Choix de l'harmonique de plus bas rang
kmax = 3  # Choix de l'harmonique de plus haut rang

t = np.linspace(0, 4 * T, 400)
s = [] # s est la sortie complète
sp = [] # sp est la sortie partielle
e = [] # e est l'entrée
for tt in t:
    sp.append(somme_sortie(f,kmin, kmax, tt))
    s.append(somme_sortie(f,0, 300, tt))
    e.append(somme_entree(f, 30, tt)) # Pour e on met plein d'harmoniques...

plt.figure("Tracé réponse, entrée")
plt.plot(t,s)
plt.plot(t,sp)
plt.plot(t,e)
plt.show()
##
# Après calcul on peut aussi faire exécuter:
for k in range(10):
    Filtrage(f, k)
# pour visualiser les résultats avec les valeurs numériques

##
# Reste à afficher le spectre et le diagramme de bode
xr=np.logspace(2, 6, 400) # xr pour pulsation réduite
H = 3 / (1 + 1j * (xr / f0 - f0 / xr))
G = np.abs(H)
GdB = 20 * np.log10(G)
Phi = np.angle(H) * 180 / pi

# On stocke dans spectre_e le spectre de raies de l'entrée
# et dans spectre_s celui de la sortie
spectre_e = []
spectre_s = []
f_spectre_e = [] # les fréquences strockées dans le spectre
f_spectre_s = [] # les fréquences strockées dans le spectre

for k in range(1,10):
    f_spectre_e.append(k * f)
    f_spectre_s.append(k * f *1.015) # On multiplie par 1.015 pour décaler un peu les raies spectrales des deux signaux
    spectre_e.append(-20 + 2 * entree_a(f,k))
    spectre_s.append(-20 + 2 * sortie_a(f,k))

fig = plt.figure('Bode', figsize=(6,6))
ax1 = fig.add_subplot(1,1,1, xscale='log',
                      title='Gain (dB)',
                      xlabel=u'fréquence (Hz)')
ax1.plot(xr, GdB)
plt.vlines(f_spectre_e,-20,spectre_e,'g')
plt.vlines(f_spectre_s,-20,spectre_s,'r')

plt.grid()
plt.show()
#ax2 = fig.add_subplot(1,2,2, sharex=ax1,  # sharex -> synchronise les échelles. Pratique pour zoomer !
#                      title=u'Phase (°)',
#                      xlabel=u'fréquence (Hz)')
#ax2.plot(xr, Phi, 'g')
## ajustement des "ticks"
##ax2.set_yticks([-90, -45, 0])
#ax2.set_yticks(np.linspace(-90,90,19))
## ou de façon plus automatique :
## ax2.set_yticks(range(0,-181,-45)) # (de 0 à -180°, tous les 45°)
#plt.grid()
#fig.tight_layout() # "étire" les graphiques dans la fenêtre
##plt.savefig("\home\papa\Documents\14-15\Exercices\Electronique\Test.eps")