Add 1D Harris Sheet

changelog/2068.feature.rst

@@ -0,0 +1 @@
+Added the `~plasmapy.plasma.equilibria1d.HarrisSheet` class to calculate magnetic field, current density, and plasma pressure for 1D Harris Sheets.

docs/api_static/plasmapy.plasma.equilibria1d.rst

@@ -0,0 +1,8 @@
+:orphan:
+
+`plasmapy.plasma.equilibria1d`
+==============================
+
+.. currentmodule:: plasmapy.plasma.equilibria1d
+
+.. automodapi::  plasmapy.plasma.equilibria1d

plasmapy/plasma/__init__.py

@@ -7,5 +7,6 @@
 """
 
 from plasmapy.plasma import exceptions, grids, sources
+from plasmapy.plasma.equilibria1d import HarrisSheet
 from plasmapy.plasma.plasma_base import BasePlasma, GenericPlasma
 from plasmapy.plasma.plasma_factory import Plasma

plasmapy/plasma/equilibria1d.py

@@ -0,0 +1,113 @@
+"""Functionality for representing one-dimensional equilibria."""
+
+__all__ = ["HarrisSheet"]
+
+import astropy.constants as const
+import astropy.units as u
+import numpy as np
+
+from plasmapy.utils.decorators.validators import validate_quantities
+
+
+class HarrisSheet:
+    r"""
+    Define a Harris Sheet Equilibrium.
+
+    Magnetic field will be in the :math:`±x` direction and the current
+    density will be in the :math:`±z` direction in a :math:`\hat{x} ×
+    \hat{y} = \hat{z}` coordinate system.
+
+    Parameters
+    ----------
+    B0 : `~astropy.units.Quantity`
+         Magnitude of magnetic field in the limit of :math:`y → ∞` in
+         units convertible to teslas.
+
+    delta : `~astropy.units.Quantity`
+        The thickness of the current sheet in units convertible to
+        meters.
+
+    P0 : `~astropy.units.Quantity`
+        The plasma pressure in the limit of :math:`y → ∞` in units
+        convertible to pascals.
+
+    Notes
+    -----
+    A current sheet is current limited to a surface.
+
+    A Harris sheet is a 1D ideal MHD equilibrium. In resistive MHD if
+    there is any resistivity, it won't be a true equilibrium since the
+    resistivity will gradually smooth the profile out over time.
+
+    A Harris sheet is often used as the initial condition for
+    simulations of magnetic reconnection.
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> harris_sheet = HarrisSheet(delta = 3 * u.m, B0 = 2 * u.T)
+    >>> harris_sheet.magnetic_field(y = 5 * u.m)
+    <Quantity 1.8622... T>
+    """
+
+    def __init__(self, B0, delta, P0=0 * u.Pa):
+        self.B0 = B0
+        self.delta = delta
+        self.P0 = P0
+
+    @validate_quantities
+    def magnetic_field(self, y: u.m) -> u.T:
+        r"""
+        Compute the magnetic field.
+
+        In this equation, :math:`B_0` is the asymptotic magnitude of the
+        magnetic field for :math:`y → ±∞` and :math:`δ` is the thickness
+        of the sheet.
+
+        .. math::
+
+            B_x(y) = B_0 \tanh \left( \frac{y}{δ} \right)
+
+        Parameters
+        ----------
+        y : `~astropy.units.Quantity`
+           Orthogonal distance from the current sheet center.
+        """
+        return self.B0 * np.tanh(u.rad * y / self.delta)
+
+    @validate_quantities
+    def current_density(self, y: u.m) -> u.A / u.m**2:
+        r"""
+        Compute the current density.
+
+        .. math::
+
+          J_z(y) = - \frac{B_0}{δ μ_0) \mathrm{sech}^2 \left( \frac{y}{δ} \right)
+
+        Parameters
+        ----------
+        y : `~astropy.units.Quantity`
+          Orthogonal distance from the current sheet center.
+        """
+        return (
+            -self.B0 / (self.delta * const.mu0) * np.cosh(u.rad * y / self.delta) ** -2
+        )
+
+    @validate_quantities
+    def plasma_pressure(self, y: u.m) -> u.Pa:
+        r"""
+        Compute plasma pressure.
+
+        .. math::
+
+            p(y) = \frac{B_0^2}{2 μ_0} \mathrm{sech}^2 \left( \frac{y}{δ} \right) + p_0
+
+        Parameters
+        ----------
+        y : `~astropy.units.Quantity`
+          Orthogonal distance from the current sheet center.
+        """
+        return (
+            self.B0**2 / (2 * const.mu0) * (np.cosh(u.rad * y / self.delta) ** -2)
+            + self.P0
+        )

plasmapy/plasma/tests/test_equilibria1d.py

@@ -0,0 +1,65 @@
+import astropy.units as u
+import numpy as np
+
+from plasmapy.formulary import magnetic_pressure
+from plasmapy.plasma.equilibria1d import HarrisSheet
+
+
+def test_HarrisSheet():
+    B0 = 1 * u.T
+    delta = 1 * u.m
+    P0 = 0 * u.Pa
+    hs = HarrisSheet(B0, delta, P0)
+    B = hs.magnetic_field(0 * u.m)
+    assert u.isclose(
+        B, 0 * u.T, atol=1e-9 * u.T
+    ), "Magnetic field is supposed to be zero at y=0"
+
+
+def test_HarrisSheet_pressure_balance():
+    B0 = 1 * u.T
+    delta = 1 * u.m
+    P0 = 0 * u.Pa
+    hs = HarrisSheet(B0, delta, P0)
+    y = [-7, -3, 0, 2, 47] * u.m
+    B = hs.magnetic_field(y)
+    P = hs.plasma_pressure(y)
+    p_b = magnetic_pressure(B)
+    total_pressure = P + p_b
+    assert u.allclose(total_pressure, total_pressure[0], atol=1e-9 * u.Pa)
+
+
+def test_HarrisSheet_current_density():
+    B0 = 1 * u.T
+    delta = 1 * u.m
+    P0 = 0 * u.Pa
+    hs = HarrisSheet(B0, delta, P0)
+    y = [-2, 0, 2] * u.m
+    J = hs.current_density(y)
+    correct_J = [-56222.1400445, -795774.715459, -56222.1400445] * u.A / u.m**2
+    assert u.allclose(J, correct_J, atol=1e-8 * u.A / u.m**2)
+
+
+def test_HarrisSheet_magnetic_field():
+    B0 = 1 * u.T
+    delta = 1 * u.m
+    P0 = 0 * u.Pa
+    hs = HarrisSheet(B0, delta, P0)
+    y = [-2, 0, 2] * u.m
+    B = hs.magnetic_field(y)
+    correct_B = [-0.96402758007, 0, 0.96402758007] * u.T
+    assert u.allclose(B, correct_B, atol=1e-9 * u.T)
+
+
+def test_HarrisSheet_limits():
+    y = [-np.inf, np.inf] * u.m
+    B0 = 1 * u.T
+    delta = 1 * u.m
+    P0 = 0 * u.Pa
+    hs = HarrisSheet(B0, delta, P0)
+    B = hs.magnetic_field(y)
+    P = hs.plasma_pressure(y)
+    J = hs.current_density(y)
+    assert u.allclose(B, [-B0, B0], atol=1e-9 * u.T)
+    assert u.allclose(P, [P0, P0], atol=1e-9 * u.Pa)
+    assert u.allclose(J, [0, 0] * u.amp / u.m**2, atol=1e-9 * u.amp / u.m**2)