1#!/usr/bin/env python3
2"""
3Whistler Wave Dispersion and Spectrogram
4
5This script demonstrates whistler wave physics including:
6- R-mode dispersion in magnetized plasma
7- Whistler regime: ωci << ω << ωce
8- Frequency-dependent group velocity
9- Simulated whistler spectrogram showing falling tone
10
11Key Physics:
12- Dispersion: n² ≈ ωpe²/(ω·ωce) in whistler regime
13- Group velocity: vg ∝ sqrt(ω) → higher frequencies travel faster
14- Observed as falling tones in magnetospheric radio emissions
15
16Author: Plasma Physics Examples
17License: MIT
18"""
19
20import numpy as np
21import matplotlib.pyplot as plt
22from matplotlib.gridspec import GridSpec
23
24# Physical constants
25EPS0 = 8.854187817e-12 # F/m
26QE = 1.602176634e-19 # C
27ME = 9.10938356e-31 # kg
28MP = 1.672621898e-27 # kg
29C = 2.99792458e8 # m/s
30
31def compute_characteristic_frequencies(ne, B0, mi=MP):
32 """
33 Compute plasma and cyclotron frequencies.
34
35 Parameters:
36 -----------
37 ne : float
38 Electron density [m^-3]
39 B0 : float
40 Magnetic field strength [T]
41 mi : float
42 Ion mass [kg]
43
44 Returns:
45 --------
46 dict : Characteristic frequencies
47 """
48 omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
49 omega_ce = QE * B0 / ME
50 omega_ci = QE * B0 / mi
51
52 f_pe = omega_pe / (2 * np.pi)
53 f_ce = omega_ce / (2 * np.pi)
54 f_ci = omega_ci / (2 * np.pi)
55
56 return {
57 'omega_pe': omega_pe,
58 'omega_ce': omega_ce,
59 'omega_ci': omega_ci,
60 'f_pe': f_pe,
61 'f_ce': f_ce,
62 'f_ci': f_ci
63 }
64
65def whistler_dispersion_full(omega, ne, B0):
66 """
67 Full R-mode dispersion relation: n² = 1 - ωpe²/(ω(ω - ωce)).
68
69 Parameters:
70 -----------
71 omega : array
72 Angular frequency [rad/s]
73 ne : float
74 Electron density [m^-3]
75 B0 : float
76 Magnetic field [T]
77
78 Returns:
79 --------
80 k : array
81 Wavenumber [rad/m]
82 """
83 params = compute_characteristic_frequencies(ne, B0)
84 omega_pe = params['omega_pe']
85 omega_ce = params['omega_ce']
86
87 # R-mode dispersion
88 n_sq = 1 - omega_pe**2 / (omega * (omega - omega_ce))
89
90 # Only valid for n² > 0
91 n_sq = np.maximum(n_sq, 0)
92
93 n = np.sqrt(n_sq)
94 k = omega * n / C
95
96 return k
97
98def whistler_dispersion_approx(omega, ne, B0):
99 """
100 Approximate whistler dispersion: n² ≈ ωpe²/(ω·ωce).
101
102 Valid in range: ωci << ω << ωce
103
104 Parameters:
105 -----------
106 omega : array
107 Angular frequency [rad/s]
108 ne : float
109 Electron density [m^-3]
110 B0 : float
111 Magnetic field [T]
112
113 Returns:
114 --------
115 k : array
116 Wavenumber [rad/m]
117 """
118 params = compute_characteristic_frequencies(ne, B0)
119 omega_pe = params['omega_pe']
120 omega_ce = params['omega_ce']
121
122 # Whistler approximation
123 n_sq = omega_pe**2 / (omega * omega_ce)
124
125 n = np.sqrt(n_sq)
126 k = omega * n / C
127
128 return k
129
130def group_velocity_whistler(omega, ne, B0):
131 """
132 Compute whistler group velocity: vg = dω/dk.
133
134 In whistler regime: vg ∝ sqrt(ω)
135
136 Parameters:
137 -----------
138 omega : array
139 Angular frequency [rad/s]
140 ne : float
141 Electron density [m^-3]
142 B0 : float
143 Magnetic field [T]
144
145 Returns:
146 --------
147 vg : array
148 Group velocity [m/s]
149 """
150 # Use numerical derivative
151 k = whistler_dispersion_full(omega, ne, B0)
152
153 # Compute dω/dk
154 vg = np.gradient(omega, k)
155
156 return vg
157
158def simulate_whistler_spectrogram(ne, B0, distance, f_start, f_end, num_freqs=100):
159 """
160 Simulate whistler spectrogram showing dispersion.
161
162 Parameters:
163 -----------
164 ne : float
165 Electron density [m^-3]
166 B0 : float
167 Magnetic field [T]
168 distance : float
169 Propagation distance [m]
170 f_start : float
171 Starting frequency [Hz]
172 f_end : float
173 Ending frequency [Hz]
174 num_freqs : int
175 Number of frequency components
176
177 Returns:
178 --------
179 t, f, spectrogram : arrays for plotting
180 """
181 # Frequency array
182 freqs = np.linspace(f_start, f_end, num_freqs)
183 omega = 2 * np.pi * freqs
184
185 # Compute group velocities
186 vg = group_velocity_whistler(omega, ne, B0)
187
188 # Arrival time at distance
189 arrival_times = distance / vg
190
191 # Normalize time to start at 0
192 arrival_times = arrival_times - arrival_times.min()
193
194 # Create time grid for spectrogram
195 t_max = arrival_times.max() * 1.2
196 t_grid = np.linspace(0, t_max, 500)
197
198 # Create spectrogram (Gaussian packets)
199 spectrogram = np.zeros((num_freqs, len(t_grid)))
200
201 for i, (f, t_arrival) in enumerate(zip(freqs, arrival_times)):
202 # Gaussian pulse centered at arrival time
203 sigma = t_max / 50 # Pulse width
204 pulse = np.exp(-(t_grid - t_arrival)**2 / (2 * sigma**2))
205 spectrogram[i, :] = pulse
206
207 return t_grid, freqs, spectrogram, arrival_times
208
209def plot_whistler_waves():
210 """
211 Create comprehensive visualization of whistler wave properties.
212 """
213 # Magnetospheric parameters (Earth's plasmasphere)
214 ne = 1e7 # m^-3 (10 cm^-3)
215 B0 = 1e-6 # Tesla (10 nT)
216 distance = 1e7 # 10,000 km propagation
217
218 params = compute_characteristic_frequencies(ne, B0)
219
220 print("=" * 70)
221 print("Whistler Wave Parameters (Magnetospheric)")
222 print("=" * 70)
223 print(f"Electron density: {ne:.2e} m^-3")
224 print(f"Magnetic field: {B0*1e9:.2f} nT")
225 print(f"Electron plasma frequency: {params['f_pe']/1e3:.2f} kHz")
226 print(f"Electron cyclotron frequency: {params['f_ce']/1e3:.2f} kHz")
227 print(f"Ion cyclotron frequency: {params['f_ci']:.2f} Hz")
228 print(f"Whistler regime: {params['f_ci']:.1f} Hz << f << {params['f_ce']/1e3:.1f} kHz")
229 print("=" * 70)
230
231 # Create figure
232 fig = plt.figure(figsize=(14, 12))
233 gs = GridSpec(3, 2, figure=fig, hspace=0.35, wspace=0.3)
234
235 # Frequency range for whistler regime
236 f_min = 100 # Hz
237 f_max = 0.5 * params['f_ce'] / (2 * np.pi) # Half electron cyclotron
238 omega = 2 * np.pi * np.linspace(f_min, f_max, 1000)
239 freqs = omega / (2 * np.pi)
240
241 # Compute dispersion relations
242 k_full = whistler_dispersion_full(omega, ne, B0)
243 k_approx = whistler_dispersion_approx(omega, ne, B0)
244
245 # Plot 1: Dispersion relation ω(k)
246 ax1 = fig.add_subplot(gs[0, :])
247
248 ax1.plot(k_full / 1e3, freqs / 1e3,
249 'b-', linewidth=2, label='Full R-mode')
250 ax1.plot(k_approx / 1e3, freqs / 1e3,
251 'r--', linewidth=2, label='Whistler approximation')
252
253 # Mark characteristic frequencies
254 ax1.axhline(y=params['f_ce'] / (1e3 * 2 * np.pi), color='green',
255 linestyle=':', linewidth=2, label=r'$f_{ce}$')
256 ax1.axhline(y=params['f_ci'] / 1e3, color='purple',
257 linestyle=':', linewidth=2, label=r'$f_{ci}$')
258
259 ax1.set_xlabel('Wavenumber k (rad/km)', fontsize=12)
260 ax1.set_ylabel('Frequency (kHz)', fontsize=12)
261 ax1.set_title('Whistler Wave Dispersion Relation', fontsize=14, fontweight='bold')
262 ax1.grid(True, alpha=0.3)
263 ax1.legend(fontsize=10)
264 ax1.set_xlim([0, np.max(k_full) / 1e3])
265
266 # Plot 2: Refractive index
267 ax2 = fig.add_subplot(gs[1, 0])
268
269 n_full = k_full * C / omega
270 n_approx = k_approx * C / omega
271
272 ax2.plot(freqs / 1e3, n_full,
273 'b-', linewidth=2, label='Full')
274 ax2.plot(freqs / 1e3, n_approx,
275 'r--', linewidth=2, label='Approximation')
276
277 ax2.set_xlabel('Frequency (kHz)', fontsize=12)
278 ax2.set_ylabel('Refractive Index n', fontsize=12)
279 ax2.set_title('Refractive Index', fontsize=12, fontweight='bold')
280 ax2.grid(True, alpha=0.3)
281 ax2.legend(fontsize=10)
282 ax2.set_yscale('log')
283
284 # Plot 3: Group velocity
285 ax3 = fig.add_subplot(gs[1, 1])
286
287 vg = group_velocity_whistler(omega, ne, B0)
288
289 ax3.plot(freqs / 1e3, vg / 1e6,
290 'b-', linewidth=2)
291
292 # Theoretical vg ∝ sqrt(ω) in whistler regime
293 vg_theory = 2 * C * omega / (params['omega_pe'] / np.sqrt(params['omega_ce'] / omega))
294 ax3.plot(freqs / 1e3, vg_theory / 1e6,
295 'r--', linewidth=2, label=r'$v_g \propto \sqrt{\omega}$')
296
297 ax3.set_xlabel('Frequency (kHz)', fontsize=12)
298 ax3.set_ylabel(r'Group Velocity (10$^6$ m/s)', fontsize=12)
299 ax3.set_title('Group Velocity (Higher f travels faster)',
300 fontsize=12, fontweight='bold')
301 ax3.grid(True, alpha=0.3)
302 ax3.legend(fontsize=10)
303
304 # Plot 4: Phase velocity
305 ax4 = fig.add_subplot(gs[2, 0])
306
307 vph = omega / k_full
308
309 ax4.plot(freqs / 1e3, vph / 1e6,
310 'b-', linewidth=2)
311
312 ax4.set_xlabel('Frequency (kHz)', fontsize=12)
313 ax4.set_ylabel(r'Phase Velocity (10$^6$ m/s)', fontsize=12)
314 ax4.set_title('Phase Velocity', fontsize=12, fontweight='bold')
315 ax4.grid(True, alpha=0.3)
316
317 # Plot 5: Whistler spectrogram
318 ax5 = fig.add_subplot(gs[2, 1])
319
320 f_start = 5e3 # 5 kHz
321 f_end = 15e3 # 15 kHz
322
323 t_grid, freqs_spec, spectrogram, arrival_times = \
324 simulate_whistler_spectrogram(ne, B0, distance, f_start, f_end, num_freqs=100)
325
326 # Plot spectrogram
327 im = ax5.pcolormesh(t_grid * 1e3, freqs_spec / 1e3, spectrogram,
328 shading='auto', cmap='hot')
329
330 # Overlay arrival time curve
331 ax5.plot(arrival_times * 1e3, freqs_spec / 1e3,
332 'c-', linewidth=2, label='Arrival time')
333
334 ax5.set_xlabel('Time (ms)', fontsize=12)
335 ax5.set_ylabel('Frequency (kHz)', fontsize=12)
336 ax5.set_title(f'Whistler Spectrogram (Distance = {distance/1e6:.0f} km)',
337 fontsize=12, fontweight='bold')
338 ax5.legend(fontsize=9, loc='upper right')
339
340 # Add colorbar
341 cbar = plt.colorbar(im, ax=ax5)
342 cbar.set_label('Intensity', fontsize=10)
343
344 plt.suptitle('Whistler Waves in Magnetized Plasma',
345 fontsize=16, fontweight='bold', y=0.997)
346
347 plt.savefig('whistler_wave.png', dpi=150, bbox_inches='tight')
348 print("\nPlot saved as 'whistler_wave.png'")
349 print(f"\nWhistler shows 'falling tone': high frequencies arrive first")
350 print(f"Frequency drop from {f_end/1e3:.1f} to {f_start/1e3:.1f} kHz")
351 print(f"over time span of {(arrival_times.max() - arrival_times.min())*1e3:.2f} ms")
352
353 plt.show()
354
355if __name__ == "__main__":
356 plot_whistler_waves()