from math import*
import matplotlib.pyplot as plt
import numpy as np
T0,U0=1E-03,0.5
w0=2*pi/T0
def F(t):# Fonction rectangle impaire
    if t<T0/2:
        return U0
    else:
        return(-U0)
def G(t):# Fonction Dents de scie paire
    if t<T0/2:
        return U0-4*U0*t/T0
    else:
        return -3*U0+4*U0*t/T0
def f(t):# Choix de F(t) ou G(t)
    return F(t)

def integrale(f,n,a,b):
    S,t,N=0,a,500
    dx=(b-a)/N
    while t<b:
        S+=f(t)*dx*np.exp(-1j*w0*n*t)
        t+=dx
    return S

#Coefficients
a0=integrale(f,0,0,T0)*1/T0
CN=[a0]
Nmax=20# Nombre d'harmoniques dans la série
for n in range(Nmax):
    cn=integrale(f,n+1,0,T0)*2/T0
    CN.append(cn)

def coeff(N):# Liste des N premiers coefficents an et bn
    for n in range(N+1):
        print("c{} = {:.2f}\n".format(n,CN[n]))

def e_(t):# Addition des Nmax premiers harmoniques = s(t)
    e=0
    for n in range(Nmax+1):
        e=e+CN[n]*np.exp(1j*w0*n*t)
    return e
def e(t):
    return np.real(e_(t))

def H_(w):
#    Q,Wc=0.5,2*w0 # Régime critique
    Q,Wc=10,2*w0 # Régime sinusoïdal
    D=1+1j*Q*(w/Wc-Wc/(w+0.001))
    return 1/D

def s_(t):
    s=0
    for n in range(Nmax+1):
        s+=CN[n]*H_(n*w0)*np.exp(1j*n*w0*t)
    return s
def s(t):
    return np.real(s_(t))

def courbe1():# Représentation graphique de e(t)
    plt.close()
    X=[(t-500)*T0/1000 for t in range(1200)]
    Y=[e(t) for t in X]
    plt.plot(X,Y)
    plt.title("Addition des {:.0f} premiers             harmoniques".format(Nmax+1))
    plt.grid()
    plt.show()

def courbe2():# Représentation graphique de e(t) et s(t)
    plt.close()
    X=[(t-500)*T0/1000 for t in range(1200)]
    Y1=[e(t) for t in X]
    Y2=[s(t) for t in X]
    plt.plot(X,Y2)
    plt.plot(X,Y1)
    plt.title("Addition des {:.0f} premiers             harmoniques".format(Nmax+1))
    plt.grid()
    plt.show()