1#!/usr/bin/env python3
2"""
3Clemmow-Mullaly-Allis (CMA) Diagram
4
5This script generates the CMA diagram showing wave propagation regions
6in a cold magnetized plasma for parallel propagation (k || B).
7
8Key Features:
9- Axes: X = ωpe²/ω² (horizontal), Y = ωce/ω (vertical)
10- Cutoff curves: R=0, L=0, P=0
11- Resonance curves: S=0 (upper/lower hybrid)
12- 13 distinct propagation regions
13- Color-coded by number of propagating modes
14
15Reference: Stix, "Waves in Plasmas" (1992)
16
17Author: Plasma Physics Examples
18License: MIT
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23from matplotlib.patches import Polygon
24from matplotlib.collections import PatchCollection
25
26def compute_stix_parameters(X, Y, Z=0):
27 """
28 Compute Stix parameters for cold plasma dispersion.
29
30 Parameters:
31 -----------
32 X : float or array
33 X = ωpe²/ω²
34 Y : float or array
35 Y = ωce/ω (electron cyclotron)
36 Z : float or array
37 Z = ωci/ω (ion cyclotron, usually neglected)
38
39 Returns:
40 --------
41 S, D, P : Stix parameters
42 """
43 # Stix parameters
44 S = 1 - X / (1 - Y**2)
45 D = -X * Y / (1 - Y**2)
46 P = 1 - X
47
48 return S, D, P
49
50def R_cutoff(Y):
51 """R-mode cutoff: R = S + D = 0."""
52 # R = 1 - X/(1-Y) = 0
53 # X = 1 - Y
54 return 1 - Y
55
56def L_cutoff(Y):
57 """L-mode cutoff: L = S - D = 0."""
58 # L = 1 - X/(1+Y) = 0
59 # X = 1 + Y
60 return 1 + Y
61
62def P_cutoff(Y):
63 """P-mode (plasma) cutoff: P = 0."""
64 # P = 1 - X = 0
65 # X = 1
66 return np.ones_like(Y)
67
68def upper_hybrid_resonance(Y):
69 """Upper hybrid resonance: S = 0, ω² = ωpe² + ωce²."""
70 # S = 1 - X/(1-Y²) = 0
71 # X = 1 - Y²
72 return 1 - Y**2
73
74def lower_hybrid_resonance(Y):
75 """
76 Lower hybrid resonance (approximate).
77
78 Exact formula involves ion cyclotron frequency.
79 Here we use simplified version.
80 """
81 # For the CMA diagram, this is the ion resonance line
82 # In practice, this would be Y → -Ωi/ω ≈ 0 for small ion effects
83 # We'll mark it as a separate region
84 return 0 * Y # Placeholder
85
86def plot_cma_diagram():
87 """
88 Create the Clemmow-Mullaly-Allis (CMA) diagram.
89 """
90 # Create grid
91 Y_array = np.linspace(-3, 3, 1000)
92
93 # Compute cutoff and resonance curves
94 X_R = R_cutoff(Y_array)
95 X_L = L_cutoff(Y_array)
96 X_P = P_cutoff(Y_array)
97 X_UH = upper_hybrid_resonance(Y_array)
98
99 # Clip negative values
100 X_R = np.maximum(X_R, 0)
101 X_L = np.maximum(X_L, 0)
102 X_UH = np.maximum(X_UH, 0)
103
104 # Create figure
105 fig, ax = plt.subplots(figsize=(12, 10))
106
107 # Plot cutoff curves
108 ax.plot(X_R, Y_array, 'r-', linewidth=2.5, label='R-cutoff (R=0)')
109 ax.plot(X_L, Y_array, 'b-', linewidth=2.5, label='L-cutoff (L=0)')
110 ax.plot(X_P, Y_array, 'g-', linewidth=2.5, label='P-cutoff (P=0)')
111
112 # Plot resonance curves
113 ax.plot(X_UH, Y_array, 'm--', linewidth=2.5, label='Upper Hybrid (S=0)')
114
115 # Add important lines
116 ax.axhline(y=0, color='black', linewidth=1.5, linestyle='-')
117 ax.axvline(x=0, color='black', linewidth=1.5, linestyle='-')
118 ax.axhline(y=1, color='orange', linewidth=1, linestyle=':', alpha=0.5)
119 ax.axhline(y=-1, color='orange', linewidth=1, linestyle=':', alpha=0.5)
120
121 # Define propagation regions and color them
122 regions = [
123 # Format: [X_range, Y_range, num_modes, label, position]
124 ([0, 0.5], [1.5, 2.5], 2, 'R+L', (0.25, 2.0)),
125 ([0, 0.5], [0.5, 1.5], 3, 'R+L+P', (0.25, 1.0)),
126 ([0.5, 1.0], [0.5, 1.5], 2, 'L+P', (0.75, 1.0)),
127 ([0.5, 1.0], [1.5, 2.5], 1, 'L', (0.75, 2.0)),
128 ([1.0, 2.0], [0.5, 1.5], 1, 'L', (1.5, 1.0)),
129 ([1.0, 2.0], [1.5, 2.5], 0, 'None', (1.5, 2.0)),
130 ([0, 0.5], [-1.5, 0.5], 3, 'R+L+P', (0.25, -0.5)),
131 ([0.5, 1.0], [-1.5, 0.5], 2, 'R+P', (0.75, -0.5)),
132 ([1.0, 2.0], [-1.5, 0.5], 1, 'R', (1.5, -0.5)),
133 ([0, 0.5], [-2.5, -1.5], 2, 'R+L', (0.25, -2.0)),
134 ([0.5, 1.0], [-2.5, -1.5], 1, 'R', (0.75, -2.0)),
135 ([1.0, 2.0], [-2.5, -1.5], 0, 'None', (1.5, -2.0)),
136 ]
137
138 # Color map based on number of modes
139 colors = {0: '#FFCCCC', 1: '#FFFFCC', 2: '#CCFFCC', 3: '#CCFFFF'}
140
141 for X_range, Y_range, num_modes, label, pos in regions:
142 # Create rectangle
143 rect = plt.Rectangle((X_range[0], Y_range[0]),
144 X_range[1] - X_range[0],
145 Y_range[1] - Y_range[0],
146 facecolor=colors[num_modes],
147 edgecolor='gray',
148 alpha=0.3,
149 linewidth=0.5)
150 ax.add_patch(rect)
151
152 # Add label
153 ax.text(pos[0], pos[1], f'{label}\n({num_modes} modes)',
154 ha='center', va='center', fontsize=9,
155 bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))
156
157 # Mark special points
158 # Point A: R-L crossover at Y=0
159 ax.plot(1, 0, 'ko', markersize=8)
160 ax.text(1.05, 0.1, 'A: ω=ωpe', fontsize=9)
161
162 # Point B: Upper hybrid at Y=0
163 ax.plot(1, 0, 'mo', markersize=8)
164
165 # Point C: Upper hybrid resonance example
166 Y_example = 0.5
167 X_example = 1 - Y_example**2
168 ax.plot(X_example, Y_example, 'mo', markersize=8)
169 ax.text(X_example + 0.05, Y_example + 0.1,
170 f'C: UH\nω²=ωpe²+ωce²', fontsize=9)
171
172 # Add annotations for important regimes
173 # Ionosphere
174 ax.annotate('Ionosphere\n(X>>1, Y<<1)',
175 xy=(3, 0.3), fontsize=10,
176 bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))
177
178 # Magnetosphere
179 ax.annotate('Magnetosphere\n(X<<1, Y~1)',
180 xy=(0.3, 1.5), fontsize=10,
181 bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.5))
182
183 # Tokamak core
184 ax.annotate('Tokamak core\n(X~1, Y~0.1)',
185 xy=(1.3, 0.15), fontsize=10,
186 bbox=dict(boxstyle='round', facecolor='pink', alpha=0.5))
187
188 # Labels and formatting
189 ax.set_xlabel(r'$X = \omega_{pe}^2 / \omega^2$', fontsize=14, fontweight='bold')
190 ax.set_ylabel(r'$Y = \omega_{ce} / \omega$', fontsize=14, fontweight='bold')
191 ax.set_title('Clemmow-Mullaly-Allis (CMA) Diagram\nParallel Propagation (k || B)',
192 fontsize=16, fontweight='bold')
193
194 ax.set_xlim([0, 3.5])
195 ax.set_ylim([-3, 3])
196 ax.grid(True, alpha=0.3, linestyle='--')
197 ax.legend(loc='upper right', fontsize=11, framealpha=0.9)
198
199 # Add text box with explanation
200 textstr = '\n'.join([
201 'Cutoffs: n² → 0',
202 'Resonances: n² → ∞',
203 '',
204 'Modes:',
205 'R: Right-hand circular',
206 'L: Left-hand circular',
207 'P: Plasma (O-mode)',
208 ])
209 ax.text(2.5, -2.3, textstr, fontsize=9,
210 verticalalignment='top',
211 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
212
213 plt.tight_layout()
214 plt.savefig('cma_diagram.png', dpi=200, bbox_inches='tight')
215 print("CMA diagram saved as 'cma_diagram.png'")
216 print("\n" + "="*70)
217 print("CMA Diagram Information")
218 print("="*70)
219 print("Axes:")
220 print(" X = ωpe²/ω² (plasma parameter)")
221 print(" Y = ωce/ω (magnetic parameter)")
222 print("\nCutoff curves (n² = 0, wave reflects):")
223 print(" R-cutoff: X = 1 - Y")
224 print(" L-cutoff: X = 1 + Y")
225 print(" P-cutoff: X = 1")
226 print("\nResonance curves (n² → ∞, wave absorbed):")
227 print(" Upper Hybrid: X = 1 - Y²")
228 print("\nPropagating modes:")
229 print(" R: Right-hand polarized (electron cyclotron)")
230 print(" L: Left-hand polarized (whistler, helicon)")
231 print(" P: Ordinary mode (electric field || B)")
232 print("="*70)
233
234 plt.show()
235
236 # Create a second figure showing specific wave types
237 plot_wave_trajectories()
238
239def plot_wave_trajectories():
240 """
241 Plot example wave trajectories through CMA diagram.
242 """
243 fig, ax = plt.subplots(figsize=(12, 8))
244
245 # Background
246 Y_array = np.linspace(-2, 2, 1000)
247 X_R = np.maximum(R_cutoff(Y_array), 0)
248 X_L = np.maximum(L_cutoff(Y_array), 0)
249 X_P = P_cutoff(Y_array)
250 X_UH = np.maximum(upper_hybrid_resonance(Y_array), 0)
251
252 ax.plot(X_R, Y_array, 'r-', linewidth=2, alpha=0.5, label='R-cutoff')
253 ax.plot(X_L, Y_array, 'b-', linewidth=2, alpha=0.5, label='L-cutoff')
254 ax.plot(X_P, Y_array, 'g-', linewidth=2, alpha=0.5, label='P-cutoff')
255 ax.plot(X_UH, Y_array, 'm--', linewidth=2, alpha=0.5, label='UH resonance')
256
257 # Example trajectories for fixed frequency, varying density
258 # Path 1: Low frequency wave entering plasma (Y = 2)
259 Y1 = 2.0
260 X_path1 = np.linspace(0, 3, 100)
261 Y_path1 = Y1 * np.ones_like(X_path1)
262
263 ax.plot(X_path1, Y_path1, 'purple', linewidth=3, label='Low freq. (Y=2)')
264 ax.annotate('', xy=(3, Y1), xytext=(0, Y1),
265 arrowprops=dict(arrowstyle='->', lw=2, color='purple'))
266
267 # Mark cutoff
268 X_cutoff1 = R_cutoff(Y1)
269 ax.plot(X_cutoff1, Y1, 'ro', markersize=10)
270 ax.text(X_cutoff1, Y1 + 0.15, 'R-cutoff', fontsize=10)
271
272 # Path 2: Intermediate frequency (Y = 0.5)
273 Y2 = 0.5
274 X_path2 = np.linspace(0, 3, 100)
275 Y_path2 = Y2 * np.ones_like(X_path2)
276
277 ax.plot(X_path2, Y_path2, 'orange', linewidth=3, label='Mid freq. (Y=0.5)')
278 ax.annotate('', xy=(3, Y2), xytext=(0, Y2),
279 arrowprops=dict(arrowstyle='->', lw=2, color='orange'))
280
281 # Mark UH resonance
282 X_UH2 = upper_hybrid_resonance(Y2)
283 ax.plot(X_UH2, Y2, 'mo', markersize=10)
284 ax.text(X_UH2, Y2 + 0.15, 'UH', fontsize=10)
285
286 # Path 3: Unmagnetized (Y = 0)
287 Y3 = 0.0
288 X_path3 = np.linspace(0, 3, 100)
289 Y_path3 = Y3 * np.ones_like(X_path3)
290
291 ax.plot(X_path3, Y_path3, 'green', linewidth=3, label='Unmagnetized (Y=0)')
292 ax.annotate('', xy=(3, Y3), xytext=(0, Y3),
293 arrowprops=dict(arrowstyle='->', lw=2, color='green'))
294
295 # Mark plasma cutoff
296 ax.plot(1, Y3, 'go', markersize=10)
297 ax.text(1, Y3 - 0.15, 'Plasma cutoff', fontsize=10, va='top')
298
299 ax.set_xlabel(r'$X = \omega_{pe}^2 / \omega^2$', fontsize=14, fontweight='bold')
300 ax.set_ylabel(r'$Y = \omega_{ce} / \omega$', fontsize=14, fontweight='bold')
301 ax.set_title('Wave Trajectories in CMA Space\n(Increasing density at fixed frequency)',
302 fontsize=14, fontweight='bold')
303 ax.set_xlim([0, 3.5])
304 ax.set_ylim([-2, 2.5])
305 ax.grid(True, alpha=0.3)
306 ax.legend(loc='upper right', fontsize=10)
307 ax.axhline(y=0, color='black', linewidth=1)
308 ax.axvline(x=0, color='black', linewidth=1)
309
310 plt.tight_layout()
311 plt.savefig('cma_wave_paths.png', dpi=150, bbox_inches='tight')
312 print("Wave trajectory plot saved as 'cma_wave_paths.png'")
313
314 plt.show()
315
316if __name__ == "__main__":
317 plot_cma_diagram()