This page was generated by
nbsphinx from
docs/notebooks/cold_plasma_tensor_elements.ipynb.
Interactive online version:
.
[1]:
%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.
[2]:
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 cold_plasma_permittivity_SDP()
and cold_plasma_permittivity_LRP()
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
[3]:
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)
[4]:
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`
Electromagnetic wave frequency in rad/s.
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]
>>> omega = 3.7e9*(2*pi)*(u.rad/u.s)
>>> 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...>
[5]:
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
[6]:
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
[7]:
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
[8]:
L, R, P = cold_plasma_permittivity_LRP(B, species, n, omega_RF)
[9]:
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
[10]:
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)