1#!/usr/bin/env python3
2"""
3Plasma Parameters Calculator and Visualization
4
5This script computes fundamental plasma parameters for various plasma regimes
6and visualizes the parameter space on a density-temperature diagram.
7
8Author: Plasma Physics Examples
9License: MIT
10"""
11
12import numpy as np
13import matplotlib.pyplot as plt
14from scipy.constants import e, m_e, m_p, epsilon_0, k, c, mu_0
15
16# Physical constants
17m_i = m_p # Assume proton mass for ions
18
19
20def compute_plasma_parameters(n_e, T_e, B=0.0, Z=1):
21 """
22 Compute fundamental plasma parameters.
23
24 Parameters:
25 -----------
26 n_e : float
27 Electron density [m^-3]
28 T_e : float
29 Electron temperature [eV]
30 B : float
31 Magnetic field strength [Tesla]
32 Z : int
33 Ion charge number
34
35 Returns:
36 --------
37 dict : Dictionary containing all plasma parameters
38 """
39 # Convert temperature to Joules
40 T_e_J = T_e * e
41
42 # Debye length [m]
43 lambda_D = np.sqrt(epsilon_0 * T_e_J / (n_e * e**2))
44
45 # Plasma frequency [rad/s]
46 omega_pe = np.sqrt(n_e * e**2 / (epsilon_0 * m_e))
47 omega_pi = np.sqrt(n_e * Z**2 * e**2 / (epsilon_0 * m_i))
48
49 # Thermal velocity [m/s]
50 v_th_e = np.sqrt(2 * T_e_J / m_e)
51 v_th_i = np.sqrt(2 * T_e_J / m_i)
52
53 # Magnetic field parameters (if B > 0)
54 if B > 0:
55 omega_ce = e * B / m_e # Electron gyrofrequency [rad/s]
56 omega_ci = Z * e * B / m_i # Ion gyrofrequency [rad/s]
57 r_Le = v_th_e / omega_ce # Electron Larmor radius [m]
58 r_Li = v_th_i / omega_ci # Ion Larmor radius [m]
59
60 # Plasma beta
61 p_thermal = n_e * k * T_e * e # Thermal pressure [Pa]
62 p_magnetic = B**2 / (2 * mu_0) # Magnetic pressure [Pa]
63 beta = p_thermal / p_magnetic if p_magnetic > 0 else np.inf
64 else:
65 omega_ce = omega_ci = 0.0
66 r_Le = r_Li = np.inf
67 beta = np.inf
68
69 # Coulomb logarithm (using simplified formula)
70 if T_e > 10: # eV
71 ln_Lambda = 24 - np.log(np.sqrt(n_e * 1e-6) / T_e)
72 else:
73 ln_Lambda = 23 - np.log(np.sqrt(n_e * 1e-6) * Z / T_e**1.5)
74 ln_Lambda = max(ln_Lambda, 2.0) # Lower bound
75
76 # Collision frequency (Spitzer, electron-ion) [Hz]
77 nu_ei = (n_e * e**4 * ln_Lambda) / (
78 12 * np.pi**1.5 * epsilon_0**2 * m_e**0.5 * (T_e_J)**1.5
79 ) if T_e > 0 else np.inf
80
81 # Mean free path [m]
82 mfp = v_th_e / nu_ei if nu_ei > 0 else np.inf
83
84 return {
85 'lambda_D': lambda_D,
86 'omega_pe': omega_pe,
87 'omega_pi': omega_pi,
88 'omega_ce': omega_ce,
89 'omega_ci': omega_ci,
90 'v_th_e': v_th_e,
91 'v_th_i': v_th_i,
92 'r_Le': r_Le,
93 'r_Li': r_Li,
94 'beta': beta,
95 'ln_Lambda': ln_Lambda,
96 'nu_ei': nu_ei,
97 'mfp': mfp
98 }
99
100
101def print_parameters_table():
102 """Print a comparison table of parameters for different plasma regimes."""
103
104 # Define plasma regimes
105 plasmas = [
106 ('Tokamak Core', 1e20, 10e3, 5.0),
107 ('Solar Wind', 5e6, 10, 5e-9),
108 ('Ionosphere', 1e11, 0.1, 50e-6),
109 ('Lightning', 1e24, 3, 0.0),
110 ('Neon Sign', 1e16, 2, 0.0)
111 ]
112
113 print("=" * 100)
114 print(f"{'Plasma Type':<15} {'n_e[m^-3]':<12} {'T_e[eV]':<10} {'λ_D[m]':<12} "
115 f"{'f_pe[Hz]':<12} {'r_L[m]':<12} {'β':<10}")
116 print("=" * 100)
117
118 for name, n_e, T_e, B in plasmas:
119 params = compute_plasma_parameters(n_e, T_e, B)
120 f_pe = params['omega_pe'] / (2 * np.pi)
121 r_L = params['r_Le']
122
123 print(f"{name:<15} {n_e:<12.2e} {T_e:<10.2e} {params['lambda_D']:<12.2e} "
124 f"{f_pe:<12.2e} {r_L:<12.2e} {params['beta']:<10.2e}")
125
126 print("=" * 100)
127
128
129def plot_parameter_space():
130 """Plot plasma parameter space diagram showing different regimes."""
131
132 # Create density and temperature arrays
133 n_e = np.logspace(4, 28, 100) # m^-3
134 T_e = np.logspace(-2, 5, 100) # eV
135
136 N, T = np.meshgrid(n_e, T_e)
137
138 # Compute Debye length
139 T_J = T * e
140 lambda_D = np.sqrt(epsilon_0 * T_J / (N * e**2))
141
142 # Compute plasma parameter
143 n_D = N * (4/3) * np.pi * lambda_D**3
144
145 # Create figure
146 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
147
148 # Plot 1: Debye length contours
149 cs1 = ax1.contourf(np.log10(N), np.log10(T), np.log10(lambda_D),
150 levels=20, cmap='viridis')
151 ax1.contour(np.log10(N), np.log10(T), np.log10(lambda_D),
152 levels=[-9, -6, -3, 0, 3], colors='white', linewidths=1.5)
153
154 # Mark different plasma regimes
155 regimes = [
156 ('Tokamak\nCore', 20, 4, 'red'),
157 ('Solar\nWind', 6.7, 1, 'yellow'),
158 ('Ionosphere', 11, -1, 'cyan'),
159 ('Lightning', 24, 0.5, 'orange'),
160 ('Neon\nSign', 16, 0.3, 'magenta')
161 ]
162
163 for name, log_n, log_T, color in regimes:
164 ax1.plot(log_n, log_T, 'o', color=color, markersize=10,
165 markeredgecolor='white', markeredgewidth=2)
166 ax1.text(log_n, log_T + 0.5, name, color=color, fontsize=9,
167 ha='center', fontweight='bold')
168
169 ax1.set_xlabel('log₁₀(n_e [m⁻³])', fontsize=12)
170 ax1.set_ylabel('log₁₀(T_e [eV])', fontsize=12)
171 ax1.set_title('Debye Length in Parameter Space', fontsize=13, fontweight='bold')
172 ax1.grid(True, alpha=0.3)
173
174 cbar1 = plt.colorbar(cs1, ax=ax1)
175 cbar1.set_label('log₁₀(λ_D [m])', fontsize=11)
176
177 # Plot 2: Plasma parameter contours
178 cs2 = ax2.contourf(np.log10(N), np.log10(T), np.log10(n_D),
179 levels=20, cmap='plasma')
180 ax2.contour(np.log10(N), np.log10(T), np.log10(n_D),
181 levels=[0, 3, 6, 9], colors='white', linewidths=1.5)
182
183 # Mark same regimes
184 for name, log_n, log_T, color in regimes:
185 ax2.plot(log_n, log_T, 'o', color=color, markersize=10,
186 markeredgecolor='white', markeredgewidth=2)
187 ax2.text(log_n, log_T + 0.5, name, color=color, fontsize=9,
188 ha='center', fontweight='bold')
189
190 # Add diagonal line for n_D = 1 (plasma validity boundary)
191 n_line = np.logspace(4, 28, 50)
192 T_line = (n_line * (4/3) * np.pi * e**6) / (epsilon_0**3 * k**3)
193 valid_idx = (T_line > 1e-2) & (T_line < 1e5)
194 ax2.plot(np.log10(n_line[valid_idx]), np.log10(T_line[valid_idx] / e),
195 'r--', linewidth=2, label='n_D = 1 (boundary)')
196
197 ax2.set_xlabel('log₁₀(n_e [m⁻³])', fontsize=12)
198 ax2.set_ylabel('log₁₀(T_e [eV])', fontsize=12)
199 ax2.set_title('Plasma Parameter (n_D) Space', fontsize=13, fontweight='bold')
200 ax2.grid(True, alpha=0.3)
201 ax2.legend(loc='lower right', fontsize=10)
202
203 cbar2 = plt.colorbar(cs2, ax=ax2)
204 cbar2.set_label('log₁₀(n_D)', fontsize=11)
205
206 plt.tight_layout()
207 plt.savefig('plasma_parameter_space.png', dpi=300, bbox_inches='tight')
208 plt.show()
209
210
211def plot_frequency_comparison():
212 """Compare characteristic frequencies for tokamak plasma."""
213
214 n_e = 1e20 # m^-3
215 T_range = np.logspace(0, 4, 100) # 1 eV to 10 keV
216 B = 5.0 # Tesla
217
218 frequencies = {
219 'ω_pe': [],
220 'ω_ce': [],
221 'ω_ci': [],
222 'ν_ei': []
223 }
224
225 for T_e in T_range:
226 params = compute_plasma_parameters(n_e, T_e, B)
227 frequencies['ω_pe'].append(params['omega_pe'])
228 frequencies['ω_ce'].append(params['omega_ce'])
229 frequencies['ω_ci'].append(params['omega_ci'])
230 frequencies['ν_ei'].append(params['nu_ei'])
231
232 # Plot
233 fig, ax = plt.subplots(figsize=(10, 7))
234
235 ax.loglog(T_range, np.array(frequencies['ω_pe']) / (2*np.pi),
236 'b-', linewidth=2, label='f_pe (plasma frequency)')
237 ax.loglog(T_range, np.array(frequencies['ω_ce']) / (2*np.pi),
238 'r-', linewidth=2, label='f_ce (electron gyro)')
239 ax.loglog(T_range, np.array(frequencies['ω_ci']) / (2*np.pi),
240 'g-', linewidth=2, label='f_ci (ion gyro)')
241 ax.loglog(T_range, frequencies['ν_ei'],
242 'm--', linewidth=2, label='ν_ei (collision freq)')
243
244 ax.set_xlabel('Electron Temperature [eV]', fontsize=13)
245 ax.set_ylabel('Frequency [Hz]', fontsize=13)
246 ax.set_title(f'Characteristic Frequencies\n(n_e = {n_e:.1e} m⁻³, B = {B} T)',
247 fontsize=14, fontweight='bold')
248 ax.grid(True, alpha=0.3, which='both')
249 ax.legend(loc='best', fontsize=11)
250
251 # Add regions
252 ax.axhspan(1e6, 1e9, alpha=0.1, color='yellow', label='Radio waves')
253 ax.axhspan(1e9, 1e12, alpha=0.1, color='cyan', label='Microwaves')
254
255 plt.tight_layout()
256 plt.savefig('plasma_frequencies.png', dpi=300, bbox_inches='tight')
257 plt.show()
258
259
260if __name__ == '__main__':
261 print("\n" + "="*100)
262 print("PLASMA PARAMETERS CALCULATOR")
263 print("="*100 + "\n")
264
265 # Print comparison table
266 print_parameters_table()
267
268 # Example calculation for tokamak
269 print("\n\nDetailed calculation for Tokamak Core:")
270 print("-" * 60)
271 n_e = 1e20 # m^-3
272 T_e = 10e3 # eV
273 B = 5.0 # Tesla
274
275 params = compute_plasma_parameters(n_e, T_e, B)
276
277 print(f"Density: {n_e:.2e} m^-3")
278 print(f"Temperature: {T_e:.2e} eV")
279 print(f"Magnetic field: {B} T\n")
280
281 print(f"Debye length: {params['lambda_D']:.4e} m")
282 print(f"Plasma frequency (e): {params['omega_pe']/(2*np.pi):.4e} Hz")
283 print(f"Plasma frequency (i): {params['omega_pi']/(2*np.pi):.4e} Hz")
284 print(f"Gyrofrequency (e): {params['omega_ce']/(2*np.pi):.4e} Hz")
285 print(f"Gyrofrequency (i): {params['omega_ci']/(2*np.pi):.4e} Hz")
286 print(f"Thermal velocity (e): {params['v_th_e']:.4e} m/s")
287 print(f"Thermal velocity (i): {params['v_th_i']:.4e} m/s")
288 print(f"Larmor radius (e): {params['r_Le']:.4e} m")
289 print(f"Larmor radius (i): {params['r_Li']:.4e} m")
290 print(f"Plasma beta: {params['beta']:.4e}")
291 print(f"Coulomb logarithm: {params['ln_Lambda']:.2f}")
292 print(f"Collision frequency: {params['nu_ei']:.4e} Hz")
293 print(f"Mean free path: {params['mfp']:.4e} m")
294
295 # Generate plots
296 print("\n\nGenerating parameter space diagrams...")
297 plot_parameter_space()
298
299 print("Generating frequency comparison plot...")
300 plot_frequency_comparison()
301
302 print("\nDone! Check plasma_parameter_space.png and plasma_frequencies.png")