paquet_donde.py

Retour Ă  la liste des codes python.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#!/usr/env python
# -*- coding: utf-8 -*-
##############################################################################
import imgen
imgen.setup()

##############################################################################
import numpy as np
######
# Paramètres du problème
######
t_min = 0.
t_max = 4.

c = 300. # Mm/s, vitesse de la lumière
w_p_disp = 60. # MHz, pulsation plasma (qq MHz pour l'ionosphère)
w_p_nondisp = 10.

# Spectre du signal
w0 = 80.  # MHz, pulsation centrale
dw = 10.  # largeur de la gaussienne en pulsation
Nw = 1000
w_min = max(0, w0 - 10*dw)
w_max = w0 + 10*dw
W = np.linspace(w_min, w_max, Nw)  # domaine de pulsations
def A(w):
    return np.exp(- ((w - w0)**2) / (2 * dw**2))  # spectre gaussien

# FenĂŞtre spatiale
x_min = -0.  # m
x_max = c * t_max  # m
Nx = 1000
X = np.linspace(x_min, x_max, Nx)

######
# Calcul de la relaiton de dispersion
######
def K(w, w_p):
    k1 = np.sqrt((0.j + w**2 - w_p**2)) / c
    return np.conjugate(k1)  # On doit prendre les parties imaginaires négatives

######
# Calcul du signal S(X, t)
######
def S(t, w_p):
    S = np.zeros_like(X, dtype=np.complex)
    for i in range(len(W)):
        w = W[i]
        k = K(w, w_p)
        a = A(w)
        S += a * np.exp(1.j * (w*t - k*X))
    return S

######
# Calcul de c*t
######
def C(t):
    return c*t

##############################################################################
import matplotlib.pyplot as plt
# Affichage

T = np.linspace(t_min, t_max, 5)
fig, axes = plt.subplots(len(T) + 1, 2)

fig.suptitle(u"Propagation d'un paquet d'ondes")
axes[0,0].set_title(u"non dispersif $\omega_0 = %i \gg \omega_p = %i$" % (w0, w_p_nondisp))
axes[0,1].set_title(u"dispersif $\omega_0 = %i > \omega_p = %i$" % (w0, w_p_disp))

for n in range(len(T)):
    t = T[n]
    # Onde
    axes[n,0].plot(X, S(t, w_p_nondisp).real)
    axes[n,1].plot(X, S(t,    w_p_disp).real)
    # Vitesse c
    axes[n,0].plot(C(t), 0, marker='o', color='r')
    axes[n,1].plot(C(t), 0, marker='o', color='r')
    # Labels
    axes[n,0].set_ylabel("t = %3.f" % t)
    axes[n,0].set_yticks([])
    axes[n,1].set_yticks([])
    # Limites
    axes[n,0].set_xlim(x_min, x_max)
    axes[n,1].set_xlim(x_min, x_max)

# Labels
axes[-2,0].set_xlabel(u"x")
axes[-2,1].set_xlabel(u"x")

# Relation de dispersion et TF du signal
axes[-1,0].plot(W, K(W, w_p_nondisp).real, color='orange')
axes[-1,0].plot(W, A(W).real, color='olive')
axes[-1,0].legend([u"Relation de dispersion", u"TF du signal"])
axes[-1,0].autoscale(enable=True)

axes[-1,1].plot(W, K(W, w_p_disp).real, color='orange')
axes[-1,1].plot(W, A(W).real, color='olive')
axes[-1,1].legend([u"Relation de dispersion", u"TF du signal"])
axes[-1,1].autoscale(enable=True)


##############################################################################
# Save
imgen.done(__file__)