1#!/usr/bin/env python3
2"""
3General Plasma Wave Dispersion Relation Solver
4
5This script solves the cold plasma dispersion relation for arbitrary
6propagation angles and plasma conditions.
7
8Features:
9- Cold plasma dielectric tensor (Stix parameters S, D, P)
10- Solve det(wave equation) = 0 for ω(k) or k(ω)
11- Wave modes: R, L, O, X for any angle θ
12- Generate ω-k diagrams and CMA diagrams
13- Multiple ion species support
14
15Author: Plasma Physics Examples
16License: MIT
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from matplotlib.gridspec import GridSpec
22from scipy.optimize import fsolve, brentq
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
31class ColdPlasmaDispersion:
32 """Cold plasma dispersion relation solver."""
33
34 def __init__(self, ne, B0, ion_species='H'):
35 """
36 Initialize plasma parameters.
37
38 Parameters:
39 -----------
40 ne : float
41 Electron density [m^-3]
42 B0 : float
43 Magnetic field [T]
44 ion_species : str
45 Ion species ('H', 'D', 'He')
46 """
47 self.ne = ne
48 self.B0 = B0
49
50 # Ion mass
51 if ion_species == 'H':
52 self.mi = MP
53 elif ion_species == 'D':
54 self.mi = 2 * MP
55 elif ion_species == 'He':
56 self.mi = 4 * MP
57 else:
58 self.mi = MP
59
60 # Compute characteristic frequencies
61 self.omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
62 self.omega_pi = np.sqrt(ne * QE**2 / (self.mi * EPS0))
63 self.omega_ce = QE * B0 / ME
64 self.omega_ci = QE * B0 / self.mi
65
66 self.f_pe = self.omega_pe / (2 * np.pi)
67 self.f_ce = self.omega_ce / (2 * np.pi)
68 self.f_ci = self.omega_ci / (2 * np.pi)
69
70 def stix_parameters(self, omega):
71 """
72 Compute Stix parameters S, D, P.
73
74 Parameters:
75 -----------
76 omega : float or array
77 Angular frequency [rad/s]
78
79 Returns:
80 --------
81 S, D, P : Stix parameters
82 """
83 # Avoid division by zero
84 eps = 1e-10
85
86 # S = 1 - Σ ωps²/(ω² - ωcs²)
87 S = 1 - self.omega_pe**2 / (omega**2 - self.omega_ce**2 + eps)
88 S -= self.omega_pi**2 / (omega**2 - self.omega_ci**2 + eps)
89
90 # D = Σ (ωcs/ω) ωps²/(ω² - ωcs²)
91 D = (self.omega_ce / omega) * self.omega_pe**2 / (omega**2 - self.omega_ce**2 + eps)
92 D += (self.omega_ci / omega) * self.omega_pi**2 / (omega**2 - self.omega_ci**2 + eps)
93
94 # P = 1 - Σ ωps²/ω²
95 P = 1 - self.omega_pe**2 / omega**2 - self.omega_pi**2 / omega**2
96
97 return S, D, P
98
99 def dispersion_general(self, omega, k, theta):
100 """
101 General dispersion relation for arbitrary angle θ.
102
103 det(M) = An⁴ - Bn² + C = 0
104
105 where n = ck/ω
106
107 Parameters:
108 -----------
109 omega : float
110 Angular frequency [rad/s]
111 k : float
112 Wavenumber [rad/m]
113 theta : float
114 Propagation angle w.r.t. B [rad]
115
116 Returns:
117 --------
118 residual : float (should be zero for valid solutions)
119 """
120 S, D, P = self.stix_parameters(omega)
121
122 n = C * k / omega # Refractive index
123
124 cos_theta = np.cos(theta)
125 sin_theta = np.sin(theta)
126
127 # Dispersion relation coefficients
128 A = S * sin_theta**2 + P * cos_theta**2
129 B = (S**2 - D**2) * sin_theta**2 + P * S * (1 + cos_theta**2)
130 C_coeff = P * (S**2 - D**2)
131
132 # Dispersion relation: An⁴ - Bn² + C = 0
133 residual = A * n**4 - B * n**2 + C_coeff
134
135 return residual
136
137 def parallel_modes(self, omega):
138 """
139 Solve for parallel propagation (θ = 0).
140
141 R-wave: n² = S + D = R
142 L-wave: n² = S - D = L
143
144 Parameters:
145 -----------
146 omega : float or array
147 Angular frequency [rad/s]
148
149 Returns:
150 --------
151 n_R, n_L : refractive indices for R and L modes
152 """
153 S, D, P = self.stix_parameters(omega)
154
155 R = S + D
156 L = S - D
157
158 n_R = np.sqrt(np.maximum(R, 0))
159 n_L = np.sqrt(np.maximum(L, 0))
160
161 return n_R, n_L
162
163 def perpendicular_modes(self, omega):
164 """
165 Solve for perpendicular propagation (θ = 90°).
166
167 O-mode: n² = P
168 X-mode: n² = (S² - D²) / S
169
170 Parameters:
171 -----------
172 omega : float or array
173 Angular frequency [rad/s]
174
175 Returns:
176 --------
177 n_O, n_X : refractive indices for O and X modes
178 """
179 S, D, P = self.stix_parameters(omega)
180
181 n_O = np.sqrt(np.maximum(P, 0))
182 n_X_sq = (S**2 - D**2) / (S + 1e-10)
183 n_X = np.sqrt(np.maximum(n_X_sq, 0))
184
185 return n_O, n_X
186
187 def find_cutoffs_resonances(self):
188 """
189 Find cutoff and resonance frequencies.
190
191 Returns:
192 --------
193 dict : Cutoff and resonance frequencies
194 """
195 # Cutoffs (n² = 0)
196 # R cutoff: R = 0
197 # L cutoff: L = 0
198 # P cutoff: P = 0
199
200 # Resonances (n² → ∞)
201 # Upper hybrid: ω² = ωpe² + ωce²
202 # Lower hybrid: ω² = ωci·ωce + ωpi²/(1 + ωpe²/ωce²)
203
204 omega_uh = np.sqrt(self.omega_pe**2 + self.omega_ce**2)
205 omega_lh = np.sqrt(self.omega_ci * self.omega_ce)
206
207 return {
208 'f_uh': omega_uh / (2 * np.pi),
209 'f_lh': omega_lh / (2 * np.pi),
210 'f_pe': self.f_pe,
211 'f_ce': self.f_ce,
212 'f_ci': self.f_ci
213 }
214
215def plot_dispersion_solver():
216 """
217 Demonstrate dispersion relation solver with multiple plots.
218 """
219 # Plasma parameters
220 ne = 1e18 # m^-3
221 B0 = 1.0 # T
222 ion_species = 'H'
223
224 solver = ColdPlasmaDispersion(ne, B0, ion_species)
225
226 print("=" * 70)
227 print("Cold Plasma Dispersion Relation Solver")
228 print("=" * 70)
229 print(f"Electron density: {ne:.2e} m^-3")
230 print(f"Magnetic field: {B0:.2f} T")
231 print(f"Ion species: {ion_species}")
232 print(f"\nCharacteristic frequencies:")
233 print(f" Electron plasma: {solver.f_pe/1e9:.3f} GHz")
234 print(f" Electron cyclotron: {solver.f_ce/1e9:.3f} GHz")
235 print(f" Ion cyclotron: {solver.f_ci/1e6:.3f} MHz")
236
237 cutoffs = solver.find_cutoffs_resonances()
238 print(f"\nCutoffs and resonances:")
239 print(f" Upper hybrid: {cutoffs['f_uh']/1e9:.3f} GHz")
240 print(f" Lower hybrid: {cutoffs['f_lh']/1e6:.3f} MHz")
241 print("=" * 70)
242
243 # Create figure
244 fig = plt.figure(figsize=(16, 12))
245 gs = GridSpec(3, 2, figure=fig, hspace=0.35, wspace=0.3)
246
247 # Plot 1: ω-k diagram for parallel propagation
248 ax1 = fig.add_subplot(gs[0, :])
249
250 omega_array = np.linspace(0.1 * solver.omega_ce, 3 * solver.omega_ce, 1000)
251
252 n_R, n_L = solver.parallel_modes(omega_array)
253 k_R = omega_array * n_R / C
254 k_L = omega_array * n_L / C
255
256 # Light line
257 k_light = omega_array / C
258
259 ax1.plot(k_R / 1e6, omega_array / (2 * np.pi * 1e9), 'r-',
260 linewidth=2, label='R-mode')
261 ax1.plot(k_L / 1e6, omega_array / (2 * np.pi * 1e9), 'b-',
262 linewidth=2, label='L-mode')
263 ax1.plot(k_light / 1e6, omega_array / (2 * np.pi * 1e9), 'k--',
264 linewidth=1, label='Light line', alpha=0.5)
265
266 # Mark characteristic frequencies
267 ax1.axhline(y=solver.f_pe / 1e9, color='g', linestyle=':', linewidth=2,
268 label=f'fpe = {solver.f_pe/1e9:.2f} GHz')
269 ax1.axhline(y=solver.f_ce / 1e9, color='m', linestyle=':', linewidth=2,
270 label=f'fce = {solver.f_ce/1e9:.2f} GHz')
271 ax1.axhline(y=cutoffs['f_uh'] / 1e9, color='orange', linestyle=':', linewidth=2,
272 label=f'fUH = {cutoffs["f_uh"]/1e9:.2f} GHz')
273
274 ax1.set_xlabel('Wavenumber k (rad/Mm)', fontsize=12)
275 ax1.set_ylabel('Frequency (GHz)', fontsize=12)
276 ax1.set_title('ω-k Diagram: Parallel Propagation (θ = 0°)',
277 fontsize=13, fontweight='bold')
278 ax1.legend(fontsize=10, loc='lower right')
279 ax1.grid(True, alpha=0.3)
280 ax1.set_xlim([0, 50])
281 ax1.set_ylim([0, 100])
282
283 # Plot 2: ω-k diagram for perpendicular propagation
284 ax2 = fig.add_subplot(gs[1, 0])
285
286 n_O, n_X = solver.perpendicular_modes(omega_array)
287 k_O = omega_array * n_O / C
288 k_X = omega_array * n_X / C
289
290 ax2.plot(k_O / 1e6, omega_array / (2 * np.pi * 1e9), 'g-',
291 linewidth=2, label='O-mode')
292 ax2.plot(k_X / 1e6, omega_array / (2 * np.pi * 1e9), 'c-',
293 linewidth=2, label='X-mode')
294 ax2.plot(k_light / 1e6, omega_array / (2 * np.pi * 1e9), 'k--',
295 linewidth=1, label='Light line', alpha=0.5)
296
297 ax2.axhline(y=solver.f_pe / 1e9, color='purple', linestyle=':', linewidth=2)
298 ax2.axhline(y=cutoffs['f_uh'] / 1e9, color='orange', linestyle=':', linewidth=2)
299
300 ax2.set_xlabel('Wavenumber k (rad/Mm)', fontsize=11)
301 ax2.set_ylabel('Frequency (GHz)', fontsize=11)
302 ax2.set_title('ω-k Diagram: Perpendicular (θ = 90°)',
303 fontsize=12, fontweight='bold')
304 ax2.legend(fontsize=10)
305 ax2.grid(True, alpha=0.3)
306 ax2.set_xlim([0, 50])
307 ax2.set_ylim([0, 100])
308
309 # Plot 3: CMA diagram
310 ax3 = fig.add_subplot(gs[1, 1])
311
312 X_array = np.linspace(0, 3, 1000)
313 Y_array = np.linspace(-2, 2, 1000)
314
315 # Current plasma parameters
316 X_plasma = solver.omega_pe**2 / solver.omega_ce**2
317 Y_plasma = 1.0 # At ω = ωce
318
319 # Cutoff curves
320 X_R = 1 - Y_array # R cutoff
321 X_L = 1 + Y_array # L cutoff
322 X_P = np.ones_like(Y_array) # P cutoff
323
324 # Resonance curves
325 X_UH = 1 - Y_array**2 # Upper hybrid
326
327 ax3.plot(X_R, Y_array, 'r-', linewidth=2, label='R cutoff')
328 ax3.plot(X_L, Y_array, 'b-', linewidth=2, label='L cutoff')
329 ax3.plot(X_P, Y_array, 'g-', linewidth=2, label='P cutoff')
330 ax3.plot(X_UH, Y_array, 'm--', linewidth=2, label='UH resonance')
331
332 # Mark current plasma
333 ax3.plot(X_plasma, Y_plasma, 'ko', markersize=10,
334 label=f'This plasma\n(X={X_plasma:.2f}, Y=1)')
335
336 ax3.set_xlabel(r'$X = \omega_{pe}^2 / \omega^2$', fontsize=11)
337 ax3.set_ylabel(r'$Y = \omega_{ce} / \omega$', fontsize=11)
338 ax3.set_title('CMA Diagram', fontsize=12, fontweight='bold')
339 ax3.legend(fontsize=9)
340 ax3.grid(True, alpha=0.3)
341 ax3.set_xlim([0, 3])
342 ax3.set_ylim([-2, 2])
343 ax3.axhline(y=0, color='k', linewidth=0.5)
344 ax3.axvline(x=0, color='k', linewidth=0.5)
345
346 # Plot 4: Refractive index vs frequency (parallel)
347 ax4 = fig.add_subplot(gs[2, 0])
348
349 ax4.plot(omega_array / (2 * np.pi * 1e9), n_R, 'r-',
350 linewidth=2, label='R-mode')
351 ax4.plot(omega_array / (2 * np.pi * 1e9), n_L, 'b-',
352 linewidth=2, label='L-mode')
353
354 ax4.axhline(y=1, color='k', linestyle='--', linewidth=1,
355 label='Vacuum (n=1)', alpha=0.5)
356 ax4.axvline(x=solver.f_pe / 1e9, color='g', linestyle=':', linewidth=2)
357 ax4.axvline(x=solver.f_ce / 1e9, color='m', linestyle=':', linewidth=2)
358
359 ax4.set_xlabel('Frequency (GHz)', fontsize=11)
360 ax4.set_ylabel('Refractive Index n', fontsize=11)
361 ax4.set_title('Refractive Index: Parallel', fontsize=12, fontweight='bold')
362 ax4.legend(fontsize=10)
363 ax4.grid(True, alpha=0.3)
364 ax4.set_xlim([0, 100])
365 ax4.set_ylim([0, 10])
366
367 # Plot 5: Refractive index vs frequency (perpendicular)
368 ax5 = fig.add_subplot(gs[2, 1])
369
370 ax5.plot(omega_array / (2 * np.pi * 1e9), n_O, 'g-',
371 linewidth=2, label='O-mode')
372 ax5.plot(omega_array / (2 * np.pi * 1e9), n_X, 'c-',
373 linewidth=2, label='X-mode')
374
375 ax5.axhline(y=1, color='k', linestyle='--', linewidth=1, alpha=0.5)
376 ax5.axvline(x=solver.f_pe / 1e9, color='purple', linestyle=':', linewidth=2)
377 ax5.axvline(x=cutoffs['f_uh'] / 1e9, color='orange', linestyle=':', linewidth=2)
378
379 ax5.set_xlabel('Frequency (GHz)', fontsize=11)
380 ax5.set_ylabel('Refractive Index n', fontsize=11)
381 ax5.set_title('Refractive Index: Perpendicular', fontsize=12, fontweight='bold')
382 ax5.legend(fontsize=10)
383 ax5.grid(True, alpha=0.3)
384 ax5.set_xlim([0, 100])
385 ax5.set_ylim([0, 5])
386
387 plt.suptitle('Cold Plasma Dispersion Relation Solver',
388 fontsize=16, fontweight='bold', y=0.995)
389
390 plt.savefig('dispersion_solver.png', dpi=150, bbox_inches='tight')
391 print("\nPlot saved as 'dispersion_solver.png'")
392
393 plt.show()
394
395if __name__ == "__main__":
396 plot_dispersion_solver()