Repository: PlasmaPy/PlasmaPy (python)

Issue:
Functions to get magnetostatic fields from known solutions (infinite wires, circular loops)
In addition to the Biot-Savart law solver proposed in #99, we should have a way of calculating the magnetic fields resulting from simple current-carrying wire setups.  This would include the magnetic field from an infinite wire (which has a simple exact solution) and the magnetic field resulting from a circular current-carrying loop (which has a complicated exact solution) for boundary conditions at infinity.  There might be other analytical solutions worth including, but these two are probably the most important ones.

This could perhaps go in the `physics` subpackage in a file called `magnetostatics.py`.  

Pull request description:
Add magnetostatics.py for calculating common magnetostatic fields in `physics` subpackage
Closes #100

Added `MagnetoStatics` class under `physics` subpackage. Implemented concrete magnetostatic fields subclasses and their `magnetic_field` methods to calculate magnetic field at any spatial point. (I use fixed-point integrals to calculate circular and general parametric coil fields. May not be the best solution)

Added `add_magnetostatics` to `Plasma` class.

Example and tests are also added. 

Patch:
diff --git a/docs/physics/index.rst b/docs/physics/index.rst
index ecf12287..342c187e 100644
--- a/docs/physics/index.rst
+++ b/docs/physics/index.rst
@@ -24,6 +24,7 @@ We thus have:
    distribution
    quantum
    relativity
+   magnetostatics
 
 The subpackage makes heavy use of `astropy.units.Quantity` for handling conversions between different unit systems.
 This is especially important for electron volts, commonly used in plasma physics to denote temperature, although
diff --git a/docs/physics/magnetostatics.rst b/docs/physics/magnetostatics.rst
new file mode 100644
index 00000000..4682bce5
--- /dev/null
+++ b/docs/physics/magnetostatics.rst
@@ -0,0 +1,11 @@
+.. _magnetostatics:
+
+**************************************************
+Magnetostatics (`plasmapy.physics.magnetostatics`)
+**************************************************
+
+.. currentmodule:: plasmapy.physics.magnetostatics
+
+.. automodapi:: plasmapy.physics.magnetostatics
+   :no-heading:
+
diff --git a/plasmapy/classes/sources/plasma3d.py b/plasmapy/classes/sources/plasma3d.py
index a2789cb1..24fcab63 100644
--- a/plasmapy/classes/sources/plasma3d.py
+++ b/plasmapy/classes/sources/plasma3d.py
@@ -5,6 +5,9 @@
 
 import numpy as np
 import astropy.units as u
+import itertools
+
+from plasmapy.physics.magnetostatics import MagnetoStatics
 
 from plasmapy.constants import (m_p,
                                 m_e,
@@ -115,3 +118,13 @@ def is_datasource_for(cls, **kwargs):
         else:
             match = False
         return match
+
+    def add_magnetostatic(self, *mstats: MagnetoStatics):
+        # for each MagnetoStatic argument
+        for mstat in mstats:
+            # loop over 3D-index (ix,iy,iz)
+            for point_index in itertools.product(*[list(range(n)) for n in self.domain_shape]):
+                # get coordinate
+                p = self.grid[(slice(None),)+point_index]  # function as [:, *index]
+                # calculate magnetic field at this point and add back
+                self.magnetic_field[(slice(None),)+point_index] += mstat.magnetic_field(p)
diff --git a/plasmapy/examples/plot_magnetic_statics.py b/plasmapy/examples/plot_magnetic_statics.py
new file mode 100644
index 00000000..5e5631fd
--- /dev/null
+++ b/plasmapy/examples/plot_magnetic_statics.py
@@ -0,0 +1,110 @@
+# coding: utf-8
+"""
+Magnetostatic Fields
+=====================
+
+An example of using PlasmaPy's `Magnetostatic` class in `physics` subpackage.
+"""
+
+import plasmapy as pp
+from plasmapy.physics import magnetostatics
+from plasmapy.classes.sources import Plasma3D
+import numpy as np
+import astropy.units as u
+import matplotlib.pyplot as plt
+import itertools
+
+############################################################
+# Some common magnetostatic fields can be generated and added to a plasma object.
+# A dipole
+
+dipole = magnetostatics.MagneticDipole(np.array([0, 0, 1])*u.A*u.m*u.m, np.array([0, 0, 0])*u.m)
+print(dipole)
+
+############################################################
+# initialize a a plasma, where the magnetic field will be calculated on
+
+plasma = Plasma3D(domain_x=np.linspace(-2, 2, 30) * u.m,
+                domain_y=np.linspace(0, 0, 1) * u.m,
+                domain_z=np.linspace(-2, 2, 20) * u.m)
+
+############################################################
+# add the dipole field to it
+plasma.add_magnetostatic(dipole)
+
+X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
+U = plasma.magnetic_field[0, :, 0, :].value.T  # because grid uses 'ij' indexing
+W = plasma.magnetic_field[2, :, 0, :].value.T  # because grid uses 'ij' indexing
+
+
+############################################################
+plt.figure()
+plt.axis('square')
+plt.xlim(-2, 2)
+plt.ylim(-2, 2)
+plt.title('Dipole field in x-z plane, generated by a dipole pointing in the z direction')
+plt.streamplot(plasma.x.value, plasma.z.value, U, W)
+
+############################################################
+# A circular current-carring wire
+
+cw = magnetostatics.CircularWire(np.array([0, 0, 1]), np.array([0, 0, 0])*u.m, 1*u.m, 1*u.A)
+print(cw)
+
+############################################################
+# initialize a a plasma, where the magnetic field will be calculated on
+plasma = Plasma3D(domain_x=np.linspace(-2, 2, 30) * u.m,
+                domain_y=np.linspace(0, 0, 1) * u.m,
+                domain_z=np.linspace(-2, 2, 20) * u.m)
+
+############################################################
+# add the circular coil field to it
+plasma.add_magnetostatic(cw)
+
+X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
+U = plasma.magnetic_field[0, :, 0, :].value.T  # because grid uses 'ij' indexing
+W = plasma.magnetic_field[2, :, 0, :].value.T  # because grid uses 'ij' indexing
+
+plt.figure()
+plt.axis('square')
+plt.xlim(-2, 2)
+plt.ylim(-2, 2)
+plt.title('Circular coil field in x-z plane, generated by a circular coil in the x-y plane')
+plt.streamplot(plasma.x.value, plasma.z.value, U, W)
+
+############################################################
+# a circular wire can be described as parametric equation and converted to GeneralWire
+
+gw_cw = cw.to_GeneralWire()
+
+# the calculated magnetic field is close
+print(gw_cw.magnetic_field([0, 0, 0]) - cw.magnetic_field([0, 0, 0]))
+
+
+############################################################
+# A infinite straight wire
+
+
+iw = magnetostatics.InfiniteStraightWire(np.array([0, 1, 0]), np.array([0, 0, 0])*u.m, 1*u.A)
+print(iw)
+
+############################################################
+# initialize a a plasma, where the magnetic field will be calculated on
+plasma = Plasma3D(domain_x=np.linspace(-2, 2, 30) * u.m,
+                domain_y=np.linspace(0, 0, 1) * u.m,
+                domain_z=np.linspace(-2, 2, 20) * u.m)
+
+# add the infinite straight wire field to it
+plasma.add_magnetostatic(iw)
+
+X, Z = plasma.grid[0, :, 0, :], plasma.grid[2, :, 0, :]
+U = plasma.magnetic_field[0, :, 0, :].value.T  # because grid uses 'ij' indexing
+W = plasma.magnetic_field[2, :, 0, :].value.T  # because grid uses 'ij' indexing
+
+plt.figure()
+plt.title('Dipole field in x-z plane, generated by a infinite straight wire '
+          'pointing in the y direction')
+plt.axis('square')
+plt.xlim(-2, 2)
+plt.ylim(-2, 2)
+plt.streamplot(plasma.x.value, plasma.z.value, U, W)
diff --git a/plasmapy/physics/__init__.py b/plasmapy/physics/__init__.py
index ce6cbd9d..05d2a29b 100644
--- a/plasmapy/physics/__init__.py
+++ b/plasmapy/physics/__init__.py
@@ -36,3 +36,5 @@
 from .relativity import Lorentz_factor
 
 from . import transport
+
+from . import magnetostatics
diff --git a/plasmapy/physics/magnetostatics.py b/plasmapy/physics/magnetostatics.py
new file mode 100644
index 00000000..9bfda4a6
--- /dev/null
+++ b/plasmapy/physics/magnetostatics.py
@@ -0,0 +1,421 @@
+"""
+Define MagneticStatics class to calculate common static magnetic fields
+as first raised in issue #100.
+"""
+
+import abc
+
+import numpy as np
+from astropy import units as u, constants
+from scipy.special import roots_legendre
+
+
+class MagnetoStatics(abc.ABC):
+    """Abstract class for all kinds of magnetic static fields"""
+
+    @abc.abstractmethod
+    def magnetic_field(self, p: u.m) -> u.T:
+        """
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        """
+
+
+class MagneticDipole(MagnetoStatics):
+    """
+    Simple magnetic dipole - two nearby opposite point charges.
+
+    Parameters
+    ----------
+    moment: `astropy.units.Quantity`
+        Magnetic moment vector, in units of A * m^2
+    p0: `astropy.units.Quantity`
+        Position of the dipole
+
+    """
+
+    @u.quantity_input()
+    def __init__(self, moment: u.A * u.m**2, p0: u.m):
+        self.moment = moment.to(u.A*u.m*u.m).value
+        self.p0 = p0.to(u.m).value
+
+    def __repr__(self):
+        return "{name}(moment={moment}, p0={p0})".format(
+            name=self.__class__.__name__,
+            moment=self.moment,
+            p0=self.p0
+        )
+
+    def magnetic_field(self, p: u.m) -> u.T:
+        r"""
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        """
+        r = p - self.p0
+        m = self.moment
+        B = constants.mu0.value/4/np.pi \
+            * (3*r*np.dot(m, r)/np.linalg.norm(r)**5 - m/np.linalg.norm(r)**3)
+        return B*u.T
+
+
+class Wire(MagnetoStatics):
+    """Abstract wire class for concrete wires to be inherited from."""
+
+
+class GeneralWire(Wire):
+    r"""
+    General wire class described by its parametric vector equation
+
+    Parameters
+    ----------
+    parametric_eq: Callable
+        A vector-valued (with units of position) function of a single real
+        parameter.
+    t1: float
+        lower bound of the parameter, smaller than t2
+    t2: float
+        upper bound of the parameter, larger than t1
+    current: `astropy.units.Quantity`
+        electric current
+
+    """
+
+    @u.quantity_input()
+    def __init__(self, parametric_eq,
+                 t1,
+                 t2,
+                 current:
+                 u.A):
+        if callable(parametric_eq):
+            self.parametric_eq = parametric_eq
+        else:
+            raise ValueError("Argument parametric_eq should be a callable")
+        if t1 < t2:
+            self.t1 = t1
+            self.t2 = t2
+        else:
+            raise ValueError(f"t1={t1} is not smaller than t2={t2}")
+        self.current = current.to(u.A).value
+
+    def magnetic_field(self, p: u.m, n: int = 1000) -> u.T:
+        r"""
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+        n : int, optional
+            Number of segments for Wire calculation
+            (defaults to 1000)
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        Notes
+        -----
+        For simplicity, we segment the wire into n equal pieces,
+        and assume each segment is straight. Default n is 1000.
+
+        .. math::
+
+            \vec B
+            \approx \frac{\mu_0 I}{4\pi} \sum_{i=1}^{n}
+            \frac{[\vec l(t_{i}) - \vec l(t_{i-1})] \times
+            \left[\vec p - \frac{\vec l(t_{i}) + \vec l(t_{i-1})}{2}\right]}
+            {\left|\vec p - \frac{\vec l(t_{i}) + \vec l(t_{i-1})}{2}\right|^3},
+            \quad \text{where}\, t_i = t_{\min}+i/n*(t_{\max}-t_{\min})
+
+        """
+
+        p1 = self.parametric_eq(self.t1)
+        step = (self.t2 - self.t1) / n
+        t = self.t1
+        B = 0
+        for i in range(n):
+            t = t + step
+            p2 = self.parametric_eq(t)
+            dl = p2 - p1
+            p1 = p2
+            R = p - (p2 + p1) / 2
+            B += np.cross(dl, R)/np.linalg.norm(R)**3
+        B = B*constants.mu0.value/4/np.pi*self.current
+        return B*u.T
+
+
+class FiniteStraightWire(Wire):
+    """
+    Finite length straight wire class.
+
+    p1 to p2 direction is the possitive current direction.
+
+    Parameters
+    ----------
+    p1: `astropy.units.Quantity`
+        three-dimensional Cartesian coordinate of one end of the straight wire
+    p2: `astropy.units.Quantity`
+        three-dimensional Cartesian coordinate of another end of the straight wire
+    current: `astropy.units.Quantity`
+        electric current
+
+    """
+
+    @u.quantity_input()
+    def __init__(self, p1: u.m, p2: u.m, current: u.A):
+        self.p1 = p1.to(u.m).value
+        self.p2 = p2.to(u.m).value
+        if np.all(p1 == p2):
+            raise ValueError("p1, p2 should not be the same point.")
+        self.current = current.to(u.A).value
+
+    def __repr__(self):
+        return "{name}(p1={p1}, p2={p2}, current={current})".format(
+            name=self.__class__.__name__,
+            p1=self.p1,
+            p2=self.p2,
+            current=self.current
+        )
+
+    def magnetic_field(self, p) -> u.T:
+        r"""
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        Notes
+        -----
+        Let :math:`P_f` be the foot of perpendicular, :math:`\theta_1`(:math:`\theta_2`) be the
+        angles between :math:`\overrightarrow{PP_1}`(:math:`\overrightarrow{PP_2}`)
+        and :math:`\overrightarrow{P_2P_1}`.
+
+        .. math:
+            \vec B = \frac{(\overrightarrow{P_2P_1}\times\overrightarrow{PP_f})^0}
+                     {|\overrightarrow{PP_f}|}
+                     \frac{\mu_0 I}{4\pi} (\cos\theta_1 - \cos\theta_2)
+
+        """
+        # foot of perpendicular
+        p1, p2 = self.p1, self.p2
+        p2_p1 = p2 - p1
+        ratio = np.dot(p - p1, p2_p1)/np.dot(p2_p1, p2_p1)
+        pf = p1 + p2_p1*ratio
+
+        # angles: theta_1 = <p - p1, p2 - p1>, theta_2 = <p - p2, p2 - p1>
+        cos_theta_1 = np.dot(p - p1, p2_p1)/np.linalg.norm(p - p1)/np.linalg.norm(p2_p1)
+        cos_theta_2 = np.dot(p - p2, p2_p1)/np.linalg.norm(p - p2)/np.linalg.norm(p2_p1)
+
+        B_unit = np.cross(p2_p1, p - pf)
+        B_unit = B_unit/np.linalg.norm(B_unit)
+
+        B = B_unit/np.linalg.norm(p-pf)*(cos_theta_1 - cos_theta_2) \
+            * constants.mu0.value/4/np.pi*self.current
+
+        return B*u.T
+
+    def to_GeneralWire(self):
+        """Convert this `Wire` into a `GeneralWire`."""
+        p1, p2 = self.p1, self.p2
+        return GeneralWire(lambda t: p1+(p2-p1)*t, 0, 1, self.current*u.A)
+
+
+class InfiniteStraightWire(Wire):
+    """
+    Infinite straight wire class.
+
+    Parameters
+    ----------
+    direction:
+        three-dimensional direction vector of the wire, also the positive current direction
+    p0: `astropy.units.Quantity`
+        one point on the wire
+    current: `astropy.units.Quantity`
+        electric current
+
+    """
+
+    @u.quantity_input()
+    def __init__(self, direction, p0: u.m, current: u.A):
+        self.direction = direction/np.linalg.norm(direction)
+        self.p0 = p0.to(u.m).value
+        self.current = current.to(u.A).value
+
+    def __repr__(self):
+        return "{name}(direction={direction}, p0={p0}, current={current})".format(
+            name=self.__class__.__name__,
+            direction=self.direction,
+            p0=self.p0,
+            current=self.current
+        )
+
+    def magnetic_field(self, p) -> u.T:
+        r"""
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        Notes
+        -----
+        .. math:
+            \vec B = \frac{\mu_0 I}{2\pi r}*(\vec l^0\times \vec{PP_0})^0,
+            \text{where}\, \vec l^0\, \text{is the unit vector of current direction},
+            r\, \text{is the perpendicular distance between} P_0 \text{and the infinite wire}
+
+        """
+        r = np.cross(self.direction, p - self.p0)
+        B_unit = r / np.linalg.norm(r)
+        r = np.linalg.norm(r)
+
+        return B_unit/r*constants.mu0.value/2/np.pi*self.current*u.T
+
+
+class CircularWire(Wire):
+    """
+    Circular wire(coil) class
+
+    Parameters
+    ----------
+    normal:
+        three-dimensional normal vector of the circular coil
+    center: `astropy.units.Quantity`
+        three-dimensional position vector of the circular coil's center
+    radius: `astropy.units.Quantity`
+        radius of the circular coil
+    current: `astropy.units.Quantity`
+        electric current
+
+    """
+
+    @u.quantity_input()
+    def __init__(self, normal, center: u.m, radius: u.m,
+                 current: u.A, n=300):
+        self.normal = normal/np.linalg.norm(normal)
+        self.center = center.to(u.m).value
+        if radius > 0:
+            self.radius = radius.to(u.m).value
+        else:
+            raise ValueError("Radius should bu larger than 0")
+        self.current = current.to(u.A).value
+
+        # parametric equation
+        # find other two axises in the disc plane
+        z = np.array([0, 0, 1])
+        axis_x = np.cross(z, self.normal)
+        axis_y = np.cross(self.normal, axis_x)
+
+        if np.linalg.norm(axis_x) == 0:
+            axis_x = np.array([1, 0, 0])
+            axis_y = np.array([0, 1, 0])
+        else:
+            axis_x = axis_x/np.linalg.norm(axis_x)
+            axis_y = axis_y/np.linalg.norm(axis_y)
+
+        self.axis_x = axis_x
+        self.axis_y = axis_y
+
+        def curve(t):
+            if isinstance(t, np.ndarray):
+                t = np.expand_dims(t, 0)
+                axis_x_mat = np.expand_dims(axis_x, 1)
+                axis_y_mat = np.expand_dims(axis_y, 1)
+                return self.radius*(np.matmul(axis_x_mat, np.cos(t))
+                                    + np.matmul(axis_y_mat, np.sin(t))) \
+                    + np.expand_dims(self.center, 1)
+            else:
+                return self.radius*(np.cos(t)*axis_x + np.sin(t)*axis_y) + self.center
+        self.curve = curve
+
+        self.roots_legendre = roots_legendre(n)
+        self.n = n
+
+    def __repr__(self):
+        return "{name}(normal={normal}, center={center}, \
+radius={radius}, current={current})".format(
+            name=self.__class__.__name__,
+            normal=self.normal,
+            center=self.center,
+            radius=self.radius,
+            current=self.current
+        )
+
+    def magnetic_field(self, p) -> u.T:
+        r"""
+        Calculate magnetic field generated by this wire at position `p`
+
+        Parameters
+        ----------
+        p : `astropy.units.Quantity`
+            three-dimensional position vector
+
+        Returns
+        -------
+        B : `astropy.units.Quantity`
+            magnetic field at the specified positon
+
+        Notes
+        -----
+        .. math:
+            \vec B
+            = \frac{\mu_0 I}{4\pi}
+            \int \frac{d\vec l\times(\vec p - \vec l(t))}{|\vec p - \vec l(t)|^3}\\
+            = \frac{\mu_0 I}{4\pi} \int_{-\pi}^{\pi} {(-r\sin\theta \hat x + r\cos\theta \hat y)}
+            \times \frac{\vec p - \vec l(t)}{|\vec p - \vec l(t)|^3} d\theta
+
+        We use n points Gauss-Legendre quadrature to compute the integral. The default n is 300.
+
+        """
+
+        x, w = self.roots_legendre
+        t = x*np.pi
+        pt = self.curve(t)
+        dl = self.radius*(
+            - np.matmul(np.expand_dims(self.axis_x, 1), np.expand_dims(np.sin(t), 0))
+            + np.matmul(np.expand_dims(self.axis_y, 1), np.expand_dims(np.cos(t), 0)))  # (3, n)
+
+        r = np.expand_dims(p, 1) - pt  # (3, n)
+        r_norm_3 = np.linalg.norm(r, axis=0)**3
+        ft = np.cross(dl, r, axisa=0, axisb=0)/np.expand_dims(r_norm_3, 1)  # (n, 3)
+
+        return np.pi*np.matmul(np.expand_dims(w, 0), ft).squeeze(0) \
+            * constants.mu0.value/4/np.pi*self.current*u.T
+
+    def to_GeneralWire(self):
+        """Convert this `Wire` into a `GeneralWire`."""
+        return GeneralWire(self.curve, -np.pi, np.pi, self.current*u.A)


Test patch:
diff --git a/plasmapy/classes/sources/tests/test_plasmablob.py b/plasmapy/classes/sources/tests/test_plasmablob.py
index f8a910e5..0b94f37e 100644
--- a/plasmapy/classes/sources/tests/test_plasmablob.py
+++ b/plasmapy/classes/sources/tests/test_plasmablob.py
@@ -2,10 +2,114 @@
 import numpy as np
 import astropy.units as u
 
-from plasmapy.classes.sources import plasmablob
+from plasmapy.classes.sources import plasma3d, plasmablob
+from plasmapy.physics import magnetostatics
 from plasmapy.utils.exceptions import InvalidParticleError
 
-
+@pytest.mark.parametrize('grid_dimensions, expected_size', [
+    ((100, 1, 1), 100),  # Test 1D setup
+    ((128, 128, 1), 16384),  # 2D
+    ((64, 64, 64), 262144),  # 3D
+])
+def test_Plasma3D_setup(grid_dimensions, expected_size):
+    r"""Function to test basic setup of the Plasma3D object.
+
+    Tests that a Plasma3D object initiated with a particular
+    specification behaves in the correct way.
+
+    Parameters
+    ----------
+    grid_dimensions : tuple of ints
+        Grid size of the Plasma3D object to test. Must be a tuple of
+        length 3, indicating length of the grid in x, y, and z
+        directions respectively. Directions not needed should have a
+        length of 1.
+
+    expected_size : int
+        Product of grid dimensions.
+
+    Examples
+    --------
+    >>> test_Plasma3D_setup((10, 10, 10), 1000)
+    >>> test_Plasma3D_setup((100, 10, 1), 1000)
+    """
+    x, y, z = grid_dimensions
+    test_plasma = plasma3d.Plasma3D(domain_x=np.linspace(0, 1, x) * u.m,
+                                    domain_y=np.linspace(0, 1, y) * u.m,
+                                    domain_z=np.linspace(0, 1, z) * u.m)
+
+    # Basic grid setup
+    assert test_plasma.x.size == x
+    assert test_plasma.y.size == y
+    assert test_plasma.z.size == z
+    assert test_plasma.grid.size == 3 * expected_size
+
+    # Core variable units and shapes
+    assert test_plasma.density.size == expected_size
+    assert test_plasma.density.si.unit == u.kg / u.m ** 3
+
+    assert test_plasma.momentum.size == 3 * expected_size
+    assert test_plasma.momentum.si.unit == u.kg / (u.m ** 2 * u.s)
+
+    assert test_plasma.pressure.size == expected_size
+    assert test_plasma.pressure.si.unit == u.Pa
+
+    assert test_plasma.magnetic_field.size == 3 * expected_size
+    assert test_plasma.magnetic_field.si.unit == u.T
+
+    assert test_plasma.electric_field.size == 3 * expected_size
+    assert test_plasma.electric_field.si.unit == u.V / u.m
+
+
+# @pytest.mark.parametrize([()])
+def test_Plasma3D_derived_vars():
+    r"""Function to test derived variables of the Plasma3D class.
+
+    Tests the shapes, units and values of variables derived from core
+    variables.  The core variables are set with arbitrary uniform
+    values.
+    """
+    test_plasma = plasma3d.Plasma3D(domain_x=np.linspace(0, 1, 64) * u.m,
+                                    domain_y=np.linspace(0, 1, 64) * u.m,
+                                    domain_z=np.linspace(0, 1, 1) * u.m)
+
+    # Set an arbitrary uniform values throughout the plasma
+    test_plasma.density[...] = 2.0 * u.kg / u.m ** 3
+    test_plasma.momentum[...] = 10.0 * u.kg / (u.m ** 2 * u.s)
+    test_plasma.pressure[...] = 1 * u.Pa
+    test_plasma.magnetic_field[...] = 0.01 * u.T
+    test_plasma.electric_field[...] = 0.01 * u.V / u.m
+
+    # Test derived variable units and shapes
+    assert test_plasma.velocity.shape == test_plasma.momentum.shape
+    assert (test_plasma.velocity == 5.0 * u.m / u.s).all()
+
+    assert test_plasma.magnetic_field_strength.shape == \
+        test_plasma.magnetic_field.shape[1:]
+    assert test_plasma.magnetic_field_strength.si.unit == u.T
+    assert np.allclose(test_plasma.magnetic_field_strength.value, 0.017320508)
+
+    assert test_plasma.electric_field_strength.shape == \
+        test_plasma.electric_field.shape[1:]
+    assert test_plasma.electric_field_strength.si.unit == u.V / u.m
+
+    assert test_plasma.alfven_speed.shape == test_plasma.density.shape
+    assert test_plasma.alfven_speed.unit.si == u.m / u.s
+    assert np.allclose(test_plasma.alfven_speed.value, 10.92548431)
+
+def test_Plasma3D_add_magnetostatics():
+    r"""Function to test add_magnetostatic function
+    """
+    dipole = magnetostatics.MagneticDipole(np.array([0, 0, 1])*u.A*u.m*u.m, np.array([0, 0, 0])*u.m)
+    cw = magnetostatics.CircularWire(np.array([0, 0, 1]), np.array([0, 0, 0])*u.m, 1*u.m, 1*u.A)
+    gw_cw = cw.to_GeneralWire()
+    iw = magnetostatics.InfiniteStraightWire(np.array([0, 1, 0]), np.array([0, 0, 0])*u.m, 1*u.A)
+    plasma = plasma3d.Plasma3D(domain_x=np.linspace(-2, 2, 30) * u.m,
+                domain_y=np.linspace(0, 0, 1) * u.m,
+                domain_z=np.linspace(-2, 2, 20) * u.m)
+
+    plasma.add_magnetostatic(dipole, cw, gw_cw, iw)
+    
 class Test_PlasmaBlobRegimes:
     def test_intermediate_coupling(self):
         r"""
diff --git a/plasmapy/physics/tests/test_magnetostatics.py b/plasmapy/physics/tests/test_magnetostatics.py
new file mode 100644
index 00000000..ceb2a162
--- /dev/null
+++ b/plasmapy/physics/tests/test_magnetostatics.py
@@ -0,0 +1,154 @@
+import numpy as np
+import pytest
+from astropy import units as u, constants
+
+from plasmapy.physics.magnetostatics import (MagnetoStatics,
+                                             MagneticDipole,
+                                             Wire,
+                                             GeneralWire,
+                                             FiniteStraightWire,
+                                             InfiniteStraightWire,
+                                             CircularWire)
+
+
+mu0_4pi = constants.mu0/4/np.pi
+
+
+class Test_MagneticDipole:
+    def setup_method(self):
+        self.moment = np.array([0, 0, 1])*u.A*u.m*u.m
+        self.p0 = np.array([0, 0, 0])*u.m
+
+    def test_value1(self):
+        "Test a known solution"
+        p = np.array([1, 0, 0])
+        B1 = MagneticDipole(self.moment, self.p0).magnetic_field(p)
+        B1_expected = np.array([0, 0, -1])*1e-7*u.T
+        assert np.all(np.isclose(B1.value, B1_expected.value))
+        assert B1.unit == u.T
+
+    def test_value2(self):
+        "Test a known solution"
+        p = np.array([0, 0, 1])
+        B2 = MagneticDipole(self.moment, self.p0).magnetic_field(p)
+        B2_expected = np.array([0, 0, 2])*1e-7*u.T
+        assert np.all(np.isclose(B2.value, B2_expected.value))
+        assert B2.unit == u.T
+
+
+class Test_GeneralWire:
+    def setup_method(self):
+        self.cw = CircularWire(np.array([0, 0, 1]), np.array([0, 0, 0])*u.m, 1*u.m, 1*u.A)
+        p1 = np.array([0., 0., 0.])*u.m
+        p2 = np.array([0., 0., 1.])*u.m
+        self.fw = FiniteStraightWire(p1, p2, 1*u.A)
+
+    def test_not_callable(self):
+        "Test that `GeneralWire` raises `ValueError` if its first argument is not callale"
+        with pytest.raises(ValueError):
+            GeneralWire("wire", 0, 1, 1*u.A)
+
+    def test_close_cw(self):
+        "Test if the GeneralWire is close to the CircularWire it converted from"
+        gw_cw = self.cw.to_GeneralWire()
+        p = np.array([0, 0, 0])
+        B_cw = self.cw.magnetic_field(p)
+        B_gw_cw = gw_cw.magnetic_field(p)
+
+        assert np.all(np.isclose(B_cw.value, B_gw_cw.value))
+        assert B_cw.unit == B_gw_cw.unit
+
+    def test_close_fw(self):
+        "Test if the GeneralWire is close to the FiniteWire it converted from"
+        gw_fw = self.fw.to_GeneralWire()
+        p = np.array([1, 0, 0])
+        B_fw = self.fw.magnetic_field(p)
+        B_gw_fw = gw_fw.magnetic_field(p)
+
+        assert np.all(np.isclose(B_fw.value, B_gw_fw.value))
+        assert B_fw.unit == B_gw_fw.unit
+    
+    def test_value_error(self):
+        "Test GeneralWire raise ValueError when argument t1>t2"
+        with pytest.raises(ValueError) as e:
+            gw_cw = GeneralWire(lambda t: [0,0,t], 2, 1, 1.*u.A)
+
+
+class Test_FiniteStraightWire:
+    def setup_method(self):
+        self.p1 = np.array([0., 0., -1.])*u.m
+        self.p2 = np.array([0., 0., 1.])*u.m
+        self.current = 1*u.A
+
+    def test_same_point(self):
+        "Test that `FintiteStraightWire` raises `ValueError` if p1 == p2 "
+        with pytest.raises(ValueError):
+            FiniteStraightWire(self.p1, self.p1, self.current)
+
+    def test_value1(self):
+        "Test a known solution"
+        fw = FiniteStraightWire(self.p1, self.p2, self.current)
+        B1 = fw.magnetic_field([1, 0, 0])
+        B1_expected = np.array([0, np.sqrt(2), 0])*1e-7*u.T
+        assert np.all(np.isclose(B1.value, B1_expected.value))
+        assert B1.unit == u.T
+
+    def test_repr(self):
+        "Test __repr__ function"
+        fw = FiniteStraightWire(self.p1, self.p2, self.current)
+        assert repr(fw) == r"FiniteStraightWire(p1=[ 0.  0. -1.], p2=[0. 0. 1.], current=1.0)"
+
+
+
+class Test_InfiniteStraightWire:
+    def setup_method(self):
+        self.direction = np.array([0, 1, 0])
+        self.p0 = np.array([0, 0, 0])*u.m
+        self.current = 1*u.A
+
+    def test_value1(self):
+        "Test a known solution"
+        iw = InfiniteStraightWire(self.direction, self.p0, self.current)
+        B1 = iw.magnetic_field([1, 0, 0])
+        B1_expected = np.array([0, 0, -2])*1e-7*u.T
+        assert np.all(np.isclose(B1.value, B1_expected.value))
+        assert B1.unit == u.T
+
+    def test_repr(self):
+        "Test __repr__ function"
+        iw = InfiniteStraightWire(self.direction, self.p0, self.current)
+        assert repr(iw) == r"InfiniteStraightWire(direction=[0. 1. 0.], p0=[0. 0. 0.], current=1.0)"
+
+class Test_CircularWire:
+    def setup_method(self):
+        self.normalz = np.array([0, 0, 1])
+        self.normalx = np.array([1, 0, 0])
+        self.center = np.array([0, 0, 0])*u.m
+        self.radius = 1*u.m
+        self.current = 1*u.A
+
+    def test_negative_radius(self):
+        "Test that `FintiteStraightWire` raises `ValueError` if radius < 0"
+        with pytest.raises(ValueError):
+            CircularWire(self.normalz, self.center, -1.*u.m, self.current)
+
+    def test_value1(self):
+        "Test a known solution"
+        cw = CircularWire(self.normalz, self.center, self.radius, self.current)
+        B1 = cw.magnetic_field([0, 0, 1])
+        B1_expected = np.array([0, 0, 1])*2*np.pi/2**1.5*1e-7*u.T
+        assert np.all(np.isclose(B1.value, B1_expected.value))
+        assert B1.unit == u.T
+
+    def test_value2(self):
+        "Test a known solution"
+        cw = CircularWire(self.normalx, self.center, self.radius, self.current)
+        B2 = cw.magnetic_field([1, 0, 0])
+        B2_expected = np.array([1, 0, 0])*2*np.pi/2**1.5*1e-7*u.T
+        assert np.all(np.isclose(B2.value, B2_expected.value))
+        assert B2.unit == u.T
+
+    def test_repr(self):
+        "Test __repr__ function"
+        cw = CircularWire(self.normalz, self.center, self.radius, self.current)
+        assert repr(cw) == r"CircularWire(normal=[0. 0. 1.], center=[0. 0. 0.], radius=1.0, current=1.0)"