1#!/usr/bin/env python3
2"""
3Velocity Distribution Functions
4
5This script plots and compares various velocity distribution functions used
6in plasma physics: Maxwellian, bi-Maxwellian, kappa, shifted Maxwellian,
7and bump-on-tail distributions.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import k, m_p, m_e
16from scipy.special import gamma
17
18
19def maxwellian_1d(v, n, m, T):
20 """
21 1D Maxwellian distribution.
22
23 f(v) = n * sqrt(m / (2πkT)) * exp(-m*v² / (2kT))
24
25 Parameters:
26 -----------
27 v : array
28 Velocity [m/s]
29 n : float
30 Density [m^-3]
31 m : float
32 Mass [kg]
33 T : float
34 Temperature [K]
35
36 Returns:
37 --------
38 f : array
39 Distribution function [s/m^4]
40 """
41 v_th = np.sqrt(2 * k * T / m)
42 return n * np.sqrt(m / (2 * np.pi * k * T)) * np.exp(-v**2 / v_th**2)
43
44
45def maxwellian_3d(v, n, m, T):
46 """
47 3D Maxwellian speed distribution (spherical).
48
49 f(v) = n * (m / (2πkT))^(3/2) * exp(-m*v² / (2kT))
50
51 Speed distribution: 4πv² * f(v)
52
53 Parameters:
54 -----------
55 v : array
56 Speed [m/s]
57 n, m, T : float
58 Density, mass, temperature
59
60 Returns:
61 --------
62 f_speed : array
63 Speed distribution [s/m^4]
64 """
65 v_th = np.sqrt(2 * k * T / m)
66 f = n * (m / (2 * np.pi * k * T))**(3/2) * np.exp(-v**2 / v_th**2)
67 return 4 * np.pi * v**2 * f
68
69
70def bi_maxwellian(v_perp, v_para, n, m, T_perp, T_para):
71 """
72 Bi-Maxwellian distribution with different perpendicular and parallel temps.
73
74 f(v_perp, v_para) = n * (m / (2π))^(3/2) / (T_perp * sqrt(T_para))
75 * exp(-m*v_perp² / (2*T_perp) - m*v_para² / (2*T_para))
76
77 Parameters:
78 -----------
79 v_perp, v_para : array
80 Perpendicular and parallel velocities [m/s]
81 n, m : float
82 Density and mass
83 T_perp, T_para : float
84 Perpendicular and parallel temperatures [K]
85
86 Returns:
87 --------
88 f : array
89 Distribution function [s/m^4]
90 """
91 v_th_perp = np.sqrt(2 * k * T_perp / m)
92 v_th_para = np.sqrt(2 * k * T_para / m)
93
94 return (n * (m / (2 * np.pi))**(3/2) /
95 (T_perp * np.sqrt(T_para)) *
96 np.exp(-v_perp**2 / v_th_perp**2 - v_para**2 / v_th_para**2))
97
98
99def kappa_distribution_1d(v, n, m, T, kappa):
100 """
101 Kappa distribution (superthermal tails).
102
103 f(v) ∝ [1 + v² / (κ * v_th²)]^(-(κ+1))
104
105 Parameters:
106 -----------
107 v : array
108 Velocity [m/s]
109 n, m, T : float
110 Density, mass, temperature
111 kappa : float
112 Kappa parameter (κ → ∞ recovers Maxwellian)
113
114 Returns:
115 --------
116 f : array
117 Distribution function [s/m^4]
118 """
119 v_th = np.sqrt(2 * k * T / m)
120
121 # Normalization factor
122 norm = (n * gamma(kappa + 1) /
123 (np.sqrt(np.pi * kappa) * gamma(kappa - 0.5) * v_th))
124
125 return norm * (1 + v**2 / (kappa * v_th**2))**(-kappa - 1)
126
127
128def shifted_maxwellian(v, n, m, T, v_drift):
129 """
130 Shifted Maxwellian (beam).
131
132 f(v) = n * sqrt(m / (2πkT)) * exp(-m*(v - v_drift)² / (2kT))
133
134 Parameters:
135 -----------
136 v : array
137 Velocity [m/s]
138 n, m, T : float
139 Density, mass, temperature
140 v_drift : float
141 Drift velocity [m/s]
142
143 Returns:
144 --------
145 f : array
146 Distribution function [s/m^4]
147 """
148 v_th = np.sqrt(2 * k * T / m)
149 return n * np.sqrt(m / (2 * np.pi * k * T)) * np.exp(-(v - v_drift)**2 / v_th**2)
150
151
152def bump_on_tail(v, n_bg, n_beam, m, T_bg, T_beam, v_beam):
153 """
154 Bump-on-tail: background Maxwellian + beam.
155
156 f = f_background + f_beam
157
158 Parameters:
159 -----------
160 v : array
161 Velocity [m/s]
162 n_bg, n_beam : float
163 Background and beam densities
164 m : float
165 Mass
166 T_bg, T_beam : float
167 Background and beam temperatures
168 v_beam : float
169 Beam velocity
170
171 Returns:
172 --------
173 f : array
174 Distribution function [s/m^4]
175 """
176 f_bg = maxwellian_1d(v, n_bg, m, T_bg)
177 f_beam = shifted_maxwellian(v, n_beam, m, T_beam, v_beam)
178 return f_bg + f_beam
179
180
181def plot_1d_distributions():
182 """Plot 1D velocity distributions."""
183
184 # Parameters
185 n = 1e20 # m^-3
186 m = m_e
187 T = 1e4 # K (~ 1 eV)
188 v_th = np.sqrt(2 * k * T / m)
189
190 # Velocity array
191 v = np.linspace(-5*v_th, 5*v_th, 500)
192
193 # Different distributions
194 f_max = maxwellian_1d(v, n, m, T)
195 f_kappa_2 = kappa_distribution_1d(v, n, m, T, kappa=2)
196 f_kappa_5 = kappa_distribution_1d(v, n, m, T, kappa=5)
197 f_kappa_10 = kappa_distribution_1d(v, n, m, T, kappa=10)
198 f_shifted = shifted_maxwellian(v, n, m, T, v_drift=2*v_th)
199
200 # Create figure
201 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
202
203 # Linear scale
204 ax1.plot(v / v_th, f_max * v_th, 'b-', linewidth=2.5, label='Maxwellian')
205 ax1.plot(v / v_th, f_kappa_2 * v_th, 'r-', linewidth=2, label='κ = 2')
206 ax1.plot(v / v_th, f_kappa_5 * v_th, 'g-', linewidth=2, label='κ = 5')
207 ax1.plot(v / v_th, f_kappa_10 * v_th, 'm-', linewidth=2, label='κ = 10')
208
209 ax1.set_xlabel('v / v_th', fontsize=12)
210 ax1.set_ylabel('f(v) × v_th [normalized]', fontsize=12)
211 ax1.set_title('Maxwellian vs Kappa Distributions', fontsize=13, fontweight='bold')
212 ax1.legend(loc='best', fontsize=10)
213 ax1.grid(True, alpha=0.3)
214 ax1.set_xlim([-5, 5])
215
216 # Log scale (shows superthermal tails)
217 ax2.semilogy(v / v_th, f_max * v_th, 'b-', linewidth=2.5, label='Maxwellian')
218 ax2.semilogy(v / v_th, f_kappa_2 * v_th, 'r-', linewidth=2, label='κ = 2')
219 ax2.semilogy(v / v_th, f_kappa_5 * v_th, 'g-', linewidth=2, label='κ = 5')
220 ax2.semilogy(v / v_th, f_kappa_10 * v_th, 'm-', linewidth=2, label='κ = 10')
221
222 ax2.set_xlabel('v / v_th', fontsize=12)
223 ax2.set_ylabel('f(v) × v_th [log scale]', fontsize=12)
224 ax2.set_title('Superthermal Tails (Log Scale)', fontsize=13, fontweight='bold')
225 ax2.legend(loc='best', fontsize=10)
226 ax2.grid(True, alpha=0.3)
227 ax2.set_xlim([-5, 5])
228 ax2.set_ylim([1e-6, 1])
229
230 plt.tight_layout()
231 plt.savefig('distribution_1d_kappa.png', dpi=300, bbox_inches='tight')
232 plt.show()
233
234
235def plot_bump_on_tail():
236 """Plot bump-on-tail distribution."""
237
238 # Parameters
239 n_bg = 1e20
240 n_beam = 0.1 * n_bg # 10% beam
241 m = m_e
242 T_bg = 1e4 # K
243 T_beam = 0.5 * T_bg # Beam is colder
244 v_beam = 3 * np.sqrt(2 * k * T_bg / m) # Beam at 3 v_th
245
246 v_th_bg = np.sqrt(2 * k * T_bg / m)
247
248 # Velocity array
249 v = np.linspace(-2*v_th_bg, 6*v_th_bg, 1000)
250
251 # Distributions
252 f_bg = maxwellian_1d(v, n_bg, m, T_bg)
253 f_beam = shifted_maxwellian(v, n_beam, m, T_beam, v_beam)
254 f_total = bump_on_tail(v, n_bg, n_beam, m, T_bg, T_beam, v_beam)
255
256 # df/dv (slope)
257 dfdv = np.gradient(f_total, v)
258
259 # Create figure
260 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
261
262 # Plot 1: Distribution
263 ax1.plot(v / v_th_bg, f_bg * v_th_bg, 'b--', linewidth=2,
264 label='Background', alpha=0.7)
265 ax1.plot(v / v_th_bg, f_beam * v_th_bg, 'r--', linewidth=2,
266 label='Beam', alpha=0.7)
267 ax1.plot(v / v_th_bg, f_total * v_th_bg, 'k-', linewidth=2.5,
268 label='Total (bump-on-tail)')
269
270 ax1.fill_between(v / v_th_bg, 0, f_total * v_th_bg, alpha=0.2, color='green',
271 label='Unstable region')
272
273 ax1.set_xlabel('v / v_th', fontsize=12)
274 ax1.set_ylabel('f(v) × v_th', fontsize=12)
275 ax1.set_title('Bump-on-Tail Distribution', fontsize=13, fontweight='bold')
276 ax1.legend(loc='best', fontsize=10)
277 ax1.grid(True, alpha=0.3)
278
279 # Plot 2: Slope df/dv
280 ax2.plot(v / v_th_bg, dfdv * v_th_bg**2, 'b-', linewidth=2.5)
281 ax2.axhline(0, color='red', linestyle='--', linewidth=2, label='df/dv = 0')
282
283 # Highlight positive slope region (unstable)
284 positive_slope = dfdv > 0
285 ax2.fill_between(v / v_th_bg, 0, dfdv * v_th_bg**2,
286 where=positive_slope, alpha=0.3, color='red',
287 label='Positive slope (unstable)')
288
289 ax2.set_xlabel('v / v_th', fontsize=12)
290 ax2.set_ylabel('df/dv × v_th²', fontsize=12)
291 ax2.set_title('Slope: Positive Slope → Instability', fontsize=13, fontweight='bold')
292 ax2.legend(loc='best', fontsize=10)
293 ax2.grid(True, alpha=0.3)
294
295 plt.tight_layout()
296 plt.savefig('bump_on_tail_distribution.png', dpi=300, bbox_inches='tight')
297 plt.show()
298
299
300def plot_2d_biMaxwellian():
301 """Plot 2D bi-Maxwellian distribution."""
302
303 # Parameters
304 n = 1e20
305 m = m_e
306 T_perp = 2e4 # K (perpendicular hotter)
307 T_para = 1e4 # K
308
309 v_th_perp = np.sqrt(2 * k * T_perp / m)
310 v_th_para = np.sqrt(2 * k * T_para / m)
311
312 # 2D velocity grid
313 v_perp = np.linspace(-4*v_th_perp, 4*v_th_perp, 200)
314 v_para = np.linspace(-4*v_th_para, 4*v_th_para, 200)
315 V_perp, V_para = np.meshgrid(v_perp, v_para)
316
317 # Bi-Maxwellian
318 f_bi = bi_maxwellian(V_perp, V_para, n, m, T_perp, T_para)
319
320 # Isotropic Maxwellian for comparison (T = average)
321 T_avg = (2 * T_perp + T_para) / 3
322 v_th_avg = np.sqrt(2 * k * T_avg / m)
323 f_iso = bi_maxwellian(V_perp, V_para, n, m, T_avg, T_avg)
324
325 # Create figure
326 fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))
327
328 # Plot 1: Bi-Maxwellian contours
329 levels = np.logspace(np.log10(f_bi.max()) - 6, np.log10(f_bi.max()), 15)
330 cs1 = ax1.contourf(V_para / v_th_para, V_perp / v_th_perp, f_bi,
331 levels=levels, cmap='hot', norm=plt.matplotlib.colors.LogNorm())
332 ax1.contour(V_para / v_th_para, V_perp / v_th_perp, f_bi,
333 levels=levels[::2], colors='white', linewidths=0.5, alpha=0.5)
334
335 ax1.set_xlabel('v_∥ / v_th,∥', fontsize=12)
336 ax1.set_ylabel('v_⊥ / v_th,⊥', fontsize=12)
337 ax1.set_title(f'Bi-Maxwellian\n(T_⊥ = {T_perp/1e4:.1f}×10⁴ K, T_∥ = {T_para/1e4:.1f}×10⁴ K)',
338 fontsize=12, fontweight='bold')
339 ax1.set_aspect('equal')
340 plt.colorbar(cs1, ax=ax1, label='f(v_⊥, v_∥)')
341
342 # Plot 2: Isotropic Maxwellian
343 cs2 = ax2.contourf(V_para / v_th_avg, V_perp / v_th_avg, f_iso,
344 levels=levels, cmap='hot', norm=plt.matplotlib.colors.LogNorm())
345 ax2.contour(V_para / v_th_avg, V_perp / v_th_avg, f_iso,
346 levels=levels[::2], colors='white', linewidths=0.5, alpha=0.5)
347
348 ax2.set_xlabel('v_∥ / v_th', fontsize=12)
349 ax2.set_ylabel('v_⊥ / v_th', fontsize=12)
350 ax2.set_title(f'Isotropic Maxwellian\n(T = {T_avg/1e4:.1f}×10⁴ K)',
351 fontsize=12, fontweight='bold')
352 ax2.set_aspect('equal')
353 plt.colorbar(cs2, ax=ax2, label='f(v_⊥, v_∥)')
354
355 # Plot 3: 1D cuts
356 idx_para_0 = np.argmin(np.abs(v_para))
357 idx_perp_0 = np.argmin(np.abs(v_perp))
358
359 f_cut_para = f_bi[:, idx_perp_0] # v_perp = 0 cut
360 f_cut_perp = f_bi[idx_para_0, :] # v_para = 0 cut
361
362 ax3.semilogy(v_para / v_th_para, f_cut_para, 'b-', linewidth=2.5,
363 label='f(v_∥, v_⊥=0)')
364 ax3.semilogy(v_perp / v_th_perp, f_cut_perp, 'r-', linewidth=2.5,
365 label='f(v_∥=0, v_⊥)')
366
367 ax3.set_xlabel('v / v_th', fontsize=12)
368 ax3.set_ylabel('f(v)', fontsize=12)
369 ax3.set_title('1D Cuts Through Distribution', fontsize=12, fontweight='bold')
370 ax3.legend(loc='best', fontsize=10)
371 ax3.grid(True, alpha=0.3)
372
373 plt.tight_layout()
374 plt.savefig('bi_maxwellian_2d.png', dpi=300, bbox_inches='tight')
375 plt.show()
376
377
378def plot_3d_speed_distribution():
379 """Plot 3D speed distributions."""
380
381 # Parameters
382 n = 1e20
383 m = m_e
384 T = 1e4 # K
385 v_th = np.sqrt(2 * k * T / m)
386
387 # Speed array (only positive)
388 v = np.linspace(0, 5*v_th, 500)
389
390 # 3D speed distributions
391 f_maxwell = maxwellian_3d(v, n, m, T)
392 f_kappa_2 = 4 * np.pi * v**2 * kappa_distribution_1d(v, n, m, T, kappa=2)
393 f_kappa_5 = 4 * np.pi * v**2 * kappa_distribution_1d(v, n, m, T, kappa=5)
394
395 # Most probable, mean, and RMS speeds for Maxwellian
396 v_mp = np.sqrt(2 * k * T / m) # Most probable
397 v_mean = np.sqrt(8 * k * T / (np.pi * m)) # Mean
398 v_rms = np.sqrt(3 * k * T / m) # RMS
399
400 # Create figure
401 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
402
403 # Linear scale
404 ax1.plot(v / v_th, f_maxwell, 'b-', linewidth=2.5, label='Maxwellian')
405 ax1.plot(v / v_th, f_kappa_2, 'r-', linewidth=2, label='κ = 2')
406 ax1.plot(v / v_th, f_kappa_5, 'g-', linewidth=2, label='κ = 5')
407
408 # Mark characteristic speeds
409 ax1.axvline(v_mp / v_th, color='cyan', linestyle='--', linewidth=2,
410 label=f'v_mp = {v_mp/v_th:.2f} v_th')
411 ax1.axvline(v_mean / v_th, color='orange', linestyle='--', linewidth=2,
412 label=f'v_mean = {v_mean/v_th:.2f} v_th')
413 ax1.axvline(v_rms / v_th, color='magenta', linestyle='--', linewidth=2,
414 label=f'v_rms = {v_rms/v_th:.2f} v_th')
415
416 ax1.set_xlabel('v / v_th', fontsize=12)
417 ax1.set_ylabel('4πv² f(v)', fontsize=12)
418 ax1.set_title('3D Speed Distribution', fontsize=13, fontweight='bold')
419 ax1.legend(loc='best', fontsize=9)
420 ax1.grid(True, alpha=0.3)
421
422 # Log scale
423 ax2.semilogy(v / v_th, f_maxwell, 'b-', linewidth=2.5, label='Maxwellian')
424 ax2.semilogy(v / v_th, f_kappa_2, 'r-', linewidth=2, label='κ = 2')
425 ax2.semilogy(v / v_th, f_kappa_5, 'g-', linewidth=2, label='κ = 5')
426
427 ax2.set_xlabel('v / v_th', fontsize=12)
428 ax2.set_ylabel('4πv² f(v) [log]', fontsize=12)
429 ax2.set_title('Speed Distribution (Log Scale)', fontsize=13, fontweight='bold')
430 ax2.legend(loc='best', fontsize=10)
431 ax2.grid(True, alpha=0.3)
432 ax2.set_ylim([1e-2, f_maxwell.max() * 2])
433
434 plt.tight_layout()
435 plt.savefig('speed_distribution_3d.png', dpi=300, bbox_inches='tight')
436 plt.show()
437
438
439if __name__ == '__main__':
440 print("\n" + "="*80)
441 print("VELOCITY DISTRIBUTION FUNCTIONS IN PLASMAS")
442 print("="*80 + "\n")
443
444 print("Generating plots...")
445 print(" 1. 1D distributions (Maxwellian vs Kappa)...")
446 plot_1d_distributions()
447
448 print(" 2. Bump-on-tail distribution...")
449 plot_bump_on_tail()
450
451 print(" 3. 2D bi-Maxwellian distribution...")
452 plot_2d_biMaxwellian()
453
454 print(" 4. 3D speed distributions...")
455 plot_3d_speed_distribution()
456
457 print("\nDone! Generated 4 figures:")
458 print(" - distribution_1d_kappa.png")
459 print(" - bump_on_tail_distribution.png")
460 print(" - bi_maxwellian_2d.png")
461 print(" - speed_distribution_3d.png")