Repository: PlasmaPy/PlasmaPy

Issue:
Add a class to describe 1D Harris sheets
### Feature description

A Harris sheet is a 1D configuration of the form $B_x(y) = B_0 \tanh(\frac{y-y_0}{\delta})$.  There are also formulas for the out-of-plane current density and plasma pressure. We should have a class that can be used to describe Harris sheet equilibria.

### Motivation

In reconnection research, Harris sheets come up all the time (and then some). They are used as standard simplified initial conditions in cartesian geometry.  Having a Harris sheet class would make the lives of reconnection researchers just a little bit easier.

### Implementation strategy

We could create a `HarrisSheet` class that would accept $y_0$, $\delta$, and $B_0$ upon creation of a `HarrisSheet` object.  The class can then have methods that would provide $B_x(y)$, $J_z(y)$, and $p(y)$.  

### Additional context

We should think about where the best location is to put this class.  I don't think it quite fits in the formulary, but maybe we want to create a `plasmapy.equilibria` or `plasmapy.configurations` subpackage. (That would also be where a Grad-Shafranov solver might go, as in #46.) Edited to add: Or perhaps we could put it in `plasmapy.plasma.equilibria`, at least for now?

We could possibly consider force-free equilibria too, where $\mathbf{J}\times\mathbf{B}=0$, which also sometimes come up in reconnection research.

This issue should be addressed before #1978, since we might learn some lessons that would help us when implementing #1978.

Discussion:
github-actions[bot]: Thank you for contributing! ✨ PlasmaPy's future depends on contributors like you, so we deeply appreciate it! 🌱 Please check out the [**contributor guide**](https://docs.plasmapy.org/en/latest/contributing/index.html) for more information on how to contribute.

  The following checklist will help you and the reviewer through the code review process.
   - Overall
     - [ ] Does the PR do what it intends to do?
     - [ ] Except for very minor changes, is a changelog entry included and consistent with the [changelog guide](https://docs.plasmapy.org/en/latest/contributing/changelog_guide.html)?
     - [ ] Are the continuous integration checks passing? (Most linter problems can be automagically fixed by commenting on this PR with `pre-commit.ci autofix`.)
   - Code
     - [ ] Is new/updated code readable and consistent with the [coding guide](https://docs.plasmapy.org/en/latest/contributing/coding_guide.html)?
     - [ ] Can any large functions be split up into shorter functions that each do one thing?
     - [ ] Are SI units used?
   - Tests
     - [ ] Are tests added/updated as required, and consistent with the [testing guide](https://docs.plasmapy.org/en/latest/contributing/testing_guide.html)?
     - [ ] Are the tests understandable?
     - [ ] Do the tests cover all important cases?
   - Docs
     - [ ] Are docs added/updated as required, and consistent with the [doc guide](https://docs.plasmapy.org/en/latest/contributing/doc_guide.html)?
     - [ ] Are the docs understandable?
     - [ ] Do the docs show up correctly in the preview, including Jupyter notebooks?
codecov[bot]: ## [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/2068?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) Report
Patch coverage: **`57.14`**% and project coverage change: **`-0.07`** :warning:
> Comparison is base [(`8c05a40`)](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/8c05a400a66acbd33af4054b8d4101049cdc43da?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) 98.35% compared to head [(`d010902`)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/2068?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) 98.29%.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #2068      +/-   ##
==========================================
- Coverage   98.35%   98.29%   -0.07%     
==========================================
  Files          96       97       +1     
  Lines        8528     8542      +14     
==========================================
+ Hits         8388     8396       +8     
- Misses        140      146       +6     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/2068?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) | Coverage Δ | |
|---|---|---|
| [plasmapy/plasma/equilibria1d.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/2068?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvcGxhc21hL2VxdWlsaWJyaWExZC5weQ==) | `57.14% <57.14%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)

</details>

[:umbrella: View full report in Codecov by Sentry](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/2068?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).

namurphy: Here are a couple of possibilities for a few more tests of `HarrisSheet`.

 - One possible test would be to pick ∼3 values of $y$, and then calculate $B_x(y)$, $J_z(y)$, and $p(y)$ both manually and using `HarrisSheet.magnetic_field`, etc.  Then the test can check that the results calculated using `HarrisSheet` match the values calculated manually using `u.isclose` or `u.allclose`.

 - A second possibility would be to make sure that $B_x$, $J_z$, and $p$ behave as expected for $y→±∞$ m.  
    - $B_x$ will approach $B_0$ as $y→\infty$ m, and $-B_0$ as $y→-\infty$ m.  
    - $J_z$ will approach $0$ A m$^{-2}$ as $y→±\infty$ m.
    - $p$ will approach $p_0$.  
   We can represent $\infty$ as `np.inf` (or in this case `np.inf * u.m`).  

Thank you again for doing this...we really appreciate it!
namurphy: pre-commit.ci autofix
namurphy: pre-commit.ci autofix

Patch:
diff --git a/changelog/2068.feature.rst b/changelog/2068.feature.rst
new file mode 100644
index 00000000..4d153f28
--- /dev/null
+++ b/changelog/2068.feature.rst
@@ -0,0 +1,1 @@
+Added the `~plasmapy.plasma.equilibria1d.HarrisSheet` class to calculate magnetic field, current density, and plasma pressure for 1D Harris Sheets.
diff --git a/docs/api_static/plasmapy.plasma.equilibria1d.rst b/docs/api_static/plasmapy.plasma.equilibria1d.rst
new file mode 100644
index 00000000..b6e7364c
--- /dev/null
+++ b/docs/api_static/plasmapy.plasma.equilibria1d.rst
@@ -0,0 +1,8 @@
+:orphan:
+
+`plasmapy.plasma.equilibria1d`
+==============================
+
+.. currentmodule:: plasmapy.plasma.equilibria1d
+
+.. automodapi::  plasmapy.plasma.equilibria1d
diff --git a/plasmapy/plasma/__init__.py b/plasmapy/plasma/__init__.py
index 5d10c103..2ad644bf 100644
--- a/plasmapy/plasma/__init__.py
+++ b/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
diff --git a/plasmapy/plasma/equilibria1d.py b/plasmapy/plasma/equilibria1d.py
new file mode 100644
index 00000000..f4541549
--- /dev/null
+++ b/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
+        )


Test patch:
diff --git a/plasmapy/plasma/tests/test_equilibria1d.py b/plasmapy/plasma/tests/test_equilibria1d.py
new file mode 100644
index 00000000..6a0d6580
--- /dev/null
+++ b/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)