Repository: PlasmaPy/PlasmaPy
Pull request #972 (2021-01-02 07:29:04)

Added diamagnetic drift 
Addresses issue #966. User inputs the pressure gradient, background field, density, and charge. I included descriptions and examples in the docstring. Also tests (similar to ExB and force drift) were made, with an addition of the case in Fig. 3.4 of Chen. I assumed the density scale n'/n were 1 / m and verified the magnitude of the drift.

StanczakDominik: Hey, thanks for doing this! I'll take a look at it tomorrow :) Sorry for the end-of-holiday-break delay!
StanczakDominik: Hah, there we go! Quick merge, then :) Thanks!

Patch:
diff --git a/changelog/966.feature.rst b/changelog/966.feature.rst
new file mode 100644
index 00000000..dcd5f2cf
--- /dev/null
+++ b/changelog/966.feature.rst
@@ -0,0 +1,2 @@
+Added the diamagnetic drift functionality to ~plasmapy.formulary.drifts and
+tests associated with this.
diff --git a/plasmapy/formulary/drifts.py b/plasmapy/formulary/drifts.py
index f40dd094..247a62c1 100644
--- a/plasmapy/formulary/drifts.py
+++ b/plasmapy/formulary/drifts.py
@@ -2,8 +2,10 @@
 Formulas for calculating particle drifts.
 """
 __all__ = [
+    "diamagnetic_drift",
     "ExB_drift",
     "force_drift",
+    "vd_",
     "veb_",
     "vfd_",
 ]
@@ -14,6 +16,65 @@
 from plasmapy.utils.decorators import validate_quantities
 
 
+@validate_quantities
+def diamagnetic_drift(dp: u.Pa / u.m, B: u.T, n: u.m ** (-3), q: u.C) -> u.m / u.s:
+    r"""
+    Calculate the diamagnetic fluid perpendicular drift.
+
+    **Aliases:** `vd_`
+
+    Parameters
+    ----------
+    dp : ~astropy.units.Quantity
+        Pressure gradient vector
+    B  : ~astropy.units.Quantity
+        Magnetic field vector
+    n  : ~astropy.units.Quantity
+        Number density
+    q  : ~astropy.units.Quantity
+        Particle charge
+
+    Returns
+    -------
+    v: ~astropy.units.Quantity
+        Drift velocity, in m/s
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> import numpy as np
+    >>> dp = np.array([0, -1, 0]) * u.Pa / u.m
+    >>> b = np.array([0, 0, 1]) * u.T
+    >>> n = 1 / u.m ** 3
+    >>> q = 1 * u.C
+    >>> diamagnetic_drift(dp, b, n, q)
+    <Quantity [ 1., -0., -0.] m / s>
+
+    Notes
+    -----
+    The diamagnetic drift is given by
+
+    .. math::
+
+        \vec{v} = -\frac{ \nabla p \times \vec{B} }{ q n |B|^2 }
+
+    This is the velocity component of a fluid element perpendicular to the
+    magnetic field.
+
+    References
+    ----------
+    - Chen, Introduction to Plasma Physics and Controlled Fusion, 3.65
+
+    """
+
+    cross = np.cross(dp, B)
+    return -cross / q / n / (B * B).sum(-1)
+
+
+vd_ = diamagnetic_drift
+""" Alias to :func:`diamagnetic_drift`. """
+
+
 @validate_quantities
 def ExB_drift(E: u.V / u.m, B: u.T) -> u.m / u.s:
     r"""
@@ -36,15 +97,14 @@ def ExB_drift(E: u.V / u.m, B: u.T) -> u.m / u.s:
     Examples
     --------
     >>> import astropy.units as u
-    >>> from astropy.constants.si import g0, e, m_e
     >>> ex = np.array([1, 0, 0])
-    >>> ez = np.array([0, 0, 1])
-    >>> force_drift(-ez*g0*m_e, ex*0.01*u.T, e)
-    <Quantity [ 0.0000000e+00, -5.5756984e-09,  0.0000000e+00] m / s>
-    >>> force_drift(-ez*g0*m_e, ez*0.01*u.T, e)
-    <Quantity [ 0., -0.,  0.] m / s>
-    >>> force_drift(-ez*g0*m_e, ex*u.T, e)
-    <Quantity [ 0.0000000e+00, -5.5756984e-11,  0.0000000e+00] m / s>
+    >>> ey = np.array([0, 1, 0])
+    >>> ExB_drift(ex * u.V/u.m, ey * u.T)
+    <Quantity [0., 0., 1.] m / s>
+    >>> ExB_drift(ex * u.V/u.m, ex * u.T)
+    <Quantity [0., 0., 0.] m / s>
+    >>> ExB_drift(ex * u.V/u.m, 100 * ey * u.T)
+    <Quantity [0.  , 0.  , 0.01] m / s>
 
     Notes
     -----
@@ -91,14 +151,15 @@ def force_drift(F: u.N, B: u.T, q: u.C) -> u.m / u.s:
     Examples
     --------
     >>> import astropy.units as u
+    >>> from astropy.constants.si import g0, e, m_e
     >>> ex = np.array([1, 0, 0])
-    >>> ey = np.array([0, 1, 0])
-    >>> ExB_drift(ex * u.V/u.m, ey * u.T)
-    <Quantity [0., 0., 1.] m / s>
-    >>> ExB_drift(ex * u.V/u.m, ex * u.T)
-    <Quantity [0., 0., 0.] m / s>
-    >>> ExB_drift(ex * u.V/u.m, 100 * ey * u.T)
-    <Quantity [0.  , 0.  , 0.01] m / s>
+    >>> ez = np.array([0, 0, 1])
+    >>> force_drift(-ez*g0*m_e, ex*0.01*u.T, e)
+    <Quantity [ 0.0000000e+00, -5.5756984e-09,  0.0000000e+00] m / s>
+    >>> force_drift(-ez*g0*m_e, ez*0.01*u.T, e)
+    <Quantity [ 0., -0.,  0.] m / s>
+    >>> force_drift(-ez*g0*m_e, ex*u.T, e)
+    <Quantity [ 0.0000000e+00, -5.5756984e-11,  0.0000000e+00] m / s>
 
     Returns
     -------


Test patch:
diff --git a/plasmapy/formulary/tests/test_drifts.py b/plasmapy/formulary/tests/test_drifts.py
index 36fedde8..f18c6b05 100644
--- a/plasmapy/formulary/tests/test_drifts.py
+++ b/plasmapy/formulary/tests/test_drifts.py
@@ -1,4 +1,5 @@
 import astropy.units as u
+import numpy as np
 import pytest
 
 from astropy.tests.helper import assert_quantity_allclose
@@ -6,6 +7,62 @@
 from plasmapy.formulary import drifts
 
 
+class Test_diamagnetic_drift:
+    q = 1 * u.C
+    n = 1 / u.m ** 3
+
+    def test_isothermal_plasma(self):
+        r"""
+            Cylindrical & isothermal plasma in Fig 3.4 of Chen.
+        """
+        T = 10 * u.eV
+        B = u.Quantity([0, 0, 100], unit=u.nT)
+        # Lambda = n' / n = 1 / m
+        dp = T * u.Quantity([0, -self.n.value, 0], unit=self.n.unit / u.m)
+
+        result = drifts.diamagnetic_drift(dp, B, self.n, self.q)
+        assert_quantity_allclose(
+            np.linalg.norm(result), T / self.q / np.linalg.norm(B) / u.m
+        )
+
+    def test_diamagnetic_drift_1d_arrays(self):
+        dp = u.Quantity([1, 0, 0], unit=u.N / u.m ** 3)
+        B = u.Quantity([0, -1, 0], unit=u.T)
+        result = drifts.diamagnetic_drift(2 * dp, 3 * B, self.n, self.q)
+        assert_quantity_allclose(result, (2 / 3) * u.Quantity([0, 0, 1], u.m / u.s))
+
+    def test_diamagnetic_drift_2d_array(self):
+        dp = u.Quantity([[1, 0, 0], [1, 0, 0], [1, 0, 0]], unit=u.N / u.m ** 3)
+        B = -u.Quantity([[0, 1, 0], [0, 1, 0], [0, 1, 0]], unit=u.T)
+
+        result = drifts.diamagnetic_drift(2 * dp, 3 * B, self.n, self.q)
+        assert_quantity_allclose(
+            result,
+            (2 / 3) * u.Quantity([[0, 0, 1], [0, 0, 1], [0, 0, 1]], unit=u.m / u.s),
+        )
+
+    def test_diamagnetic_drift_3d_array(self):
+        dp = u.Quantity([[[1, 0, 0]]], unit=u.N / u.m ** 3)
+        B = u.Quantity([[[0, -1, 0]]], unit=u.T)
+
+        result = drifts.diamagnetic_drift(2 * dp, 3 * B, self.n, self.q)
+        assert_quantity_allclose(
+            result, (2 / 3) * u.Quantity([[[0, 0, 1]]], unit=u.m / u.s)
+        )
+
+    def test_nonsensical_units(self):
+        dp = u.Quantity([[1, 0, 0], [1, 0, 0], [1, 0, 0]], unit=u.C)
+        B = u.Quantity([[0, 1, 0], [0, 1, 0], [0, 1, 0]], unit=u.kg)
+        q = 1 * u.C
+        n = 1 / u.m ** 3
+
+        with pytest.raises(u.UnitTypeError):
+            drifts.diamagnetic_drift(dp, B, n, q)
+
+    def test_alias(self):
+        assert drifts.vd_ is drifts.diamagnetic_drift
+
+
 class Test_ExB_drift:
     def test_E_x_B_1d_arrays(self):
         E = u.Quantity([1, 0, 0], unit=u.V / u.m)