:

%matplotlib inline


# 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.

:

import matplotlib.pyplot as plt
import numpy as np

from astropy import units as u

from plasmapy.formulary import (
cold_plasma_permittivity_LRP,
cold_plasma_permittivity_SDP,
)


Take a look at the docs to and for more information on these.

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(cold_plasma_permittivity_SDP)

Help on function cold_plasma_permittivity_SDP in module plasmapy.formulary.dielectric:

cold_plasma_permittivity_SDP(B: Unit("T"), species, n, omega: Unit("rad / s"))
Magnetized cold plasma dielectric permittivity tensor elements.

Elements (S, D, P) are given in the "Stix" frame, i.e. with
:math:B ∥ \hat{z} :cite:p:stix:1992.

The :math:\exp(-i ω t) time-harmonic convention is assumed.

Parameters
----------
B : ~astropy.units.Quantity
Magnetic field magnitude in units convertible to tesla.

species : list of str
List of the plasma particle species,
e.g.: ['e', 'D+'] or ['e', 'D+', 'He+'].

n : list of ~astropy.units.Quantity
list of species density in units convertible to per cubic meter
The order of the species densities should follow species.

omega : ~astropy.units.Quantity

Returns
-------
sum : ~astropy.units.Quantity
S ("Sum") dielectric tensor element.

difference : ~astropy.units.Quantity
D ("Difference") dielectric tensor element.

plasma : ~astropy.units.Quantity
P ("Plasma") dielectric tensor element.

Notes
-----
The dielectric permittivity tensor is expressed in the Stix frame with
the :math:\exp(-i ω t) time-harmonic convention as
:math:ε = ε_0 A, with :math:A being

.. math::

ε = ε_0 \left(\begin{matrix}  S & -i D & 0 \\
+i D & S & 0 \\
0 & 0 & P \end{matrix}\right)

where:

.. math::
S = 1 - \sum_s \frac{ω_{p,s}^2}{ω^2 - Ω_{c,s}^2}

D = \sum_s \frac{Ω_{c,s}}{ω}
\frac{ω_{p,s}^2}{ω^2 - Ω_{c,s}^2}

P = 1 - \sum_s \frac{ω_{p,s}^2}{ω^2}

where :math:ω_{p,s} is the plasma frequency and
:math:Ω_{c,s} is the signed version of the cyclotron frequency
for the species :math:s.

Examples
--------
>>> from astropy import units as u
>>> from numpy import pi
>>> B = 2*u.T
>>> species = ['e', 'D+']
>>> n = [1e18*u.m**-3, 1e18*u.m**-3]
>>> permittivity = S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
>>> S
<Quantity 1.02422...>
>>> permittivity.sum   # namedtuple-style access
<Quantity 1.02422...>
>>> D
<Quantity 0.39089...>
>>> P
<Quantity -4.8903...>


:

S, D, P = 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 = 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)