adding Spheromak solutions

plasmapy/plasma/cylndrical_equilibria.py

@@ -0,0 +1,19 @@
+import scipy.special
+
+
+class ForceFreeFluxRope:
+    def __init__(self, B0, a, lamb):
+        self.B0 = B0
+        self.a = a
+        self.lamb = lamb
+        self.r = r
+        self.z = z
+
+    def B_radial_cyln():
+        return 0
+
+    def B_phi_cyln(self, B, r):
+        return B * scipy.special.j1(self.lamb * r)
+
+    def B_z_cyln(self, B, r):
+        return B * scipy.special.j1(self.lamb * r)

plasmapy/plasma/spheromak_solutions.py

@@ -0,0 +1,113 @@
+import math
+import numpy as np
+import scipy.special
+import sympy.functions.special.bessel
+
+from sympy import Derivative
+
+
+class solution:
+    """
+    Define Analytical solution for spheromak equilibria.
+
+    Parameters
+    ----------
+    B0: `float`
+    magnetic field.
+
+    j1 : `float`
+    spherical bessel function.
+
+    r : 'float'
+    A surface where the radial magnetic field vanishes.
+
+    lamb : 'float'
+    eigenvalue to make j cross B = 0.
+
+    Notes
+    -----
+    A spheromak is an arrangement of plasma that is formed smilar to a smoke ring. The plasma uses its own properties to creat a torodial shape.
+
+
+    """
+
+    def __init__(self, B0, a, lamb):
+        self.B0 = B0
+        self.a = a
+        self.lamb = lamb
+        self.r = r
+        self.z = z
+
+    def B_radial(self, r, theta):
+        r"""
+        Compute the magnetic field in the radial direction.
+
+        .. math::
+
+            2*B_0*(a/r)*j1*(\lambda*r)*\cos(\theta)
+
+        Parameters
+        ----------
+        lamb : `float`
+          eigenvalue to make J cross B = 0.
+
+        """
+
+        return (
+            2
+            * self.B0
+            * (self.a / r)
+            * scipy.special.j1(self.lamb * r)
+            * math.cos(theta)
+        )
+
+    def B_theta(self, r, j1, theta):
+        r"""
+        Compute the magnetic field for theta.
+
+        .. math::
+
+            -1*B_0*(a/r)*Derivative[r*j1*(\lambda*r)]*\sin(\theta)
+
+        Parameters
+        ----------
+        lamb : `float`
+          eigenvalue to make J cross B = 0.
+
+        """
+        return (
+            -1
+            * self.B0
+            * (self.a / r)
+            * Derivative[r * sympy.functions.special.bessel.ji(self.lamb * r), r]
+            * math.sin(theta)
+        )
+
+    def B_phi(self, r, theta):
+        r"""
+        Compute the magnetic field for phi.
+
+        .. math::
+
+        lamb*a*B_0*j1*(\lambda*r)*\sin(\theta)
+
+        Parameters
+        ----------
+        lamb : `float`
+          eigenvalue to make J cross B = 0.
+
+        """
+        return (
+            self.lamb
+            * self.a
+            * self.B0
+            * scipy.special.j1(self.lamb * r)
+            * math.sin(theta)
+        )
+
+    def pshi(self, r, z):
+        beta = self.lamb * (np.sqrt(r**2 + z**2))
+
+        return ((2 * np.pi * (r) ** 2) / (self.lamb * (r**2 + z**2))) * (
+            (np.sin(beta) / beta) - np.cos(beta)
+        )

plasmapy/plasma/test_equilibria_spheromak.py

@@ -0,0 +1,66 @@
+import numpy as np
+
+from astropy import units as u
+
+from plasmapy.formulary import magnetic_pressure
+from plasmapy.plasma.equilibria1d import HarrisSheet
+
+
+def test_Spheromak():
+    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_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_currentDensity():
+    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_magneticField():
+    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_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)