Repository: PlasmaPy/PlasmaPy

Issue:
Improve project description on PyPI
[Project description](https://pypi.org/project/plasmapy/#description) is currently bland. We should detail some functionality, problems that the project is trying to address, mention openAstronomy affiliation, etc.


Discussion:

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)
diff --git a/setup.cfg b/setup.cfg
index 504ccf85..fc4a9ff7 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,7 +1,7 @@
 [metadata]
 package_name = plasmapy
-description = "Python package for plasma physics"
-long_description = "PlasmaPy is a community-developed and community-driven Python package for plasma physics."
+description = Python package for plasma physics
+long_description = PlasmaPy is a community-developed and community-driven generalist Python package for plasma physics tools.
 author = PlasmaPy Community
 license = BSD 3-Clause
 url = http://www.plasmapy.org
@@ -9,10 +9,9 @@ edit_on_github = True
 github_project = PlasmaPy/PlasmaPy
 # install_requires should be formatted as a comma-separated list, e.g.:
 # install_requires = astropy, scipy, matplotlib
-setup_requires = numpy
-install_requires = numpy, scipy, astropy, cython, lmfit, matplotlib, pytest, colorama, roman, mpmath
+install_requires = numpy, scipy, astropy, cython, lmfit, matplotlib, colorama, roman, mpmath
 # version should be PEP386 compatible (http://www.python.org/dev/peps/pep-0386)
-version = 0.1.1
+version = 0.1.1.dev
 
 [build_sphinx]
 source-dir = docs
diff --git a/setup.py b/setup.py
index 780eb6db..46c9ca79 100644
--- a/setup.py
+++ b/setup.py
@@ -23,6 +23,7 @@
 
 from astropy_helpers.setup_helpers import (register_commands, get_debug_option,
                                            get_package_info)
+from astropy_helpers.distutils_helpers import is_distutils_display_option
 from astropy_helpers.git_helpers import get_git_devstr
 from astropy_helpers.version_helpers import generate_version_py
 
@@ -124,6 +125,18 @@
                     os.path.relpath(root, PACKAGENAME), filename))
 package_info['package_data'][PACKAGENAME].extend(c_files)
 
+setup_requires = ['numpy']
+
+# Make sure to have the packages needed for building PlasmaPy, but do not require them
+# when installing from an sdist as the c files are included there.
+if not os.path.exists(os.path.join(os.path.dirname(__file__), 'PKG-INFO')):
+    setup_requires.extend(['cython>=0.27.2'])
+
+# Avoid installing setup_requires dependencies if the user just
+# queries for information
+if is_distutils_display_option():
+    setup_requires = []
+
 # Note that requires and provides should not be included in the call to
 # ``setup``, since these are now deprecated. See this link for more details:
 # https://groups.google.com/forum/#!topic/astropy-dev/urYO8ckB2uM
@@ -132,18 +145,33 @@
       version=VERSION,
       description=DESCRIPTION,
       scripts=scripts,
-      setup_requires=metadata.get("setup_requires", None),
+      setup_requires=[s.strip() for s in metadata.get('install_requires', 'astropy').split(',')],
       install_requires=[s.strip() for s in metadata.get('install_requires', 'astropy').split(',')],
       author=AUTHOR,
       author_email=AUTHOR_EMAIL,
       license=LICENSE,
       url=URL,
       long_description=LONG_DESCRIPTION,
+      keywords=['plasma', 'physics', 'transport', 'collisions', 'science',
+                'atomic', 'particle', 'simulation', 'langmuir', 'tokamak',
+                'instability', 'modeling'],
+      classifiers=[
+          'Intended Audience :: Science/Research',
+          'License :: OSI Approved :: BSD-2-Clause-Patent',
+          'Operating System :: OS Independent',
+          'Programming Language :: C',
+          'Programming Language :: Cython',
+          'Programming Language :: Python :: 3',
+          'Programming Language :: Python :: Implementation :: CPython',
+          'Topic :: Scientific/Engineering :: Astronomy',
+          'Topic :: Scientific/Engineering :: Physics'
+      ],
       cmdclass=cmdclassd,
       zip_safe=False,
       use_2to3=False,
       include_package_data=True,
       entry_points=entry_points,
       python_requires='>={}'.format("3.6"),
+      tests_require=["pytest", "pytest-astropy"],
       **package_info
 )


Test patch:
diff --git a/plasmapy/classes/sources/tests/__init__.py b/plasmapy/classes/sources/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/plasmapy/classes/sources/tests/test_plasma3d.py b/plasmapy/classes/sources/tests/test_plasma3d.py
new file mode 100644
index 00000000..92de1b38
--- /dev/null
+++ b/plasmapy/classes/sources/tests/test_plasma3d.py
@@ -0,0 +1,98 @@
+import pytest
+import numpy as np
+import astropy.units as u
+
+from plasmapy.classes.sources import plasma3d
+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)
diff --git a/plasmapy/classes/tests/test_Plasma.py b/plasmapy/classes/sources/tests/test_plasmablob.py
similarity index 94%
rename from plasmapy/classes/tests/test_Plasma.py
rename to plasmapy/classes/sources/tests/test_plasmablob.py
index e0e901db..0b94f37e 100644
--- a/plasmapy/classes/tests/test_Plasma.py
+++ b/plasmapy/classes/sources/tests/test_plasmablob.py
@@ -3,6 +3,7 @@
 import astropy.units as u
 
 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', [
@@ -96,7 +97,19 @@ def test_Plasma3D_derived_vars():
     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)"