This page was generated by nbsphinx from docs/notebooks/cold_plasma_tensor_elements.ipynb.
Interactive online version: .

%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 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

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

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)