/
plot_cold_plasma_tensor_elements.py
109 lines (98 loc) · 3.64 KB
/
plot_cold_plasma_tensor_elements.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
"""
Cold Magnetized Plasma Waves Tensor Elements (S, D, P in Stix's notation)
=========================================================================
This example shows how to calculate the values of the cold plasma tensor
elements for various electromagnetic wave frequencies.
"""
# First, import some basics (and `PlasmaPy`!)
import numpy as np
import plasmapy as pp
import matplotlib.pyplot as plt
from astropy import units as u
#######################################################################
# Let's define some parameters, such as the magnetic field magnitude,
# the plasma species and densities and the frequency band of interest
B = 2 * u.T
species = ['e', 'D+']
n = [1e18 * u.m ** -3, 1e18 * u.m ** -3]
f = np.logspace(start=6, stop=11.3, num=3001) # 1 MHz to 200 GHz
omega_RF = f * (2 * np.pi) * (u.rad / u.s)
#######################################################################
help(pp.physics.cold_plasma_permittivity_SDP)
#######################################################################
S, D, P = pp.physics.cold_plasma_permittivity_SDP(B, species, n, omega_RF)
#######################################################################
# Filter positive and negative values, for display purposes only.
# Still for display purposes, replace 0 by NaN to NOT plot 0 values
S_pos = S * (S > 0)
D_pos = D * (D > 0)
P_pos = P * (P > 0)
S_neg = S * (S < 0)
D_neg = D * (D < 0)
P_neg = P * (P < 0)
S_pos[S_pos == 0] = np.NaN
D_pos[D_pos == 0] = np.NaN
P_pos[P_pos == 0] = np.NaN
S_neg[S_neg == 0] = np.NaN
D_neg[D_neg == 0] = np.NaN
P_neg[P_neg == 0] = np.NaN
#######################################################################
plt.figure(figsize=(12, 6))
plt.semilogx(f, abs(S_pos),
f, abs(D_pos),
f, abs(P_pos), lw=2)
plt.semilogx(f, abs(S_neg), '#1f77b4',
f, abs(D_neg), '#ff7f0e',
f, abs(P_neg), '#2ca02c', lw=2, ls='--')
plt.yscale('log')
plt.grid(True, which='major')
plt.grid(True, which='minor')
plt.ylim(1e-4, 1e8)
plt.xlim(1e6, 200e9)
plt.legend(('S > 0', 'D > 0', 'P > 0', 'S < 0', 'D < 0', 'P < 0'),
fontsize=16, ncol=2)
plt.xlabel('RF Frequency [Hz]', size=16)
plt.ylabel('Absolute value', size=16)
plt.tick_params(labelsize=14)
#######################################################################
# Cold Plasma tensor elements in the rotating basis
L, R, P = pp.physics.cold_plasma_permittivity_LRP(B, species, n, omega_RF)
#######################################################################
L_pos = L * (L > 0)
R_pos = R * (R > 0)
L_neg = L * (L < 0)
R_neg = R * (R < 0)
L_pos[L_pos == 0] = np.NaN
R_pos[R_pos == 0] = np.NaN
L_neg[L_neg == 0] = np.NaN
R_neg[R_neg == 0] = np.NaN
plt.figure(figsize=(12, 6))
plt.semilogx(f, abs(L_pos),
f, abs(R_pos),
f, abs(P_pos), lw=2)
plt.semilogx(f, abs(L_neg), '#1f77b4',
f, abs(R_neg), '#ff7f0e',
f, abs(P_neg), '#2ca02c', lw=2, ls='--')
plt.yscale('log')
plt.grid(True, which='major')
plt.grid(True, which='minor')
plt.xlim(1e6, 200e9)
plt.legend(('L > 0', 'R > 0', 'P > 0', 'L < 0', 'R < 0', 'P < 0'),
fontsize=16, ncol=2)
plt.xlabel('RF Frequency [Hz]', size=16)
plt.ylabel('Absolute value', size=16)
plt.tick_params(labelsize=14)
#######################################################################
# Checks if the values obtained are coherent. They should satisfy
# S = (R+L)/2 and D = (R-L)/2
try:
np.testing.assert_allclose(S, (R + L) / 2)
np.testing.assert_allclose(D, (R - L) / 2)
except AssertionError as e:
print(e)
# Checks for R=S+D and L=S-D
try:
np.testing.assert_allclose(R, S + D)
np.testing.assert_allclose(L, S - D)
except AssertionError as e:
print(e)