# Capacité 4 : PROPAGATION D'UN PAQUET D'ONDES
#############################################################################################################
# Exemple  5 : PROPAGATION D'UN PAQUETS D'ONDES GAUSSIEN DANS UN MILIEU DISPERSIF
#              Relation de dispersion développée à l'ordre 2 en w - w0
#              Utilisation de la fonction animation   

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Spectre en pulsation w du paquet d'ondes gaussien 

def A(w):
    return np.exp(-(w-w0)**2/(2*sigma**2))
    
wmin = 0                # Pulsation minimale
wmax = 10               # Pulsation maximale
w0 = 5                  # Pulsation centrale
sigma = 0.35            # Paramètre caractéristique de la largeur du spectre

N = 1000                # Choix du nombre de valeurs de w           

w = np.linspace(wmin,wmax,N)    # Plage de variation de w

plt.plot(w,A(w))
plt.xlabel('pulsation w')
plt.ylabel('amplitude',rotation=0)  
plt.title("Spectre du paquet d'ondes")
plt.grid()
plt.show()
plt.pause(5)            # Temps d'affichage de la courbe  
plt.clf()               # Effacement de la courbe

# Expression du paquet d'ondes

def paquet(x,t):
    sum = 0
    for w in np.linspace(wmin,wmax,N):
        sum = sum + A(w)*np.cos(w*t-k(w)*x)
    return sum

# Relation de dispersion développée à l'ordre 2 en w-w0 

def k(w):
    return k0 + dk0*(w-w0) + d2k0*(w-w0)**2/2   # k(w) = k(w0) + (dk/dw)(w0)(w-w0) + (d²k/dw²)(w0)(w-w0)²/2
    
k0 = 5                         # k(w0)
dk0 = 1.5                      # (dk/dw)(w0)
d2k0 = 0.5                     # (d²k/dw²)(w0)
    
x=np.linspace(-10,100,1000)

# Tracé du paquet d'ondes à t = 0

plt.axis([-10,100,-100,100])
plt.xlabel('abscisse x')
plt.ylabel('s(x,0)',rotation=0)  
plt.title("Paquet d'ondes s(x,0)")
plt.plot(x,paquet(x,0),'b')
plt.grid()
plt.show()
plt.pause(5)   

# Propagation du paquet d'ondes

dt = 0.1                  # Pas de temps

fig = plt.figure()        # Nom de la fenêtre graphique contenant la courbe à animer        
plt.xlim(-10, 100)
plt.ylim(-100,100)
plt.grid()

def init():               # Initialise la figure
    return () 
    
def animate(i): 
    t = i * dt
    return (plt.plot(x,paquet(x,t),'b'))

# Il y a 1000 images (frames) tracées avec un pas dt = 0,1 s donc la plage de temps du tracé est 100 s
ani = animation.FuncAnimation(fig,animate, init_func=init, frames=1000, blit=True, interval=100, repeat=False)

plt.show()
##############################################################################################################