Repository: PlasmaPy/PlasmaPy
Pull request #224 (2018-01-04 13:54:38)

Magnetized Cold Plasma tensor elements functions
Added functions to calculate Magnetized Cold Plasma tensor elements (aka "Stix" tensor), which is useful for waves in plasma physics. Also added an example which plot a graph for a wide range of frequencies.



Patch:
diff --git a/plasmapy/examples/plot_cold_plasma_tensor_elements.py b/plasmapy/examples/plot_cold_plasma_tensor_elements.py
new file mode 100644
index 00000000..5b09b61f
--- /dev/null
+++ b/plasmapy/examples/plot_cold_plasma_tensor_elements.py
@@ -0,0 +1,110 @@
+# coding: utf-8
+"""
+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)
diff --git a/plasmapy/physics/__init__.py b/plasmapy/physics/__init__.py
index ebe38e3d..c0024cda 100644
--- a/plasmapy/physics/__init__.py
+++ b/plasmapy/physics/__init__.py
@@ -6,3 +6,4 @@
 from .relativity import *
 from .transport import *
 from .distribution import *
+from .dielectric import *
\ No newline at end of file
diff --git a/plasmapy/physics/dielectric.py b/plasmapy/physics/dielectric.py
new file mode 100644
index 00000000..b634c07b
--- /dev/null
+++ b/plasmapy/physics/dielectric.py
@@ -0,0 +1,197 @@
+"""Functions to calculate plasma dielectric paramaters"""
+
+from astropy import units
+import plasmapy.utils as utils
+
+from ..constants import (pi, m_e, c, mu0, e, eps0)
+from .parameters import gyrofrequency, plasma_frequency
+
+r"""
+Values should be returned as an Astropy Quantity in SI units.
+"""
+
+
+@utils.check_quantity({
+    'B': {'units': units.T, 'can_be_negative': False},
+    'omega': {'units': units.rad/units.s, 'can_be_negative': False},
+})
+def cold_plasma_permittivity_SDP(B, species, n, omega):
+    r"""
+    Magnetized Cold Plasma Dielectric Permittivity Tensor Elements.
+    Elements (S, D, P) are given in the "Stix" frame, ie. with B // z.
+
+    The :math:`\exp(-i \omega t)` time-harmonic convention is assumed.
+
+
+    Parameters
+    ----------
+    B : ~astropy.units.Quantity
+        Magnetic field magnitude in units convertible to tesla.
+
+    species : List of string
+        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
+    -------
+    S : ~astropy.units.Quantity
+        S ("Sum") dielectric tensor element
+
+    D : ~astropy.units.Quantity
+        D ("Difference") dielectric tensor element
+
+    P : ~astropy.units.Quantity
+        P ("Plasma") dielectric tensor element
+
+    Notes
+    -----
+    The dielectric permittivity tensor is expressed in the Stix frame with
+    the :math:`\exp(-i \omega t)` time-harmonic convention as
+    :math:`\varepsilon = \varepsilon_0 A`, with :math:`A` being
+
+    .. math::
+
+        \varepsilon = \varepsilon_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{\omega_{p,s}^2}{\omega^2 - \Omega_{c,s}^2}
+
+        D = \sum_s \frac{\Omega_{c,s}}{\omega}
+            \frac{\omega_{p,s}^2}{\omega^2 - \Omega_{c,s}^2}
+
+        P = 1 - \sum_s \frac{\omega_{p,s}^2}{\omega^2}
+
+    where :math:`\omega_{p,s}` is the plasma frequency and
+    :math:`\Omega_{c,s}` is the signed version of the cyclotron frequency
+    for the species :math:`s`.
+
+    References
+    ----------
+    - T.H. Stix, Waves in Plasma, 1992.
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from plasmapy.constants 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)
+    >>> S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
+    >>> S
+    <Quantity 1.0242290172193407>
+    >>> D
+    <Quantity 0.3908935245261273>
+    >>> P
+    <Quantity -4.890310398611263>
+    """
+    S, D, P = 1, 0, 1
+
+    for s, n_s in zip(species, n):
+        omega_c = gyrofrequency(B=B, particle=s, signed=True)
+        omega_p = plasma_frequency(n=n_s, particle=s)
+
+        S += - omega_p ** 2 / (omega ** 2 - omega_c ** 2)
+        D += omega_c/omega * omega_p ** 2 / (omega ** 2 - omega_c ** 2)
+        P += - omega_p ** 2 / omega ** 2
+    return S, D, P
+
+
+def cold_plasma_permittivity_LRP(B, species, n, omega):
+    r"""
+    Magnetized Cold Plasma Dielectric Permittivity Tensor Elements.
+    Elements (L, R, P) are given in the "rotating" basis, ie. in the basis
+    :math:`(\mathbf{u}_{+}, \mathbf{u}_{-}, \mathbf{u}_z)`,
+    where the tensor is diagonal and with B // z.
+
+    The :math:`\exp(-i \omega t)` time-harmonic convention is assumed.
+
+
+    Parameters
+    ----------
+    B : ~astropy.units.Quantity
+        Magnetic field magnitude in units convertible to tesla.
+
+    species : List of string
+        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
+    -------
+    L : ~astropy.units.Quantity
+        L ("Left") Left-handed circularly polarization tensor element
+
+    R : ~astropy.units.Quantity
+        R ("Right") Right-handed circularly polarization tensor element
+
+    P : ~astropy.units.Quantity
+        P ("Plasma") dielectric tensor element
+
+    Notes
+    -----
+    In the rotating frame defined by
+     :math:`(\mathbf{u}_{+}, \mathbf{u}_{-}, \mathbf{u}_z)`
+    with :math:`\mathbf{u}_{\pm}=(\mathbf{u}_x \pm \mathbf{u}_y)/\sqrt{2}`,
+     the dielectric tensor takes a diagonal form with elements L, R, P with:
+
+    .. math::
+        L = 1 - \sum_s
+                \frac{\omega_{p,s}^2}{\omega\left(\omega - \Omega_{c,s}\right)}
+
+        R = 1 - \sum_s
+                \frac{\omega_{p,s}^2}{\omega\left(\omega + \Omega_{c,s}\right)}
+
+        P = 1 - \sum_s \frac{\omega_{p,s}^2}{\omega^2}
+
+    where :math:`\omega_{p,s}` is the plasma frequency and
+    :math:`\Omega_{c,s}` is the signed version of the cyclotron frequency
+    for the species :math:`s`.
+
+    References
+    ----------
+    - T.H. Stix, Waves in Plasma, 1992.
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from plasmapy.constants 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)
+    >>> L, R, P = cold_plasma_permittivity_LRP(B, species, n, omega)
+    >>> L
+    <Quantity 0.6333354926932133>
+    >>> R
+    <Quantity 1.4151225417454678>
+    >>> P
+    <Quantity -4.890310398611263>
+    """
+    L, R, P = 1, 1, 1
+
+    for s, n_s in zip(species, n):
+        omega_c = gyrofrequency(B=B, particle=s, signed=True)
+        omega_p = plasma_frequency(n=n_s, particle=s)
+
+        L += - omega_p ** 2 / (omega*(omega - omega_c))
+        R += - omega_p ** 2 / (omega*(omega + omega_c))
+        P += - omega_p ** 2 / omega ** 2
+    return L, R, P


Test patch:
diff --git a/plasmapy/physics/tests/test_dielectric.py b/plasmapy/physics/tests/test_dielectric.py
new file mode 100644
index 00000000..f1e0f512
--- /dev/null
+++ b/plasmapy/physics/tests/test_dielectric.py
@@ -0,0 +1,77 @@
+"""Tests for functions that calculate plasma dielectric parameters in
+dielectry.py"""
+
+import numpy as np
+from astropy import units as u
+
+from ..dielectric import (cold_plasma_permittivity_LRP,
+                          cold_plasma_permittivity_SDP)
+
+from ..parameters import (plasma_frequency, gyrofrequency)
+
+B = 1.0 * u.T
+n = [1e18/u.m**3]
+omega = 55e6*u.rad/u.s
+
+single_species = ['e']
+two_species = ['e', 'D+']
+three_species = ['e', 'D+', 'H+']
+
+
+class Test_ColdPlasmaPermittivity(object):
+    def test_proton_electron_plasma(self):
+        """
+        Test proton-electron plasma against the (approximate) 
+        analytical formulas
+        """
+        B = 1*u.T
+        n = [1, 1]*1/u.m**3
+        omega = 1*u.rad/u.s
+        omega_ce = gyrofrequency(B, particle='e', signed=True)
+        omega_pe = plasma_frequency(n[0], particle='e')
+        omega_cp = abs(omega_ce)/1860
+        omega_pp = omega_pe/43
+
+        S_analytical = 1 \
+            - omega_pe ** 2 / (omega ** 2 - omega_ce ** 2) \
+            - omega_pp ** 2 / (omega ** 2 - omega_cp ** 2)
+
+        D_analytical = \
+            + omega_ce / omega * omega_pe ** 2 / (omega ** 2 - omega_ce ** 2) \
+            + omega_cp / omega * omega_pp ** 2 / (omega ** 2 - omega_cp ** 2)
+
+        P_analytical = 1 - (omega_pe ** 2 + omega_pp ** 2) / omega ** 2
+
+        species = ['e', 'p']
+        S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
+
+        assert np.isclose(S, S_analytical)
+        assert np.isclose(D, D_analytical)
+        assert np.isclose(P, P_analytical)
+
+    def test_three_species(self):
+        """
+        Test with three species (2 ions): D plasma with 5%H minority fraction
+        """
+        n_3 = np.array([1, 1, 5/100])*1e19/u.m**3
+        S, D, P = cold_plasma_permittivity_SDP(B, three_species, n_3, omega)
+        assert np.isclose(S, -11753.3)
+        assert np.isclose(D, 13408.99181054283)
+        assert np.isclose(P, -10524167.9)
+
+    def test_SD_to_LR_relationships(self):
+        """
+        Test the relationships between (S, D, P) notation in Stix basis and
+        (L, R, P) notation in the rotating basis, ie test:
+         S = (R+L)/2 and D = (R-L)/2
+        and
+         R = S+D and L = S-D
+        """
+        # test with a single species
+        S, D, _ = cold_plasma_permittivity_SDP(B, single_species, n, omega)
+        L, R, _ = cold_plasma_permittivity_LRP(B, single_species, n, omega)
+
+        assert np.isclose(R, S+D)
+        assert np.isclose(L, S-D)
+        assert np.isclose(S, (R+L)/2)
+        assert np.isclose(D, (R-L)/2)