Repository: PlasmaPy/PlasmaPy
Pull request #459 (2018-05-16 06:09:58)

[WIP] Basic plasma factory implementation
We can now instantiate the expected plasma subclass based on input keyword arguments.

```python
>>> import astropy.units as u
>>> import numpy as np
>>> import plasmapy.classes

>>> T_e = 25 * 15e3 * u.K
>>> n_e = 1e26 * u.cm ** -3
>>> Z = 2.0
>>> particle = 'p'
>>> blob = plasmapy.classes.Plasma(T_e=T_e,
...                                n_e=n_e,
...                                Z=Z,
...                                particle=particle)

>>> type(blob)
plasmapy.classes.sources.plasmablob.PlasmaBlob

>>> three_dims = plasmapy.classes.Plasma(domain_x=np.linspace(0, 1, 3) * u.m,
...                                      domain_y=np.linspace(0, 1, 3) * u.m,
...                                      domain_z=np.linspace(0, 1, 3) * u.m)

>>> type(three_dims)
plasmapy.classes.sources.plasma3d.Plasma3D
```

I am not sure if relying solely on keyword arguments is good enough to assume the plasma subclass but I don't have a better idea either (SunPy uses different parameters for `data` and `header`, I don't know how well will it work in our case).

We also need to add common plasma data to `GenericPlasma` and as of now `PlasmaFactory` does nothing but just inherits from `BasicRegistrationFactory`.


Patch:
diff --git a/CITATION.md b/CITATION.md
new file mode 100644
index 00000000..a2ca3f1d
--- /dev/null
+++ b/CITATION.md
@@ -0,0 +1,9 @@
+If you use PlasmaPy for work or research resulting in a publication, please use the following acknowledgement:
+  
+* This project made use of PlasmaPy, a community-developed core Python package for plasma physics.
+
+The earliest citation for PlasmaPy is:
+
+* Murphy, Nicholas A, Huang, Yi-Min, & PlasmaPy Community. (2016, October). PlasmaPy: beginning a community developed Python package for plasma physics. Zenodo. http://doi.org/10.5281/zenodo.163752
+    
+We highly encourage researchers to acknowledge the packages that PlasmaPy depends on, including but not limited to Astropy, NumPy, and SciPy.
diff --git a/CITATION.rst b/CITATION.rst
deleted file mode 100644
index 0d07b9a1..00000000
--- a/CITATION.rst
+++ /dev/null
@@ -1,23 +0,0 @@
-Acknowledging and Citing PlasmaPy
-=================================
-
-If you use PlasmaPy for a project resulting in a publication, we ask
-that you cite the following reference
-(`BibTeX <https://zenodo.org/record/1238132/export/hx#.WvMkQK0cChc>`_):
-
-* PlasmaPy Community, Nicholas A. Murphy, Andrew J. Leonard, Dominik
-  Stańczak, Pawel M. Kozlowski, Samuel J. Langendorf, Colby C. Haggerty,
-  Jasper P. Beckers, Stuart J. Mumford, Tulasi N. Parashar, and Yi-Min
-  Huang. (2018, April). PlasmaPy: an open source community-developed
-  Python package for plasma physics. Zenodo.
-  http://doi.org/10.5281/zenodo.1238132
-
-We provide the following standard acknowledgment with this citation that
-you may use instead of a citation in the body of a paper.
-
-* This research made use of PlasmaPy, a community-developed open source
-  core Python package for plasma physics (PlasmaPy Community 2018).
-
-We highly encourage researchers to acknowledge the packages that
-PlasmaPy depends on, including but not limited to Astropy, NumPy, and
-SciPy.
diff --git a/README.md b/README.md
index f10c0f1a..6e1f3f9b 100644
--- a/README.md
+++ b/README.md
@@ -1,63 +1,49 @@
 # PlasmaPy
-
-[![PyPI version](https://badge.fury.io/py/plasmapy.svg)](https://badge.fury.io/py/plasmapy)
 [![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](./LICENSE.md)
-[![astropy](http://img.shields.io/badge/powered%20by-AstroPy-orange.svg?style=flat)](http://www.astropy.org/)
-[![Matrix](https://matrix.to/img/matrix-badge.svg)](https://riot.im/app/#/room/#plasmapy:matrix.org)
-[![Gitter](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/PlasmaPy/Lobby)
-
 [![Build Status](https://travis-ci.org/PlasmaPy/PlasmaPy.svg?branch=master)](https://travis-ci.org/PlasmaPy/PlasmaPy)
 [![Build status](https://ci.appveyor.com/api/projects/status/hbduy62sqrvy8rn7?svg=true)](https://ci.appveyor.com/project/namurphy/plasmapy)
+[![Documentation Status](https://readthedocs.org/projects/plasmapy/badge/?version=latest)](http://docs.plasmapy.org/en/latest/?badge=latest)
 [![Coverage Status](https://coveralls.io/repos/github/PlasmaPy/PlasmaPy/badge.svg?branch=master)](https://coveralls.io/github/PlasmaPy/PlasmaPy?branch=master)
-[![Documentation Status](https://readthedocs.org/projects/plasmapy/badge/?version=latest)](http://plasmapy.readthedocs.io/en/latest/?badge=latest)
-
-PlasmaPy is an open source community developed Python 3.6+ package for
-plasma physics in the early stages of development.  PlasmaPy intends to
-be for plasmas what [Astropy](https://github.com/astropy/astropy) is for
-astronomy - a collection of functionality commonly used and shared
-between plasma physicists and researchers globally, running within and
-leveraging the open source scientific Python ecosystem.  The goals of
-this project are more thoroughly described in our [vision
-statement](http://docs.plasmapy.org/en/stable/about/vision_statement.html)
-and [this recent reference](https://doi.org/10.5281/zenodo.1238132).
+[![Matrix](https://matrix.to/img/matrix-badge.svg)](https://riot.im/app/#/room/#plasmapy:matrix.org)
+[![Gitter](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/PlasmaPy/Lobby)
+
+An open source community developed Python 3.6+ package for plasma physics in
+the early stages of development. PlasmaPy intends to be for plasmas what
+[Astropy](https://github.com/astropy/astropy) is for astronomy - a
+collection of functionality commonly used and shared between plasma physicists
+and researchers globally, running within and leveraging the open source
+scientific Python ecosystem. The goals of this project are more thoroughly described in our
+[vision statement](http://docs.plasmapy.org/en/stable/about/vision_statement.html)
+and [an earlier conference poster](http://doi.org/10.5281/zenodo.163752).
 We are in the process of writing [online
 documentation](http://docs.plasmapy.org/en/latest/).
 
-We created a guide on [contributing to
-PlasmaPy](http://plasmapy.readthedocs.io/en/master/CONTRIBUTING.html)
-and have a [Code of
-Conduct](http://plasmapy.readthedocs.io/en/master/CODE_OF_CONDUCT.html).
+We created a guide on
+[contributing to PlasmaPy](http://docs.plasmapy.org/en/master/CONTRIBUTING.html) and have a
+[Code of Conduct](http://docs.plasmapy.org/en/master/CODE_OF_CONDUCT.html).
 New contributors are very welcome!
 
-If you have any questions, the quickest way to get a response is to ask
-on our
+If you have any questions, the quickest way to get a response is to ask on our
 [Matrix](https://riot.im/app/#/room/#plasmapy:matrix.org)/[Gitter](https://gitter.im/PlasmaPy/Lobby)
-channel.  For a slower response, you may also contact Nick Murphy at
-<namurphy@cfa.harvard.edu> and Drew Leonard at
-<andy.j.leonard@gmail.com>.  You may also [sign up for PlasmaPy's email
-list](https://groups.google.com/forum/#!forum/plasmapy).
+channel.  You may also contact
+Nick Murphy at <namurphy@cfa.harvard.edu> and
+Drew Leonard at <andy.j.leonard@gmail.com>.
 
 ## Installation
 
-You can get PlasmaPy from pip via `pip install plasmapy`. To contribute
-to the package, check out [our instructions on installing PlasmaPy from
-source](INSTALL.md).
+You can get PlasmaPy from pip via `pip install plasmapy`. To contribute to the package, check
+out [our instructions on installing PlasmaPy from source](INSTALL.md).
 
 We're not on Conda yet, but we're working on it!
 
 Like most scientific Python packages, PlasmaPy probably runs best on the
 [Anaconda distribution](https://www.anaconda.com/downloads).
 
-PlasmaPy requires Python 3.6+ and is [not compatible with
-Python 2](https://pythonclock.org/).
+PlasmaPy requires Python 3.6+ and is
+[not compatible with Python 2](https://pythonclock.org/).
 
 ## License
 
 PlasmaPy is licensed under a 3-clause BSD license with added protections
 against software patents - see the [``LICENSE.md``](LICENSE.md) file in
 the top-level directory.
-
-## Acknowledging and Citing PlasmaPy
-
-Please see the [``CITATION.rst``](CITATION.rst) file for instructions on
-how to cite PlasmaPy.
diff --git a/docs/about/citation.rst b/docs/about/citation.rst
deleted file mode 100644
index 70504a4b..00000000
--- a/docs/about/citation.rst
+++ /dev/null
@@ -1,25 +0,0 @@
-.. _citation:
-
-Acknowledging and Citing PlasmaPy
-=================================
-
-If you use PlasmaPy for a project resulting in a publication, we ask
-that you cite the following reference
-(`BibTeX <https://zenodo.org/record/1238132/export/hx#.WvMkQK0cChc>`_):
-
-* PlasmaPy Community, Nicholas A. Murphy, Andrew J. Leonard, Dominik
-  Stańczak, Pawel M. Kozlowski, Samuel J. Langendorf, Colby C. Haggerty,
-  Jasper P. Beckers, Stuart J. Mumford, Tulasi N. Parashar, and Yi-Min
-  Huang. (2018, April). PlasmaPy: an open source community-developed
-  Python package for plasma physics. Zenodo.
-  http://doi.org/10.5281/zenodo.1238132
-
-We provide the following standard acknowledgment with this citation that
-you may use instead of a citation in the body of a paper.
-
-* This research made use of PlasmaPy, a community-developed open source
-  core Python package for plasma physics (PlasmaPy Community 2018).
-
-We highly encourage researchers to acknowledge the packages that
-PlasmaPy depends on, including but not limited to Astropy, NumPy, and
-SciPy.
diff --git a/docs/about/index.rst b/docs/about/index.rst
index 7bcdcfbd..c7255f3c 100644
--- a/docs/about/index.rst
+++ b/docs/about/index.rst
@@ -12,4 +12,3 @@ About PlasmaPy
    release_notes
    change_log
    vision_statement
-   citation
diff --git a/docs/index.rst b/docs/index.rst
index cb4a0700..f5844265 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -11,6 +11,7 @@ PlasmaPy Documentation
    :align: right
    :scale: 80%
 
+
 `PlasmaPy <http://www.plasmapy.org/>`_ is an open source
 community-developed core `Python <https://www.python.org/>`_ 3.6+
 package for plasma physics in the early stages of development.
@@ -30,7 +31,6 @@ Getting Started
 * `PlasmaPy website
   <http://www.plasmapy.org/>`_
 * `Using astropy.units <http://docs.astropy.org/en/stable/units/>`_
-* :ref:`citation`
 
 .. _toplevel-user-documentation:
 
@@ -126,15 +126,12 @@ documentation.
 Project Details
 ***************
 
-.. The about PlasmaPy section has some important information that would
-   be helpful to have more readily accessible from the main doc index
-   page.
-
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
 
    about/index
 
+
 .. _toplevel-index:
 
 *****
diff --git a/docs/plasma/index.rst b/docs/plasma/index.rst
index be47a4d0..bce0189a 100644
--- a/docs/plasma/index.rst
+++ b/docs/plasma/index.rst
@@ -1,32 +1,123 @@
-.. _plasma:
+.. _classes:
 
-****************************************
-Plasma class (`plasmapy.classes.plasma`)
-****************************************
+===============
+PlasmaPy Plasma
+===============
 
-Introduction
-============
+.. module:: plasmapy.classes
+.. currentmodule:: plasmapy.classes
 
-The :class:`Plasma3D` class is a basic structure to contain spatial
-information about a plasma.  To initialize a Plasma3D system, first
-create an instance of the :class:`Plasma3D` class and then set the
-:attr:`Plasma3D.density`, :attr:`Plasma3D.momentum`,
-:attr:`Plasma3D.pressure` and the :attr:`Plasma3D.magnetic_field`.
+Overview
+--------
+One of core classes in PlasmaPy is a `Plasma`. In order to make it easy to work
+with different plasma data in PlasmaPy, the Plasma object provides a number of
+methods for commonly existing plasmas in nature.
 
-This feature is currently under development.
+All Plasma objects are created using the Plasma factory `~plasmapy.classes.Plasma`.
 
-The :class:`PlasmaBlob` class is a basic structure to contain just
-plasma parameter information about a plasma with no associated
-spatial or temporal scales.  To initialize a PlasmaBlob system, call
-it with arguments: electron temperature :attr:`PlasmaBlob.T_e`, 
-and electron density :attr:`PlasmaBlob.n_e`. You may also optionally
-define the ionization, :attr:`PlasmaBlob.Z`, and relevant plasma 
-particle, :attr:`PlasmaBlob.particle`.
+A number of plasma data structures are supported by subclassing this base object. See
+:ref:`plasma-sources` to see a list of all of them.
 
-This feature is currently under heavy development.
 
-Reference/API
-=============
+Creating Plasma Objects
+-----------------------
+Plasma objects are constructed using the special factory
+class `~plasmapy.classes.Plasma`: ::
 
-.. automodapi:: plasmapy.classes.plasma
-   :no-heading:
+    >>> x = plasmapy.classes.Plasma(T_e=T_e,
+                                    n_e=n_e,
+                                    Z=Z,
+                                    particle=particle)  # doctest: +SKIP
+
+The result of a call to `~plasmapy.classes.Plasma` will be either a
+`~plasmapy.classes.GenericPlasma` object, or a subclass of `~plasmapy.classes.GenericPlasma`
+which deals with a specific type of data,
+e.g. `~plasmapy.classes.sources.PlasmaBlob` or `~plasmapy.classes.sources.Plasma3D`
+(see :ref:`plasma-sources` to see a list of all of them).
+
+
+.. autoclass:: plasmapy.classes.plasma_factory.PlasmaFactory
+
+
+Using Plasma Objects
+--------------------
+
+Once a Plasma object has been created using `~plasmapy.classes.Plasma` it will be a instance
+or a subclass of the `~plasmapy.classes.GenericPlasma` class. The documentation of
+`~plasmapy.classes.GenericPlasma` lists the attributes and methods that are available
+on all Plasma objects.
+
+.. _plasma-classes:
+
+Plasma Classes
+--------------
+Defined in `plasmapy.classes.sources` are a set of `~plasmapy.classes.GenericPlasma` subclasses
+which convert the keyword arguments data to the standard `~plasmapy.classes.GenericPlasma` interface.
+These subclasses also provide a method, which describes to the `Plasma <plasmapy.classes.plasma_factory.PlasmaFactory>` factory
+which the data match its plasma data structure.
+
+.. automodapi:: plasmapy.classes
+    :no-main-docstr:
+    :no-heading:
+
+.. _plasma-sources:
+
+Plasma Subclasses
+-----------------
+
+.. toctree::
+    sources
+
+
+Writing a new Plasma subclass
+--------------------------------
+
+Any subclass of `~plasmapy.classes.GenericPlasma` which defines a method named
+`~plasmapy.classes.GenericPlasma.is_datasource_for` will automatically be registered with
+the `Plasma <plasmapy.classes.plasma_factory.PlasmaFactory>` factory.
+The ``is_datasource_for`` method describes the form of the data for which the
+`~plasmapy.classes.GenericPlasma` subclass is valid. For
+example it might check the number and types of keyword arguments.
+This makes it straightforward to define your own
+`~plasmapy.classes.GenericPlasma` subclass for a new data structure or a custom data source
+like simulated data. These classes only have to be imported for this to work, as
+demonstrated by the following example.
+
+.. code-block:: python
+
+    import plasmapy.classes
+    import astropy.units as u
+
+    class FuturePlasma(plasmapy.classes.GenericPlasma):
+        def __init__(self, **kwargs):
+
+            super(FuturePlasma, self).__init__(**kwargs)
+
+        # Specify a classmethod that determines if the input data matches
+        # this new subclass
+        @classmethod
+        def is_datasource_for(cls, **kwargs):
+            """
+            Determines if any of keyword arguments have a dimensionless value.
+            """
+            for _, value in kwargs.items():
+                try:
+                    if value.unit == u.dimensionless_unscaled:
+                        return True
+                except AttributeError:
+                    pass
+
+            return False
+
+
+This class will now be available through the `Plasma <plasmapy.classes.plasma_factory.PlasmaFactory>` factory as long as this
+class has been defined, i.e. imported into the current session.
+
+If you do not want to create a method named ``is_datasource_for`` you can
+manually register your class and matching method using the following method
+
+.. code-block:: python
+
+    import plasmapy.classes
+
+    plasmapy.classes.Plasma.register(FuturePlasma, FuturePlasma.some_matching_method)
diff --git a/docs/plasma/sources.rst b/docs/plasma/sources.rst
new file mode 100644
index 00000000..9f0f1b52
--- /dev/null
+++ b/docs/plasma/sources.rst
@@ -0,0 +1,32 @@
+.. _plasma:
+
+**********************************************
+Plasma Subclasses (`plasmapy.classes.sources`)
+**********************************************
+
+Introduction
+============
+
+The :class:`Plasma3D` class is a basic structure to contain spatial
+information about a plasma.  To initialize a Plasma3D system, first
+create an instance of the :class:`Plasma3D` class and then set the
+:attr:`Plasma3D.density`, :attr:`Plasma3D.momentum`,
+:attr:`Plasma3D.pressure` and the :attr:`Plasma3D.magnetic_field`.
+
+This feature is currently under development.
+
+The :class:`PlasmaBlob` class is a basic structure to contain just
+plasma parameter information about a plasma with no associated
+spatial or temporal scales.  To initialize a PlasmaBlob system, call
+it with arguments: electron temperature :attr:`PlasmaBlob.T_e`, 
+and electron density :attr:`PlasmaBlob.n_e`. You may also optionally
+define the ionization, :attr:`PlasmaBlob.Z`, and relevant plasma 
+particle, :attr:`PlasmaBlob.particle`.
+
+This feature is currently under heavy development.
+
+Reference/API
+=============
+
+.. automodapi:: plasmapy.classes.sources
+   :no-heading:
diff --git a/plasmapy/__init__.py b/plasmapy/__init__.py
index 11b9b593..c15676c2 100644
--- a/plasmapy/__init__.py
+++ b/plasmapy/__init__.py
@@ -15,7 +15,6 @@
 __doc__ = ("A community-developed and community-driven open source "
            "core Python package for plasma physics.")
 
-
 class UnsupportedPythonError(Exception):
     pass
 
@@ -23,6 +22,7 @@ class UnsupportedPythonError(Exception):
 if sys.version_info < tuple((int(val) for val in "3.6".split('.'))):
     raise UnsupportedPythonError("plasmapy does not support Python < {}".format(3.6))
 
+
 if not _ASTROPY_SETUP_:
     # For egg_info test builds to pass, put package imports here.
     from . import atomic
@@ -33,45 +33,3 @@ class UnsupportedPythonError(Exception):
     from . import physics
     from . import utils
 
-def online_help(query):
-    """
-    Search the online PlasmaPy documentation for the given query from plasmapy.org
-    Opens the results in the default web browser.
-    Requires an active Internet connection.
-    Redirects to Astropy.units in case of query 'unit' or 'units'
-
-    Parameters
-    ----------
-    query : str
-        The search query.
-    """
-    from urllib.parse import urlencode
-    import webbrowser
-
-    url = 'http://docs.plasmapy.org/en/stable/search.html?\
-    {0}&check_keywords=yes&area=default'.format(urlencode({'q': query}))
-
-    if(query.lower() in ('unit', 'units')):
-        url = 'http://docs.astropy.org/en/stable/units/'
-
-    webbrowser.open(url)
-
-__citation__ = """@misc{plasmapy_community_2018_1238132,
-  author       = {PlasmaPy Community and
-                  Murphy, Nicholas A. and
-                  Leonard, Andrew J. and
-                  Sta\'nczak, Dominik and
-                  Kozlowski, Pawel M. and
-                  Langendorf, Samuel J. and
-                  Haggerty, Colby C. and
-                  Beckers, Jasper P. and
-                  Mumford, Stuart J. and
-                  Parashar, Tulasi N. and
-                  Huang, Yi-Min},
-  title        = {{PlasmaPy: an open source community-developed 
-                   Python package for plasma physics}},
-  month        = apr,
-  year         = 2018,
-  doi          = {10.5281/zenodo.1238132},
-  url          = {https://doi.org/10.5281/zenodo.1238132}
-}"""
\ No newline at end of file
diff --git a/plasmapy/classes/__init__.py b/plasmapy/classes/__init__.py
index 7b2c2a0e..6b2e2527 100644
--- a/plasmapy/classes/__init__.py
+++ b/plasmapy/classes/__init__.py
@@ -1,2 +1,4 @@
-from .plasma import Plasma3D, PlasmaBlob
+from .plasma_base import BasePlasma, GenericPlasma
+from . import sources
+from .plasma_factory import Plasma
 from .species import Species
diff --git a/plasmapy/classes/plasma_base.py b/plasmapy/classes/plasma_base.py
new file mode 100644
index 00000000..dcb65b41
--- /dev/null
+++ b/plasmapy/classes/plasma_base.py
@@ -0,0 +1,80 @@
+from abc import ABC, abstractmethod, abstractproperty
+
+__all__ = [
+    "BasePlasma",
+    "GenericPlasma",
+]
+
+class BasePlasma(ABC):
+    """
+    Registration class for `~plasmapy.classes.GenericPlasma` and declares
+    some abstract methods for data common in different kinds of plasmas.
+
+    This class checks for the existance of a method named ``is_datasource_for``
+    when a subclass of `GenericPlasma` is defined. If it exists it will add that
+    class to the registry.
+    """
+    # GenericPlasma subclass registry
+    _registry = dict()
+
+    def __init_subclass__(cls, **kwargs):
+        super().__init_subclass__(**kwargs)
+        if hasattr(cls, 'is_datasource_for'):
+            cls._registry[cls] = cls.is_datasource_for
+
+    # This class is supposed to declare abstract methods (@abstractmethod or
+    # @abstractproperty as appropriate) that are common in most plasmas
+    # (like `electronTemperature`, `ionTemperature`,
+    # `electronDensitry`, `ionDensity`, `averageIonization`, etc.)
+    # where as `GenericPlasma` class will hold the definitions for these
+    # abstract methods.
+
+    # For reference, see
+    # https://github.com/sunpy/ndcube/blob/master/ndcube/ndcube.py#L26
+
+    @abstractproperty
+    def electron_temperature(self):
+        pass
+
+    @abstractproperty
+    def ion_temperature(self):
+        pass
+
+    @abstractproperty
+    def electron_density(self):
+        pass
+
+    @abstractproperty
+    def ion_density(self):
+        pass
+
+    @abstractproperty
+    def average_ionization(self):
+        pass
+
+
+class GenericPlasma(BasePlasma):
+    """
+    A Generic Plasma class. This class contains definitions for abstract
+    methods declared in the `~plasmapy.classes.plasma_base.BaseClass`.
+    """
+    def __init__(self, **kwargs):
+        pass
+
+    # The definitions for the abstract methods declared in `BasePlasma`
+    # goes here.
+
+    def electron_temperature(self):
+        pass
+
+    def ion_temperature(self):
+        pass
+
+    def electron_density(self):
+        pass
+
+    def ion_density(self):
+        pass
+
+    def average_ionization(self):
+        pass
diff --git a/plasmapy/classes/plasma_factory.py b/plasmapy/classes/plasma_factory.py
new file mode 100644
index 00000000..d4cc79e4
--- /dev/null
+++ b/plasmapy/classes/plasma_factory.py
@@ -0,0 +1,24 @@
+from plasmapy.classes.plasma_base import GenericPlasma
+
+from plasmapy.utils.datatype_factory_base import BasicRegistrationFactory
+from plasmapy.utils.datatype_factory_base import NoMatchError
+from plasmapy.utils.datatype_factory_base import MultipleMatchError
+from plasmapy.utils.datatype_factory_base import ValidationFunctionError
+
+__all__ = [
+    "PlasmaFactory",
+    "Plasma",
+]
+
+class PlasmaFactory(BasicRegistrationFactory):
+    """
+    Plasma factory class. Used to create a variety of Plasma objects.
+    Valid plasma structures are specified by registering them with the
+    factory.
+    """
+    pass
+
+
+Plasma = PlasmaFactory(default_widget_type=GenericPlasma,
+                       registry=GenericPlasma._registry,
+                       additional_validation_functions=['is_datasource_for'])
diff --git a/plasmapy/classes/sources/__init__.py b/plasmapy/classes/sources/__init__.py
new file mode 100644
index 00000000..8f01d98d
--- /dev/null
+++ b/plasmapy/classes/sources/__init__.py
@@ -0,0 +1,2 @@
+from .plasma3d import Plasma3D
+from .plasmablob import PlasmaBlob
diff --git a/plasmapy/classes/sources/plasma3d.py b/plasmapy/classes/sources/plasma3d.py
new file mode 100644
index 00000000..a2789cb1
--- /dev/null
+++ b/plasmapy/classes/sources/plasma3d.py
@@ -0,0 +1,117 @@
+"""
+Defines the core Plasma class used by PlasmaPy to represent plasma properties.
+"""
+import warnings
+
+import numpy as np
+import astropy.units as u
+
+from plasmapy.constants import (m_p,
+                                m_e,
+                                c,
+                                mu0,
+                                k_B,
+                                e,
+                                eps0,
+                                pi,
+                                )
+
+from plasmapy.classes import GenericPlasma
+
+__all__ = [
+    "Plasma3D"
+]
+
+
+class Plasma3D(GenericPlasma):
+    """
+    Core class for describing and calculating plasma parameters with
+    spatial dimensions.
+
+    Attributes
+    ----------
+    x : `astropy.units.Quantity`
+        x-coordinates within the plasma domain. Equal to the
+        `domain_x` input parameter.
+    y : `astropy.units.Quantity`
+        y-coordinates within the plasma domain. Equal to the
+        `domain_y` input parameter.
+    z : `astropy.units.Quantity`
+        z-coordinates within the plasma domain. Equal to the
+        `domain_z` input parameter.
+    grid : `astropy.units.Quantity`
+        (3, x, y, z) array containing the values of each coordinate at
+        every point in the domain.
+    domain_shape : tuple
+        Shape of the plasma domain.
+    density : `astropy.units.Quantity`
+        (x, y, z) array of mass density at every point in the domain.
+    momentum : `astropy.units.Quantity`
+        (3, x, y, z) array of the momentum vector at every point in
+        the domain.
+    pressure : `astropy.units.Quantity`
+        (x, y, z) array of pressure at every point in the domain.
+    magnetic_field : `astropy.units.Quantity`
+        (3, x, y, z) array of the magnetic field vector at every point
+        in the domain.
+
+    Parameters
+    ----------
+    domain_x : `astropy.units.Quantity`
+        1D array of x-coordinates for the plasma domain. Must have
+        units convertable to length.
+    domain_y : `astropy.units.Quantity`
+        1D array of y-coordinates for the plasma domain. Must have
+        units convertable to length.
+    domain_z : `astropy.units.Quantity`
+        1D array of z-coordinates for the plasma domain. Must have
+        units convertable to length.
+
+    """
+    @u.quantity_input(domain_x=u.m, domain_y=u.m, domain_z=u.m)
+    def __init__(self, domain_x, domain_y, domain_z):
+        # Define domain sizes
+        self.x = domain_x
+        self.y = domain_y
+        self.z = domain_z
+
+        self.grid = np.array(np.meshgrid(self.x, self.y, self.z,
+                                         indexing='ij'))
+        self.domain_shape = (len(self.x), len(self.y), len(self.z))
+
+        # Initiate core plasma variables
+        self.density = np.zeros(self.domain_shape) * u.kg / u.m**3
+        self.momentum = np.zeros((3, *self.domain_shape)) * u.kg / (u.m ** 2 * u.s)
+        self.pressure = np.zeros(self.domain_shape) * u.Pa
+        self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T
+        self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m
+
+    @property
+    def velocity(self):
+        return self.momentum / self.density
+
+    @property
+    def magnetic_field_strength(self):
+        B = self.magnetic_field
+        return np.sqrt(np.sum(B * B, axis=0))
+
+    @property
+    def electric_field_strength(self):
+        E = self.electric_field
+        return np.sqrt(np.sum(E * E, axis=0))
+
+    @property
+    def alfven_speed(self):
+        B = self.magnetic_field
+        rho = self.density
+        return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))
+
+    @classmethod
+    def is_datasource_for(cls, **kwargs):
+        if len(kwargs) == 3:
+            match = kwargs.get('domain_x') and \
+                    kwargs.get('domain_y') and \
+                    kwargs.get('domain_z')
+        else:
+            match = False
+        return match
diff --git a/plasmapy/classes/plasma.py b/plasmapy/classes/sources/plasmablob.py
similarity index 63%
rename from plasmapy/classes/plasma.py
rename to plasmapy/classes/sources/plasmablob.py
index 28f279c1..60d3ee0e 100644
--- a/plasmapy/classes/plasma.py
+++ b/plasmapy/classes/sources/plasmablob.py
@@ -24,97 +24,14 @@
 
 from plasmapy.utils import call_string, CouplingWarning
 
+from plasmapy.classes import GenericPlasma
+
 __all__ = [
-    "Plasma3D",
     "PlasmaBlob"
 ]
 
 
-class Plasma3D:
-    """
-    Core class for describing and calculating plasma parameters with
-    spatial dimensions.
-
-    Attributes
-    ----------
-    x : `astropy.units.Quantity`
-        x-coordinates within the plasma domain. Equal to the
-        `domain_x` input parameter.
-    y : `astropy.units.Quantity`
-        y-coordinates within the plasma domain. Equal to the
-        `domain_y` input parameter.
-    z : `astropy.units.Quantity`
-        z-coordinates within the plasma domain. Equal to the
-        `domain_z` input parameter.
-    grid : `astropy.units.Quantity`
-        (3, x, y, z) array containing the values of each coordinate at
-        every point in the domain.
-    domain_shape : tuple
-        Shape of the plasma domain.
-    density : `astropy.units.Quantity`
-        (x, y, z) array of mass density at every point in the domain.
-    momentum : `astropy.units.Quantity`
-        (3, x, y, z) array of the momentum vector at every point in
-        the domain.
-    pressure : `astropy.units.Quantity`
-        (x, y, z) array of pressure at every point in the domain.
-    magnetic_field : `astropy.units.Quantity`
-        (3, x, y, z) array of the magnetic field vector at every point
-        in the domain.
-
-    Parameters
-    ----------
-    domain_x : `astropy.units.Quantity`
-        1D array of x-coordinates for the plasma domain. Must have
-        units convertable to length.
-    domain_y : `astropy.units.Quantity`
-        1D array of y-coordinates for the plasma domain. Must have
-        units convertable to length.
-    domain_z : `astropy.units.Quantity`
-        1D array of z-coordinates for the plasma domain. Must have
-        units convertable to length.
-
-    """
-    @u.quantity_input(domain_x=u.m, domain_y=u.m, domain_z=u.m)
-    def __init__(self, domain_x, domain_y, domain_z):
-        # Define domain sizes
-        self.x = domain_x
-        self.y = domain_y
-        self.z = domain_z
-
-        self.grid = np.array(np.meshgrid(self.x, self.y, self.z,
-                                         indexing='ij'))
-        self.domain_shape = (len(self.x), len(self.y), len(self.z))
-
-        # Initiate core plasma variables
-        self.density = np.zeros(self.domain_shape) * u.kg / u.m**3
-        self.momentum = np.zeros((3, *self.domain_shape)) * u.kg / (u.m ** 2 * u.s)
-        self.pressure = np.zeros(self.domain_shape) * u.Pa
-        self.magnetic_field = np.zeros((3, *self.domain_shape)) * u.T
-        self.electric_field = np.zeros((3, *self.domain_shape)) * u.V / u.m
-
-    @property
-    def velocity(self):
-        return self.momentum / self.density
-
-    @property
-    def magnetic_field_strength(self):
-        B = self.magnetic_field
-        return np.sqrt(np.sum(B * B, axis=0))
-
-    @property
-    def electric_field_strength(self):
-        E = self.electric_field
-        return np.sqrt(np.sum(E * E, axis=0))
-
-    @property
-    def alfven_speed(self):
-        B = self.magnetic_field
-        rho = self.density
-        return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))
-
-
-class PlasmaBlob:
+class PlasmaBlob(GenericPlasma):
     """
     Class for describing and calculating plasma parameters without
     spatial/temporal description.
@@ -245,3 +162,8 @@ def quantum_theta(self):
         """
         theta = quantum_theta(self.T_e, self.n_e)
         return theta
+
+    @classmethod
+    def is_datasource_for(cls, **kwargs):
+        match = kwargs.get('T_e') and kwargs.get('n_e')
+        return match
diff --git a/plasmapy/examples/plot_distribution.py b/plasmapy/examples/plot_distribution.py
index d366d9c2..33bd6639 100644
--- a/plasmapy/examples/plot_distribution.py
+++ b/plasmapy/examples/plot_distribution.py
@@ -29,10 +29,8 @@
 # plasma at a temperature of 30 000 K:
 
 
-p_dens = Maxwellian_1D(v=1 * u.m / u.s,
-                       T=30000 * u.K,
-                       particle='e',
-                       v_drift=0 * u.m / u.s)
+p_dens = Maxwellian_1D(v=1 * u.m / u.s, T=30000 * u.K,
+                       particle='e', V_drift=0 * u.m / u.s)
 print(p_dens)
 
 ############################################################
@@ -60,8 +58,8 @@
 # temperature:
 
 
-std = np.sqrt((Maxwellian_1D(v, T=T, particle='e') * v ** 2 * dv).sum())
-T_theo = (std ** 2 / k_B * m_e).to(u.K)
+std = np.sqrt((Maxwellian_1D(v, T=T, particle='e') * v**2 * dv).sum())
+T_theo = (std**2 / k_B * m_e).to(u.K)
 
 print('T from standard deviation:', T_theo)
 print('Initial T:', T)
diff --git a/plasmapy/examples/plot_physics.py b/plasmapy/examples/plot_physics.py
index bb7b42b4..85cfe746 100644
--- a/plasmapy/examples/plot_physics.py
+++ b/plasmapy/examples/plot_physics.py
@@ -55,6 +55,3 @@
 plt.scatter(
     electron_concentration,
     physics.plasma_frequency(electron_concentration))
-plt.xlabel("Electron Concentration (m^-3)")
-plt.ylabel("Langmuir Wave Plasma Frequency (rad/s)")
-plt.show()
diff --git a/plasmapy/physics/__init__.py b/plasmapy/physics/__init__.py
index ce6cbd9d..010eb5db 100644
--- a/plasmapy/physics/__init__.py
+++ b/plasmapy/physics/__init__.py
@@ -17,10 +17,8 @@
                          cold_plasma_permittivity_SDP)
 
 from .distribution import (Maxwellian_1D,
-                           Maxwellian_velocity_2D,
                            Maxwellian_velocity_3D,
                            Maxwellian_speed_1D,
-                           Maxwellian_speed_2D,
                            Maxwellian_speed_3D,
                            kappa_velocity_3D,
                            kappa_velocity_1D)
diff --git a/plasmapy/physics/dielectric.py b/plasmapy/physics/dielectric.py
index f7b3c8d7..81db6627 100644
--- a/plasmapy/physics/dielectric.py
+++ b/plasmapy/physics/dielectric.py
@@ -1,7 +1,7 @@
 """Functions to calculate plasma dielectric parameters"""
 
 from astropy import units as u
-from plasmapy import utils
+import plasmapy.utils as utils
 from plasmapy.physics import parameters
 from plasmapy.constants import (pi, m_e, c, mu0, e, eps0)
 
diff --git a/plasmapy/physics/distribution.py b/plasmapy/physics/distribution.py
index acb88319..d3176d96 100644
--- a/plasmapy/physics/distribution.py
+++ b/plasmapy/physics/distribution.py
@@ -9,14 +9,14 @@
 import numpy as np
 from scipy.special import gamma
 
-__all__ = ["Maxwellian_1D",
-           "Maxwellian_velocity_2D",
-           "Maxwellian_velocity_3D",
-           "Maxwellian_speed_1D",
-           "Maxwellian_speed_2D",
-           "Maxwellian_speed_3D",
-           "kappa_velocity_1D",
-           "kappa_velocity_3D"]
+__all__ = [
+    "Maxwellian_1D",
+    "Maxwellian_velocity_3D",
+    "Maxwellian_speed_1D",
+    "Maxwellian_speed_3D",
+    "kappa_velocity_1D",
+    "kappa_velocity_3D",
+]
 
 def _v_drift_units(v_drift):
     # Helper method to assign units to  v_drift if it takes a default value
@@ -30,14 +30,13 @@ def _v_drift_units(v_drift):
 def Maxwellian_1D(v,
                   T,
                   particle="e",
-                  v_drift=0,
+                  V_drift=0,
                   vTh=np.nan,
                   units="units"):
     r"""
-    Probability distribution function of velocity for a Maxwellian 
-    distribution in 1D.
+    pdf of velocity for a Maxwellian distribution.
 
-    Returns the probability density function at the velocity `v` in m/s
+    Returns the probability at the velocity `v` in m/s
     to find a particle `particle` in a plasma of temperature `T`
     following the Maxwellian distribution function.
 
@@ -54,7 +53,7 @@ def Maxwellian_1D(v,
         ``'D+'`` for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
         (singly ionized helium-4)), which defaults to electrons.
 
-    v_drift: ~astropy.units.Quantity, optional
+    V_drift: ~astropy.units.Quantity, optional
         The drift velocity in units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
@@ -71,7 +70,7 @@ def Maxwellian_1D(v,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in units of Velocity^-1, normalized so that
+        Probability in units of Velocity^-1, normalized so that
         :math:`\int_{-\infty}^{+\infty} f(v) dv = 1`.
 
     Raises
@@ -104,7 +103,7 @@ def Maxwellian_1D(v,
     >>> from plasmapy.physics import Maxwellian_1D
     >>> from astropy import units as u
     >>> v=1*u.m/u.s
-    >>> Maxwellian_1D(v=v, T=30000 * u.K, particle='e', v_drift=0 * u.m / u.s)
+    >>> Maxwellian_1D(v=v, T= 30000*u.K, particle='e',V_drift=0*u.m/u.s)
     <Quantity 5.91632969e-07 s / m>
     """
 
@@ -114,7 +113,7 @@ def Maxwellian_1D(v,
         v = v.to(u.m / u.s)
         # Catching case where drift velocities have default values, they
         # need to be assigned units
-        v_drift = _v_drift_units(v_drift)
+        V_drift = _v_drift_units(V_drift)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -133,7 +132,7 @@ def Maxwellian_1D(v,
     # Get thermal velocity squared
     vThSq = vTh ** 2
     # Get square of relative particle velocity
-    vSq = (v - v_drift) ** 2
+    vSq = (v - V_drift) ** 2
     # calculating distribution function
     coeff = (vThSq * np.pi) ** (-1 / 2)
     expTerm = np.exp(-vSq / vThSq)
@@ -144,163 +143,23 @@ def Maxwellian_1D(v,
         return distFunc
 
 
-def Maxwellian_velocity_2D(vx,
-                           vy,
-                           T,
-                           particle="e",
-                           vx_drift=0,
-                           vy_drift=0,
-                           vTh=np.nan,
-                           units="units"):
-    r"""
-    Probability distribution function of velocity for a Maxwellian 
-    distribution in 2D.
-
-    Return the probability density function for finding a particle with 
-    velocity components `vx` and `vy` in m/s in an equilibrium plasma of 
-    temperature `T` which follows the 2D Maxwellian distribution function. 
-    This function assumes Cartesian coordinates.
-
-    Parameters
-    ----------
-    vx: ~astropy.units.Quantity
-        The velocity in x-direction units convertible to m/s.
-
-    vy: ~astropy.units.Quantity
-        The velocity in y-direction units convertible to m/s.
-
-    T: ~astropy.units.Quantity
-        The temperature, preferably in Kelvin.
-
-    particle: str, optional
-        Representation of the particle species (e.g., ``'p'`` for protons,
-        ``'D+'`` for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
-        (singly ionized helium-4)), which defaults to electrons.
-
-    vx_drift: ~astropy.units.Quantity, optional
-        The drift velocity in x-direction units convertible to m/s.
-
-    vy_drift: ~astropy.units.Quantity, optional
-        The drift velocity in y-direction units convertible to m/s.
-
-    vTh: ~astropy.units.Quantity, optional
-        Thermal velocity (most probable) in m/s. This is used for
-        optimization purposes to avoid re-calculating `vTh`, for example
-        when integrating over velocity-space.
-
-    units: str, optional
-        Selects whether to run function with units and unit checks (when
-        equal to "units") or to run as unitless (when equal to "unitless").
-        The unitless version is substantially faster for intensive
-        computations.
-
-    Returns
-    -------
-    f : ~astropy.units.Quantity
-        Probability density in Velocity^-1, normalized so that
-        :math:`\iiint_{0}^{\infty} f(\vec{v}) d\vec{v} = 1`.
-
-    Raises
-    ------
-    TypeError
-        A parameter argument is not a `~astropy.units.Quantity` and
-        cannot be converted into a `~astropy.units.Quantity`.
-
-    ~astropy.units.UnitConversionError
-        If the parameters is not in appropriate units.
-
-    ValueError
-        If the temperature is negative, or the particle mass or charge state
-        cannot be found.
-
-    Notes
-    -----
-    In 2D, the Maxwellian velocity distribution function describing
-    the distribution of particles with speed :math:`v` in a plasma with
-    temperature :math:`T` is given by:
-
-    .. math::
-
-        f = (\pi v_{Th}^2)^{-1} \exp \left [-(\vec{v} -
-        \vec{V}_{drift})^2 / v_{Th}^2 \right ]
-
-    where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
-
-    See also
-    --------
-    Maxwellian_1D
-
-    Example
-    -------
-    >>> from plasmapy.physics import Maxwellian_velocity_2D
-    >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
-    >>> Maxwellian_velocity_2D(vx=v,
-    ... vy=v,
-    ... T=30000*u.K,
-    ... particle='e',
-    ... vx_drift=0 * u.m / u.s,
-    ... vy_drift=0 * u.m / u.s)
-    <Quantity 3.5002957e-13 s2 / m2>
-
-
-    """
-    if units == "units":
-        # unit checks and conversions
-        # checking velocity units
-        vx = vx.to(u.m / u.s)
-        vy = vy.to(u.m / u.s)
-        # catching case where drift velocities have default values, they
-        # need to be assigned units
-        vx_drift = _v_drift_units(vx_drift)
-        vy_drift = _v_drift_units(vy_drift)
-        # convert temperature to Kelvins
-        T = T.to(u.K, equivalencies=u.temperature_energy())
-        if np.isnan(vTh):
-            # get thermal velocity and thermal velocity squared
-            vTh = (parameters.thermal_speed(T,
-                                            particle=particle,
-                                            method="most_probable"))
-        elif not np.isnan(vTh):
-            # check units of thermal velocity
-            vTh = vTh.to(u.m / u.s)
-    elif np.isnan(vTh) and units == "unitless":
-        # assuming unitless temperature is in Kelvins
-        vTh = parameters.thermal_speed(T * u.K,
-                                       particle=particle,
-                                       method="most_probable").si.value
-    # accounting for thermal velocity in 2D
-    vThSq = vTh ** 2
-    # Get square of relative particle velocity
-    vSq = ((vx - vx_drift) ** 2 + (vy - vy_drift) ** 2)
-    # calculating distribution function
-    coeff = (vThSq * np.pi) ** (-1)
-    expTerm = np.exp(-vSq / vThSq)
-    distFunc = coeff * expTerm
-    if units == "units":
-        return distFunc.to((u.s / u.m)**2)
-    elif units == "unitless":
-        return distFunc
-
-
 def Maxwellian_velocity_3D(vx,
                            vy,
                            vz,
                            T,
                            particle="e",
-                           vx_drift=0,
-                           vy_drift=0,
-                           vz_drift=0,
+                           Vx_drift=0,
+                           Vy_drift=0,
+                           Vz_drift=0,
                            vTh=np.nan,
                            units="units"):
     r"""
-    Probability distribution function of velocity for a Maxwellian 
-    distribution in 3D.
+    pdf of velocity for a 3D Maxwellian distribution.
 
-    Return the probability density function for finding a particle with 
-    velocity components `vx`, `vy`, and `vz` in m/s in an equilibrium 
-    plasma of temperature `T` which follows the 3D Maxwellian distribution 
-    function. This function assumes Cartesian coordinates.
+    Return the probability of finding a particle with velocity components
+    `vx`, `vy`, and `vz` in m/s in an equilibrium plasma of temperature
+    `T` which follows the 3D Maxwellian distribution function. This
+    function assumes Cartesian coordinates.
 
     Parameters
     ----------
@@ -321,13 +180,13 @@ def Maxwellian_velocity_3D(vx,
         ``'D+'`` for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
         (singly ionized helium-4)), which defaults to electrons.
 
-    vx_drift: ~astropy.units.Quantity, optional
+    Vx_drift: ~astropy.units.Quantity, optional
         The drift velocity in x-direction units convertible to m/s.
 
-    vy_drift: ~astropy.units.Quantity, optional
+    Vy_drift: ~astropy.units.Quantity, optional
         The drift velocity in y-direction units convertible to m/s.
 
-    vz_drift: ~astropy.units.Quantity, optional
+    Vz_drift: ~astropy.units.Quantity, optional
         The drift velocity in z-direction units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
@@ -344,7 +203,7 @@ def Maxwellian_velocity_3D(vx,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in Velocity^-1, normalized so that
+        Probability in Velocity^-1, normalized so that
         :math:`\iiint_{0}^{\infty} f(\vec{v}) d\vec{v} = 1`.
 
     Raises
@@ -362,7 +221,7 @@ def Maxwellian_velocity_3D(vx,
 
     Notes
     -----
-    In 3D, the Maxwellian speed distribution function describing
+    In one dimension, the Maxwellian speed distribution function describing
     the distribution of particles with speed :math:`v` in a plasma with
     temperature :math:`T` is given by:
 
@@ -381,16 +240,16 @@ def Maxwellian_velocity_3D(vx,
     -------
     >>> from plasmapy.physics import Maxwellian_velocity_3D
     >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
+    >>> v=1*u.m/u.s
     >>> Maxwellian_velocity_3D(vx=v,
     ... vy=v,
     ... vz=v,
-    ... T=30000 * u.K,
+    ... T=30000*u.K,
     ... particle='e',
-    ... vx_drift=0 * u.m / u.s,
-    ... vy_drift=0 * u.m / u.s,
-    ... vz_drift=0 * u.m / u.s)
-    <Quantity 2.07089033e-19 s3 / m3>
+    ... Vx_drift=0*u.m/u.s,
+    ... Vy_drift=0*u.m/u.s,
+    ... Vz_drift=0*u.m/u.s)
+    <Quantity 3.98543031e-20 s3 / m3>
 
 
     """
@@ -402,9 +261,9 @@ def Maxwellian_velocity_3D(vx,
         vz = vz.to(u.m / u.s)
         # catching case where drift velocities have default values, they
         # need to be assigned units
-        vx_drift = _v_drift_units(vx_drift)
-        vy_drift = _v_drift_units(vy_drift)
-        vz_drift = _v_drift_units(vz_drift)
+        Vx_drift = _v_drift_units(Vx_drift)
+        Vy_drift = _v_drift_units(Vy_drift)
+        Vz_drift = _v_drift_units(Vz_drift)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -421,9 +280,9 @@ def Maxwellian_velocity_3D(vx,
                                        particle=particle,
                                        method="most_probable").si.value
     # accounting for thermal velocity in 3D
-    vThSq = vTh ** 2
+    vThSq = 3 * vTh ** 2
     # Get square of relative particle velocity
-    vSq = ((vx - vx_drift) ** 2 + (vy - vy_drift) ** 2 + (vz - vz_drift) ** 2)
+    vSq = ((vx - Vx_drift) ** 2 + (vy - Vy_drift) ** 2 + (vz - Vz_drift) ** 2)
     # calculating distribution function
     coeff = (vThSq * np.pi) ** (-3 / 2)
     expTerm = np.exp(-vSq / vThSq)
@@ -437,16 +296,15 @@ def Maxwellian_velocity_3D(vx,
 def Maxwellian_speed_1D(v,
                         T,
                         particle="e",
-                        v_drift=0,
+                        V_drift=0,
                         vTh=np.nan,
                         units="units"):
     r"""
-    Probability distribution function of speed for a Maxwellian distribution
-    in 1D.
+    pdf of speed for a Maxwellian distribution.
 
-    Return the probability density function for finding a particle with 
-    speed `v` in m/s in an equilibrium plasma of temperature `T` which 
-    follows the Maxwellian distribution function.
+    Return the probability of finding a particle with speed `v` in m/s
+    in an equilibrium plasma of temperature `T` which follows the
+    Maxwellian distribution function.
 
     Parameters
     ----------
@@ -461,7 +319,7 @@ def Maxwellian_speed_1D(v,
         for deuterium, or `'He-4 +1'` for :math:`He_4^{+1}`
         (singly ionized helium-4)), which defaults to electrons.
 
-    v_drift: ~astropy.units.Quantity
+    V_drift: ~astropy.units.Quantity
         The drift speed in units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
@@ -478,7 +336,7 @@ def Maxwellian_speed_1D(v,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in speed^-1, normalized so that
+        Probability in speed^-1, normalized so that
         :math:`\int_{0}^{\infty} f(v) dv = 1`.
 
     Raises
@@ -502,7 +360,7 @@ def Maxwellian_speed_1D(v,
 
     .. math::
 
-       f(v) = 2 \frac{1}{(\pi v_{Th}^2)^{1/2}} \exp(-(v - V_{drift})^2 / v_{Th}^2 )
+       f(v) = 4\pi v^2 \frac{1}{(\pi v_{Th}^2)^{3/2}} \exp(-(v - V_{drift})^2 / v_{Th}^2 )
 
     where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
 
@@ -510,143 +368,18 @@ def Maxwellian_speed_1D(v,
     -------
     >>> from plasmapy.physics import Maxwellian_speed_1D
     >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
-    >>> Maxwellian_speed_1D(v=v, T=30000 * u.K, particle='e', v_drift=0 * u.m / u.s)
-    <Quantity 1.18326594e-06 s / m>
-
-    """
-    if units == "units":
-        # unit checks and conversions
-        # checking velocity units
-        v = v.to(u.m / u.s)
-        # Catching case where drift velocities have default values, they
-        # need to be assigned units
-        v_drift = _v_drift_units(v_drift)
-        # convert temperature to Kelvins
-        T = T.to(u.K, equivalencies=u.temperature_energy())
-        if np.isnan(vTh):
-            # get thermal velocity and thermal velocity squared
-            vTh = (parameters.thermal_speed(T,
-                                            particle=particle,
-                                            method="most_probable"))
-        elif not np.isnan(vTh):
-            # check units of thermal velocity
-            vTh = vTh.to(u.m / u.s)
-    elif np.isnan(vTh) and units == "unitless":
-        # assuming unitless temperature is in Kelvins
-        vTh = (parameters.thermal_speed(T * u.K,
-                                        particle=particle,
-                                        method="most_probable")).si.value
-    # Get thermal velocity squared
-    vThSq = vTh ** 2
-    # Get square of relative particle velocity
-    vSq = (v - v_drift) ** 2
-    # calculating distribution function
-    coeff = 2 * (vThSq * np.pi) ** (-1 / 2)
-    expTerm = np.exp(-vSq / vThSq)
-    distFunc = coeff * expTerm
-    if units == "units":
-        return distFunc.to(u.s / u.m)
-    elif units == "unitless":
-        return distFunc
-
-
-def Maxwellian_speed_2D(v,
-                        T,
-                        particle="e",
-                        v_drift=0,
-                        vTh=np.nan,
-                        units="units"):
-    r"""
-    Probability distribution function of speed for a Maxwellian distribution
-    in 2D.
-
-    Return the probability density function of finding a particle with speed components
-    `vx` and `vy` in m/s in an equilibrium plasma of temperature
-    `T` which follows the 2D Maxwellian distribution function. This
-    function assumes Cartesian coordinates.
-
-    Parameters
-    ----------
-    v: ~astropy.units.Quantity
-        The speed in units convertible to m/s.
-
-    T: ~astropy.units.Quantity
-        The temperature, preferably in Kelvin.
-
-    particle: str, optional
-        Representation of the particle species(e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for :math:`He_4^{+1}`
-        (singly ionized helium-4)), which defaults to electrons.
-
-    v_drift: ~astropy.units.Quantity
-        The drift speed in units convertible to m/s.
-
-    vTh: ~astropy.units.Quantity, optional
-        Thermal velocity (most probable) in m/s. This is used for
-        optimization purposes to avoid re-calculating vTh, for example
-        when integrating over velocity-space.
-
-    units: str, optional
-        Selects whether to run function with units and unit checks (when
-        equal to "units") or to run as unitless (when equal to "unitless").
-        The unitless version is substantially faster for intensive
-        computations.
-
-    Returns
-    -------
-    f : ~astropy.units.Quantity
-        Probability density in speed^-1, normalized so that:
-        :math:`\iiint_{0}^{\infty} f(\vec{v}) d\vec{v} = 1`.
-
-    Raises
-    ------
-    TypeError
-        A parameter argument is not a `~astropy.units.Quantity` and
-        cannot be converted into a `~astropy.units.Quantity`.
-
-    ~astropy.units.UnitConversionError
-        If the parameters is not in appropriate units.
-
-    ValueError
-        If the temperature is negative, or the particle mass or charge state
-        cannot be found.
-
-    Notes
-    -----
-    In 2D, the Maxwellian speed distribution function describing
-    the distribution of particles with speed :math:`v` in a plasma with
-    temperature :math:`T` is given by:
-
-    .. math::
-
-       f = 2 \pi \vec{v} (\pi v_{Th}^2)^{-1} \exp(-(\vec{v} -
-       \vec{V}_{drift})^2 / v_{Th}^2)
-
-    where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
-
-    See also
-    --------
-    Maxwellian_speed_1D
-
-    Example
-    -------
-    >>> from plasmapy.physics import Maxwellian_speed_2D
-    >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
-    >>> Maxwellian_speed_2D(v=v, T=30000 * u.K, particle='e', v_drift=0 * u.m / u.s)
-    <Quantity 2.19930065e-12 s / m>
+    >>> v=1*u.m/u.s
+    >>> Maxwellian_speed_1D(v=v, T= 30000*u.K, particle='e',V_drift=0*u.m/u.s)
+    <Quantity 2.60235754e-18 s / m>
 
     """
-    if v_drift != 0:
-        raise NotImplementedError("Non-zero drift speed is work in progress.")
     if units == "units":
         # unit checks and conversions
         # checking velocity units
         v = v.to(u.m / u.s)
         # Catching case where drift velocity has default value, and
         # needs to be assigned units
-        v_drift = _v_drift_units(v_drift)
+        V_drift = _v_drift_units(V_drift)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -665,10 +398,10 @@ def Maxwellian_speed_2D(v,
     # getting square of thermal speed
     vThSq = vTh ** 2
     # get square of relative particle speed
-    vSq = (v - v_drift) ** 2
+    vSq = (v - V_drift) ** 2
     # calculating distribution function
-    coeff1 = (np.pi * vThSq) ** (-1)
-    coeff2 = 2 * np.pi * (v - v_drift)
+    coeff1 = (np.pi * vThSq) ** (-3 / 2)
+    coeff2 = 4 * np.pi * vSq
     expTerm = np.exp(-vSq / vThSq)
     distFunc = coeff1 * coeff2 * expTerm
     if units == "units":
@@ -677,36 +410,51 @@ def Maxwellian_speed_2D(v,
         return distFunc
 
 
-def Maxwellian_speed_3D(v,
+def Maxwellian_speed_3D(vx,
+                        vy,
+                        vz,
                         T,
                         particle="e",
-                        v_drift=0,
+                        Vx_drift=0,
+                        Vy_drift=0,
+                        Vz_drift=0,
                         vTh=np.nan,
                         units="units"):
     r"""
-    Probability distribution function of speed for a Maxwellian
-    distribution in 3D.
+    pdf of speed for a 3D Maxwellian distribution.
 
-    Return the probability density function for finding a particle with 
-    speed components `vx`, `vy`, and `vz` in m/s in an equilibrium 
-    plasma of temperature `T` which follows the 3D Maxwellian 
-    distribution function. This function assumes Cartesian coordinates.
+    Return the probability of finding a particle with speed components
+    `vx`, `vy`, and `vz` in m/s in an equilibrium plasma of temperature
+    `T` which follows the 3D Maxwellian distribution function. This
+    function assumes Cartesian coordinates.
 
     Parameters
     ----------
-    v: ~astropy.units.Quantity
-        The speed in units convertible to m/s.
+    vx: ~astropy.units.Quantity
+        The speed in x-direction units convertible to m/s.
+
+    vy: ~astropy.units.Quantity
+        The speed in y-direction units convertible to m/s.
+
+    vz: ~astropy.units.Quantity
+        The speed in z-direction units convertible to m/s.
 
     T: ~astropy.units.Quantity
         The temperature, preferably in Kelvin.
 
     particle: str, optional
-        Representation of the particle species(e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for :math:`He_4^{+1}`
+        Representation of the particle species(e.g., 'p' for protons, 'D+'
+        for deuterium, or 'He-4 +1' for :math:`He_4^{+1}`
         (singly ionized helium-4)), which defaults to electrons.
 
-    v_drift: ~astropy.units.Quantity
-        The drift speed in units convertible to m/s.
+    Vx_drift: ~astropy.units.Quantity
+        The drift speed in x-direction units convertible to m/s.
+
+    Vy_drift: ~astropy.units.Quantity
+        The drift speed in y-direction units convertible to m/s.
+
+    Vz_drift: ~astropy.units.Quantity
+        The drift speed in z-direction units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
         Thermal velocity (most probable) in m/s. This is used for
@@ -722,7 +470,7 @@ def Maxwellian_speed_3D(v,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in speed^-1, normalized so that:
+        Probability in speed^-1, normalized so that:
         :math:`\iiint_{0}^{\infty} f(\vec{v}) d\vec{v} = 1`.
 
     Raises
@@ -740,7 +488,7 @@ def Maxwellian_speed_3D(v,
 
     Notes
     -----
-    In 3D, the Maxwellian speed distribution function describing
+    In one dimension, the Maxwellian speed distribution function describing
     the distribution of particles with speed :math:`v` in a plasma with
     temperature :math:`T` is given by:
 
@@ -759,20 +507,29 @@ def Maxwellian_speed_3D(v,
     -------
     >>> from plasmapy.physics import Maxwellian_speed_3D
     >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
-    >>> Maxwellian_speed_3D(v=v, T=30000*u.K, particle='e', v_drift=0 * u.m / u.s)
-    <Quantity 2.60235754e-18 s / m>
+    >>> v=1*u.m/u.s
+    >>> Maxwellian_speed_3D(vx=v,
+    ... vy=v,
+    ... vz=v,
+    ... T=30000*u.K,
+    ... particle='e',
+    ... Vx_drift=0*u.m/u.s,
+    ... Vy_drift=0*u.m/u.s,
+    ... Vz_drift=0*u.m/u.s)
+    <Quantity 1.76238544e-53 s3 / m3>
 
     """
-    if v_drift != 0:
-        raise NotImplementedError("Non-zero drift speed is work in progress.")
     if units == "units":
         # unit checks and conversions
         # checking velocity units
-        v = v.to(u.m / u.s)
-        # Catching case where drift velocity has default value, and
-        # needs to be assigned units
-        v_drift = _v_drift_units(v_drift)
+        vx = vx.to(u.m / u.s)
+        vy = vy.to(u.m / u.s)
+        vz = vz.to(u.m / u.s)
+        # Catching case where drift velocities have default values, they
+        # need to be assigned units
+        Vx_drift = _v_drift_units(Vx_drift)
+        Vy_drift = _v_drift_units(Vy_drift)
+        Vz_drift = _v_drift_units(Vz_drift)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -788,17 +545,29 @@ def Maxwellian_speed_3D(v,
         vTh = (parameters.thermal_speed(T * u.K,
                                         particle=particle,
                                         method="most_probable")).si.value
-    # getting square of thermal speed
-    vThSq = vTh ** 2
-    # get square of relative particle speed
-    vSq = (v - v_drift) ** 2
-    # calculating distribution function
-    coeff1 = (np.pi * vThSq) ** (-3 / 2)
-    coeff2 = 4 * np.pi * vSq
-    expTerm = np.exp(-vSq / vThSq)
-    distFunc = coeff1 * coeff2 * expTerm
+    # getting distribution functions along each axis
+    fx = Maxwellian_speed_1D(vx,
+                             T,
+                             particle=particle,
+                             V_drift=Vx_drift,
+                             vTh=vTh,
+                             units=units)
+    fy = Maxwellian_speed_1D(vy,
+                             T,
+                             particle=particle,
+                             V_drift=Vy_drift,
+                             vTh=vTh,
+                             units=units)
+    fz = Maxwellian_speed_1D(vz,
+                             T,
+                             particle=particle,
+                             V_drift=Vz_drift,
+                             vTh=vTh,
+                             units=units)
+    # multiplying probabilities in each axis to get 3D probability
+    distFunc = fx * fy * fz
     if units == "units":
-        return distFunc.to(u.s / u.m)
+        return distFunc.to((u.s / u.m)**3)
     elif units == "unitless":
         return distFunc
 
@@ -807,13 +576,13 @@ def kappa_velocity_1D(v,
                       T,
                       kappa,
                       particle="e",
-                      v_drift=0,
+                      V_drift=0,
                       vTh=np.nan,
                       units="units"):
     r"""
-    Return the probability density at the velocity `v` in m/s
+    Return the probability at the velocity `v` in m/s
     to find a particle `particle` in a plasma of temperature `T`
-    following the Kappa distribution function in 1D. The slope of the
+    following the Kappa distribution function. The slope of the
     tail of the Kappa distribution function is set by 'kappa', which
     must be greater than :math:`1/2`.
 
@@ -836,7 +605,7 @@ def kappa_velocity_1D(v,
         for deuterium, or `'He-4 +1'` for :math:`He_4^{+1}`
         (singly ionized helium-4)), which defaults to electrons.
 
-    v_drift: ~astropy.units.Quantity, optional
+    V_drift: ~astropy.units.Quantity, optional
         The drift velocity in units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
@@ -853,7 +622,7 @@ def kappa_velocity_1D(v,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in Velocity^-1, normalized so that
+        probability in Velocity^-1, normalized so that
         :math:`\int_{-\infty}^{+\infty} f(v) dv = 1`.
 
     Raises
@@ -893,8 +662,8 @@ def kappa_velocity_1D(v,
     --------
     >>> from plasmapy.physics import kappa_velocity_1D
     >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
-    >>> kappa_velocity_1D(v=v, T=30000*u.K, kappa=4, particle='e', v_drift=0 * u.m / u.s)
+    >>> v=1*u.m/u.s
+    >>> kappa_velocity_1D(v=v, T=30000*u.K, kappa=4, particle='e',V_drift=0*u.m/u.s)
     <Quantity 6.75549854e-07 s / m>
 
     See Also
@@ -911,11 +680,11 @@ def kappa_velocity_1D(v,
         v = v.to(u.m / u.s)
         # catching case where drift velocities have default values, they
         # need to be assigned units
-        if v_drift == 0:
-            if not isinstance(v_drift, astropy.units.quantity.Quantity):
-                v_drift = v_drift * u.m / u.s
+        if V_drift == 0:
+            if not isinstance(V_drift, astropy.units.quantity.Quantity):
+                V_drift = V_drift * u.m / u.s
         # checking units of drift velocities
-        v_drift = v_drift.to(u.m / u.s)
+        V_drift = V_drift.to(u.m / u.s)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -934,7 +703,7 @@ def kappa_velocity_1D(v,
     # Get thermal velocity squared and accounting for 1D instead of 3D
     vThSq = vTh ** 2
     # Get square of relative particle velocity
-    vSq = (v - v_drift) ** 2
+    vSq = (v - V_drift) ** 2
     # calculating distribution function
     expTerm = (1 + vSq / (kappa * vThSq)) ** (-kappa)
     coeff1 = 1 / (np.sqrt(np.pi) * kappa ** (3 / 2) * vTh)
@@ -952,17 +721,16 @@ def kappa_velocity_3D(vx,
                       T,
                       kappa,
                       particle="e",
-                      vx_drift=0,
-                      vy_drift=0,
-                      vz_drift=0,
+                      Vx_drift=0,
+                      Vy_drift=0,
+                      Vz_drift=0,
                       vTh=np.nan,
                       units="units"):
     r"""
-    Return the probability density function for finding a particle with 
-    velocity components `v_x`, `v_y`, and `v_z`in m/s in a suprathermal 
-    plasma of temperature `T` and parameter 'kappa' which follows the 
-    3D Kappa distribution function. This function assumes Cartesian 
-    coordinates.
+    Return the probability of finding a particle with velocity components
+    `v_x`, `v_y`, and `v_z`in m/s in a suprathermal plasma of temperature
+    `T` and parameter 'kappa' which follows the 3D Kappa distribution function.
+    This function assumes Cartesian coordinates.
 
     Parameters
     ----------
@@ -989,13 +757,13 @@ def kappa_velocity_3D(vx,
         for deuterium, or 'He-4 +1' for :math:`He_4^{+1}` : singly ionized
         helium-4)), which defaults to electrons.
 
-    vx_drift: ~astropy.units.Quantity, optional
+    Vx_drift: ~astropy.units.Quantity, optional
         The drift velocity in x-direction units convertible to m/s.
 
-    vy_drift: ~astropy.units.Quantity, optional
+    Vy_drift: ~astropy.units.Quantity, optional
         The drift velocity in y-direction units convertible to m/s.
 
-    vz_drift: ~astropy.units.Quantity, optional
+    Vz_drift: ~astropy.units.Quantity, optional
         The drift velocity in z-direction units convertible to m/s.
 
     vTh: ~astropy.units.Quantity, optional
@@ -1012,7 +780,7 @@ def kappa_velocity_3D(vx,
     Returns
     -------
     f : ~astropy.units.Quantity
-        Probability density in Velocity^-1, normalized so that:
+        probability in Velocity^-1, normalized so that:
         :math:`\iiint_{0}^{\infty} f(\vec{v}) d\vec{v} = 1`
 
     Raises
@@ -1056,16 +824,16 @@ def kappa_velocity_3D(vx,
     -------
     >>> from plasmapy.physics import kappa_velocity_3D
     >>> from astropy import units as u
-    >>> v=1 * u.m / u.s
+    >>> v=1*u.m/u.s
     >>> kappa_velocity_3D(vx=v,
     ... vy=v,
     ... vz=v,
-    ... T=30000 * u.K,
+    ... T=30000*u.K,
     ... kappa=4,
     ... particle='e',
-    ... vx_drift=0 * u.m / u.s,
-    ... vy_drift=0 * u.m / u.s,
-    ... vz_drift=0 * u.m / u.s)
+    ... Vx_drift=0*u.m/u.s,
+    ... Vy_drift=0*u.m/u.s,
+    ... Vz_drift=0*u.m/u.s)
     <Quantity 3.7833988e-19 s3 / m3>
     """
     # must have kappa > 3/2 for distribution function to be valid
@@ -1079,9 +847,9 @@ def kappa_velocity_3D(vx,
         vz = vz.to(u.m / u.s)
         # Catching case where drift velocities have default values, they
         # need to be assigned units
-        vx_drift = _v_drift_units(vx_drift)
-        vy_drift = _v_drift_units(vy_drift)
-        vz_drift = _v_drift_units(vz_drift)
+        Vx_drift = _v_drift_units(Vx_drift)
+        Vy_drift = _v_drift_units(Vy_drift)
+        Vz_drift = _v_drift_units(Vz_drift)
         # convert temperature to Kelvins
         T = T.to(u.K, equivalencies=u.temperature_energy())
         if np.isnan(vTh):
@@ -1100,7 +868,7 @@ def kappa_velocity_3D(vx,
     # getting square of thermal velocity
     vThSq = vTh ** 2
     # Get square of relative particle velocity
-    vSq = ((vx - vx_drift) ** 2 + (vy - vy_drift) ** 2 + (vz - vz_drift) ** 2)
+    vSq = ((vx - Vx_drift) ** 2 + (vy - Vy_drift) ** 2 + (vz - Vz_drift) ** 2)
     # calculating distribution function
     expTerm = (1 + vSq / (kappa * vThSq)) ** (-(kappa + 1))
     coeff1 = 1 / (2 * np.pi * (kappa * vThSq) ** (3 / 2))
diff --git a/plasmapy/utils/checks.py b/plasmapy/utils/checks.py
index 703c7548..8e6e6245 100644
--- a/plasmapy/utils/checks.py
+++ b/plasmapy/utils/checks.py
@@ -409,15 +409,13 @@ def _check_relativistic(V, funcname, betafrac=0.05):
     if not isinstance(V, u.Quantity):
         raise TypeError(errmsg)
 
-    try:
-        V_over_c = (V / c).to_value(u.dimensionless_unscaled)
-    except Exception:
+    if V.si.unit != u.m / u.s:
         raise u.UnitConversionError(errmsg)
 
     if np.any(np.isnan(V.value)):
         raise ValueError(f"V includes NaNs in {funcname}")
 
-    beta = np.max(np.abs((V_over_c)))
+    beta = np.max(np.abs((V.si / c).value))
 
     if beta == np.inf:
         raise RelativityError(f"{funcname} is yielding an infinite velocity.")
diff --git a/plasmapy/utils/datatype_factory_base.py b/plasmapy/utils/datatype_factory_base.py
new file mode 100644
index 00000000..7e2a1946
--- /dev/null
+++ b/plasmapy/utils/datatype_factory_base.py
@@ -0,0 +1,192 @@
+# SunPy is released under a BSD-style open source licence:
+
+# Copyright (c) 2013-2018 The SunPy developers
+# All rights reserved.
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+
+# * Redistributions of source code must retain the above copyright
+#   notice, this list of conditions and the following disclaimer.
+
+# * Redistributions in binary form must reproduce the above copyright
+#   notice, this list of conditions and the following disclaimer in the
+#   documentation and/or other materials provided with the distribution.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+# A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT
+# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+# Project     : https://github.com/sunpy/sunpy
+# File        : sunpy/util/datatype_factory_base.py
+# Commit hash : f6330eea602ea796b5b004dee283b8877b24da23
+
+
+import inspect
+
+
+class BasicRegistrationFactory:
+    """
+    Generalized registerable factory type.
+
+    Widgets (classes) can be registered with an instance of this class.
+    Arguments to the factory's `__call__` method are then passed to a function
+    specified by the registered factory, which validates the input and returns
+    a instance of the class that best matches the inputs.
+
+    Attributes
+    ----------
+
+    registry : dict
+        Dictionary mapping classes (key) to function (value) which validates
+        input.
+
+    default_widget_type : type
+        Class of the default widget.  Defaults to None.
+
+    validation_functions : list of strings
+        List of function names that are valid validation functions.
+
+    Parameters
+    ----------
+
+    default_widget_type : type, optional
+
+    additional_validation_functions : list of strings, optional
+        List of strings corresponding to additional validation function names.
+
+    Notes
+    -----
+
+    * A valid validation function must be a classmethod of the registered widget
+      and it must return True or False.
+
+    """
+
+    def __init__(self, default_widget_type=None,
+                 additional_validation_functions=[], registry=None):
+
+        if registry is None:
+            self.registry = dict()
+        else:
+            self.registry = registry
+
+        self.default_widget_type = default_widget_type
+
+        self.validation_functions = (['_factory_validation_function'] +
+                                     additional_validation_functions)
+
+    def __call__(self, *args, **kwargs):
+        """ Method for running the factory.
+
+        Arguments args and kwargs are passed through to the validation
+        function and to the constructor for the final type.
+        """
+
+        # Any preprocessing and massaging of inputs can happen here
+
+        return self._check_registered_widget(*args, **kwargs)
+
+    def _check_registered_widget(self, *args, **kwargs):
+
+        candidate_widget_types = list()
+
+        for key in self.registry:
+
+            # Call the registered validation function for each registered class
+            if self.registry[key](*args, **kwargs):
+                candidate_widget_types.append(key)
+
+        n_matches = len(candidate_widget_types)
+
+        if n_matches == 0:
+            if self.default_widget_type is None:
+                raise NoMatchError("No types match specified arguments and no default is set.")
+            else:
+                candidate_widget_types = [self.default_widget_type]
+        elif n_matches > 1:
+            print(candidate_widget_types)
+            raise MultipleMatchError("Too many candidate types identified ({0})."
+                                     "Specify enough keywords to guarantee unique type "
+                                     "identification.".format(n_matches))
+
+        # Only one is found
+        WidgetType = candidate_widget_types[0]
+
+        return WidgetType(*args, **kwargs)
+
+    def register(self, WidgetType, validation_function=None, is_default=False):
+        """ Register a widget with the factory.
+
+        If `validation_function` is not specified, tests `WidgetType` for
+        existence of any function in in the list `self.validation_functions`,
+        which is a list of strings which must be callable class attribute
+
+        Parameters
+        ----------
+
+        WidgetType : type
+            Widget to register.
+
+        validation_function : function, optional
+            Function to validate against.  Defaults to None, which indicates
+            that a classmethod in validation_functions is used.
+
+        is_default : bool, optional
+            Sets WidgetType to be the default widget.
+
+        """
+        if is_default:
+            self.default_widget_type = WidgetType
+
+        elif validation_function is not None:
+            if not callable(validation_function):
+                raise AttributeError("Keyword argument 'validation_function' must be callable.")
+
+            self.registry[WidgetType] = validation_function
+
+        else:
+            found = False
+            for vfunc_str in self.validation_functions:
+                if hasattr(WidgetType, vfunc_str):
+                    vfunc = getattr(WidgetType, vfunc_str)
+
+                    # check if classmethod: stackoverflow #19227724
+                    _classmethod = inspect.ismethod(vfunc) and vfunc.__self__ is WidgetType
+
+                    if _classmethod:
+                        self.registry[WidgetType] = vfunc
+                        found = True
+                        break
+                    else:
+                        raise ValidationFunctionError("{0}.{1} must be a classmethod."
+                                                      .format(WidgetType.__name__, vfunc_str))
+
+            if not found:
+                raise ValidationFunctionError("No proper validation function for class {0} "
+                                              "found.".format(WidgetType.__name__))
+
+    def unregister(self, WidgetType):
+        """ Remove a widget from the factory's registry."""
+        self.registry.pop(WidgetType)
+
+
+class NoMatchError(Exception):
+    """Exception for when no candidate class is found."""
+
+
+class MultipleMatchError(Exception):
+    """Exception for when too many candidate classes are found."""
+
+
+class ValidationFunctionError(AttributeError):
+    """Exception for when no candidate class is found."""


Test patch:
diff --git a/plasmapy/classes/tests/test_Plasma.py b/plasmapy/classes/tests/test_Plasma.py
index 206db539..e0e901db 100644
--- a/plasmapy/classes/tests/test_Plasma.py
+++ b/plasmapy/classes/tests/test_Plasma.py
@@ -2,7 +2,7 @@
 import numpy as np
 import astropy.units as u
 
-from plasmapy.classes import plasma
+from plasmapy.classes.sources import plasma3d, plasmablob
 from plasmapy.utils.exceptions import InvalidParticleError
 
 @pytest.mark.parametrize('grid_dimensions, expected_size', [
@@ -33,9 +33,9 @@ def test_Plasma3D_setup(grid_dimensions, expected_size):
     >>> test_Plasma3D_setup((100, 10, 1), 1000)
     """
     x, y, z = grid_dimensions
-    test_plasma = plasma.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)
+    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
@@ -68,9 +68,9 @@ def test_Plasma3D_derived_vars():
     variables.  The core variables are set with arbitrary uniform
     values.
     """
-    test_plasma = plasma.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)
+    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
@@ -113,10 +113,10 @@ def test_intermediate_coupling(self):
         n_e = 1e26 * u.cm ** -3
         Z = 2.0
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Intermediate coupling regime: Gamma = 10.585076050938532.'
         regime, _ = blob.regimes()
@@ -140,10 +140,10 @@ def test_strongly_coupled(self):
         n_e = 1e26 * u.cm ** -3
         Z = 3.0
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Strongly coupled regime: Gamma = 104.02780112828943.'
 
@@ -168,10 +168,10 @@ def test_weakly_coupled(self):
         n_e = 1e15 * u.cm ** -3
         Z = 2.5
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Weakly coupled regime: Gamma = 0.0075178096952688445.'
 
@@ -196,10 +196,10 @@ def test_thermal_kinetic_energy_dominant(self):
         n_e = 1e20 * u.cm ** -3
         Z = 2.5
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Thermal kinetic energy dominant: Theta = 120.65958493847927'
 
@@ -224,10 +224,10 @@ def test_fermi_quantum_energy_dominant(self):
         n_e = 1e26 * u.cm ** -3
         Z = 3.0
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Fermi quantum energy dominant: Theta = 0.009872147858602853'
 
@@ -252,10 +252,10 @@ def test_both_fermi_and_thermal_energy_important(self):
         n_e = 1e25 * u.cm ** -3
         Z = 2.0
         particle = 'p'
-        blob = plasma.PlasmaBlob(T_e=T_e,
-                                 n_e=n_e,
-                                 Z=Z,
-                                 particle=particle)
+        blob = plasmablob.PlasmaBlob(T_e=T_e,
+                                     n_e=n_e,
+                                     Z=Z,
+                                     particle=particle)
 
         expect_regime = 'Both Fermi and thermal energy important: Theta = 0.03818537605355442'
 
@@ -274,10 +274,10 @@ def setup_class(self):
         self.n_e = 1e23 * u.cm ** -3
         self.Z = 2.5
         self.particle = 'p'
-        self.blob = plasma.PlasmaBlob(T_e=self.T_e,
-                                      n_e=self.n_e,
-                                      Z=self.Z,
-                                      particle=self.particle)
+        self.blob = plasmablob.PlasmaBlob(T_e=self.T_e,
+                                          n_e=self.n_e,
+                                          Z=self.Z,
+                                          particle=self.particle)
         self.couplingVal = 10.468374460435724
         self.thetaVal = 0.6032979246923964
 
@@ -286,10 +286,10 @@ def test_invalid_particle(self):
         Checks if function raises error for invalid particle.
         """
         with pytest.raises(InvalidParticleError):
-            plasma.PlasmaBlob(T_e=self.T_e,
-                              n_e=self.n_e,
-                              Z=self.Z,
-                              particle="cupcakes")
+            plasmablob.PlasmaBlob(T_e=self.T_e,
+                                  n_e=self.n_e,
+                                  Z=self.Z,
+                                  particle="cupcakes")
 
     def test_electron_temperature(self):
         """Testing if we get the same electron temperature we put in """
diff --git a/plasmapy/classes/tests/test_Species.py b/plasmapy/classes/tests/test_Species.py
index 0bdff05b..bbcc78f2 100644
--- a/plasmapy/classes/tests/test_Species.py
+++ b/plasmapy/classes/tests/test_Species.py
@@ -4,7 +4,8 @@
 from astropy.modeling import models, fitting
 from scipy.optimize import curve_fit
 
-from plasmapy.classes import Plasma3D, Species
+from plasmapy.classes import Species
+from plasmapy.classes.sources import Plasma3D
 
 
 @pytest.fixture()
diff --git a/plasmapy/classes/tests/test_plasma_base.py b/plasmapy/classes/tests/test_plasma_base.py
new file mode 100644
index 00000000..eb6be2d2
--- /dev/null
+++ b/plasmapy/classes/tests/test_plasma_base.py
@@ -0,0 +1,52 @@
+from plasmapy.classes import BasePlasma, GenericPlasma
+# Get rid of any previously registered classes.
+BasePlasma._registry = {}
+
+
+class NoDataSource(BasePlasma):
+    pass
+
+
+class IsDataSource(BasePlasma):
+    @classmethod
+    def is_datasource_for(cls, **kwargs):
+        return True
+
+
+class IsNotDataSource(BasePlasma):
+    @classmethod
+    def is_datasource_for(cls, **kwargs):
+        return False
+
+
+class TestRegistrar:
+    def test_no_data_source(self):
+        """
+        NoDataSource class should not be registered since it has
+        no method named ``is_datasource_for``.
+        """
+        assert not BasePlasma._registry.get(NoDataSource)
+
+    def test_is_data_source(self):
+        """
+        IsDataSource class should be registered since it has a
+        method named ``is_datasource_for`` and must return True.
+        """
+        assert BasePlasma._registry.get(IsDataSource)
+        assert BasePlasma._registry[IsDataSource]()
+        # Delete the class from _registry once test is done
+        # to not interfere with plasma factory tests
+        del BasePlasma._registry[IsDataSource]
+
+    def test_is_not_data_source(self):
+        """
+        IsNotDataSource class should be registered since it has a
+        method named ``is_datasource_for`` but must return False.
+        """
+        assert BasePlasma._registry.get(IsNotDataSource)
+        assert not BasePlasma._registry[IsNotDataSource]()
+        del BasePlasma._registry[IsNotDataSource]
+
+
+def test_subclasses():
+    assert issubclass(GenericPlasma, BasePlasma)
diff --git a/plasmapy/classes/tests/test_plasma_factory.py b/plasmapy/classes/tests/test_plasma_factory.py
new file mode 100644
index 00000000..1957d0df
--- /dev/null
+++ b/plasmapy/classes/tests/test_plasma_factory.py
@@ -0,0 +1,32 @@
+import astropy.units as u
+import numpy as np
+
+import plasmapy.classes
+
+
+class TestPlasma(object):
+    def test_patters(self):
+        # Input data whose specific subclass cannot be known
+        generic = plasmapy.classes.Plasma(blablobleh='spam')
+        assert isinstance(generic, plasmapy.classes.GenericPlasma)
+
+    def test_Plasma3D(self):
+        # Input data for Plasma3D
+        three_dimensional = plasmapy.classes.Plasma(domain_x=np.linspace(0, 1, 3) * u.m,
+                                                    domain_y=np.linspace(0, 1, 3) * u.m,
+                                                    domain_z=np.linspace(0, 1, 3) * u.m)
+        assert isinstance(three_dimensional, plasmapy.classes.sources.Plasma3D)
+        assert isinstance(three_dimensional, plasmapy.classes.BasePlasma)
+
+    def test_PlasmaBlob(self):
+        # Input data for PlasmaBlob
+        T_e = 25 * 15e3 * u.K
+        n_e = 1e26 * u.cm ** -3
+        Z = 2.0
+        particle = 'p'
+        blob = plasmapy.classes.Plasma(T_e=T_e,
+                                       n_e=n_e,
+                                       Z=Z,
+                                       particle=particle)
+        assert isinstance(blob, plasmapy.classes.sources.PlasmaBlob)
+        assert isinstance(blob, plasmapy.classes.BasePlasma)
diff --git a/plasmapy/physics/tests/test_distribution.py b/plasmapy/physics/tests/test_distribution.py
index f99c5113..35e7ffb0 100644
--- a/plasmapy/physics/tests/test_distribution.py
+++ b/plasmapy/physics/tests/test_distribution.py
@@ -8,8 +8,6 @@
 from ...constants import (m_p, m_e, c, mu0, k_B, e, eps0, pi)
 from ..distribution import (Maxwellian_1D,
                             Maxwellian_speed_1D,
-                            Maxwellian_velocity_2D,
-                            Maxwellian_speed_2D,
                             Maxwellian_velocity_3D,
                             Maxwellian_speed_3D,
                             kappa_velocity_1D,
@@ -26,9 +24,9 @@ def setup_class(self):
         """initializing parameters for tests """
         self.T_e = 30000 * u.K
         self.v = 1e5 * u.m / u.s
-        self.v_drift = 1000000 * u.m / u.s
-        self.v_drift2 = 0 * u.m / u.s
-        self.v_drift3 = 1e5 * u.m / u.s
+        self.V_drift = 1000000 * u.m / u.s
+        self.V_drift2 = 0 * u.m / u.s
+        self.V_drift3 = 1e5 * u.m / u.s
         self.start = -5000
         self.stop = - self.start
         self.dv = 10000 * u.m / u.s
@@ -49,7 +47,7 @@ def test_max_noDrift(self):
         max_index = Maxwellian_1D(self.v_vect,
                                   T=self.T_e,
                                   particle=self.particle,
-                                  v_drift=0 * u.m / u.s
+                                  V_drift=0 * u.m / u.s
                                   ).argmax()
         assert np.isclose(self.v_vect[max_index].value, 0.0)
 
@@ -61,9 +59,9 @@ def test_max_drift(self):
         max_index = Maxwellian_1D(self.v_vect,
                                   T=self.T_e,
                                   particle=self.particle,
-                                  v_drift=self.v_drift
+                                  V_drift=self.V_drift
                                   ).argmax()
-        assert np.isclose(self.v_vect[max_index].value, self.v_drift.value)
+        assert np.isclose(self.v_vect[max_index].value, self.V_drift.value)
 
     def test_norm(self):
         """
@@ -182,7 +180,7 @@ def test_zero_drift_units(self):
         distFunc = Maxwellian_1D(v=self.v,
                                  T=self.T_e,
                                  particle=self.particle,
-                                 v_drift=self.v_drift2,
+                                 V_drift=self.V_drift2,
                                  units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -199,7 +197,7 @@ def test_value_drift_units(self):
         distFunc = Maxwellian_1D(v=self.v,
                                  T=self.T_e,
                                  particle=self.particle,
-                                 v_drift=self.v_drift3,
+                                 V_drift=self.V_drift3,
                                  units="units")
         errStr = (f"Distribution function should be {testVal} "
                   f"and not {distFunc}.")
@@ -221,10 +219,9 @@ def setup_class(self):
                                  particle=self.particle,
                                  method="most_probable")
         self.v = 1e5 * u.m / u.s
-        self.v_drift = 0 * u.m / u.s
-        self.v_drift2 = 1e5 * u.m / u.s
-        self.distFuncTrue = 1.72940389716217e-27
-        self.distFuncDrift = 2 * (self.vTh ** 2 * np.pi) ** (-1 / 2)
+        self.V_drift = 0 * u.m / u.s
+        self.V_drift2 = 1e5 * u.m / u.s
+        self.distFuncTrue = 1.8057567503860518e-25
 
     def test_norm(self):
         """
@@ -316,7 +313,7 @@ def test_zero_drift_units(self):
         distFunc = Maxwellian_speed_1D(v=self.v,
                                        T=self.T,
                                        particle=self.particle,
-                                       v_drift=self.v_drift,
+                                       V_drift=self.V_drift,
                                        units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -332,18 +329,18 @@ def test_value_drift_units(self):
         distFunc = Maxwellian_speed_1D(v=self.v,
                                        T=self.T,
                                        particle=self.particle,
-                                       v_drift=self.v_drift2,
+                                       V_drift=self.V_drift2,
                                        units="units")
         errStr = (f"Distribution function should be 0.0 "
                   f"and not {distFunc}.")
         assert np.isclose(distFunc.value,
-                          self.distFuncDrift.value,
+                          0.0,
                           rtol=1e-8,
                           atol=0.0), errStr
 
 
-# test class for Maxwellian_velocity_2D function
-class Test_Maxwellian_velocity_2D(object):
+# test class for Maxwellian_velocity_3D function
+class Test_Maxwellian_velocity_3D(object):
     @classmethod
     def setup_class(self):
         """initializing parameters for tests """
@@ -355,11 +352,14 @@ def setup_class(self):
                                  method="most_probable")
         self.vx = 1e5 * u.m / u.s
         self.vy = 1e5 * u.m / u.s
-        self.vx_drift = 0 * u.m / u.s
-        self.vy_drift = 0 * u.m / u.s
-        self.vx_drift2 = 1e5 * u.m / u.s
-        self.vy_drift2 = 1e5 * u.m / u.s
-        self.distFuncTrue = 7.477094598799251e-55
+        self.vz = 1e5 * u.m / u.s
+        self.Vx_drift = 0 * u.m / u.s
+        self.Vy_drift = 0 * u.m / u.s
+        self.Vz_drift = 0 * u.m / u.s
+        self.Vx_drift2 = 1e5 * u.m / u.s
+        self.Vy_drift2 = 1e5 * u.m / u.s
+        self.Vz_drift2 = 1e5 * u.m / u.s
+        self.distFuncTrue = 2.7654607627522045e-37
 
     def test_norm(self):
         """
@@ -370,21 +370,24 @@ def test_norm(self):
         # setting up integration from -10*vTh to 10*vTh, which is close to Inf
         infApprox = (10 * vTh)
         # integrating, this should be close to 1
-        integ = spint.dblquad(Maxwellian_velocity_2D,
+        integ = spint.tplquad(Maxwellian_velocity_3D,
                               -infApprox,
                               infApprox,
-                              lambda y: -infApprox,
-                              lambda y: infApprox,
+                              lambda z: -infApprox,
+                              lambda z: infApprox,
+                              lambda z, y: -infApprox,
+                              lambda z, y: infApprox,
                               args=(self.T,
                                     self.particle,
                                     0,
                                     0,
+                                    0,
                                     vTh,
                                     "unitless"),
                               epsabs=1e0,
                               epsrel=1e0,
                               )
-        # value returned from dblquad is (integral, error), we just need
+        # value returned from tplquad is (integral, error), we just need
         # the 1st
         integVal = integ[0]
         exceptStr = ("Integral of distribution function should be 1 "
@@ -398,8 +401,9 @@ def test_units_no_vTh(self):
         """
         Tests distribution function with units, but not passing vTh.
         """
-        distFunc = Maxwellian_velocity_2D(vx=self.vx,
+        distFunc = Maxwellian_velocity_3D(vx=self.vx,
                                           vy=self.vy,
+                                          vz=self.vz,
                                           T=self.T,
                                           particle=self.particle,
                                           units="units")
@@ -414,8 +418,9 @@ def test_units_vTh(self):
         """
         Tests distribution function with units and passing vTh.
         """
-        distFunc = Maxwellian_velocity_2D(vx=self.vx,
+        distFunc = Maxwellian_velocity_3D(vx=self.vx,
                                           vy=self.vy,
+                                          vz=self.vz,
                                           T=self.T,
                                           vTh=self.vTh,
                                           particle=self.particle,
@@ -434,8 +439,9 @@ def test_unitless_no_vTh(self):
         # converting T to SI then stripping units
         T = self.T.to(u.K, equivalencies=u.temperature_energy())
         T = T.si.value
-        distFunc = Maxwellian_velocity_2D(vx=self.vx.si.value,
+        distFunc = Maxwellian_velocity_3D(vx=self.vx.si.value,
                                           vy=self.vy.si.value,
+                                          vz=self.vz.si.value,
                                           T=T,
                                           particle=self.particle,
                                           units="unitless")
@@ -453,8 +459,9 @@ def test_unitless_vTh(self):
         # converting T to SI then stripping units
         T = self.T.to(u.K, equivalencies=u.temperature_energy())
         T = T.si.value
-        distFunc = Maxwellian_velocity_2D(vx=self.vx.si.value,
+        distFunc = Maxwellian_velocity_3D(vx=self.vx.si.value,
                                           vy=self.vy.si.value,
+                                          vz=self.vz.si.value,
                                           T=T,
                                           vTh=self.vTh.si.value,
                                           particle=self.particle,
@@ -471,12 +478,14 @@ def test_zero_drift_units(self):
         Testing inputting drift equal to 0 with units. These should just
         get passed and not have extra units applied to them.
         """
-        distFunc = Maxwellian_velocity_2D(vx=self.vx,
+        distFunc = Maxwellian_velocity_3D(vx=self.vx,
                                           vy=self.vy,
+                                          vz=self.vz,
                                           T=self.T,
                                           particle=self.particle,
-                                          vx_drift=self.vx_drift,
-                                          vy_drift=self.vy_drift,
+                                          Vx_drift=self.Vx_drift,
+                                          Vy_drift=self.Vy_drift,
+                                          Vz_drift=self.Vz_drift,
                                           units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -489,13 +498,15 @@ def test_value_drift_units(self):
         """
         Testing vdrifts with values
         """
-        testVal = ((self.vTh**2 * np.pi) ** (-1)).si.value
-        distFunc = Maxwellian_velocity_2D(vx=self.vx,
+        testVal = ((3 * self.vTh**2 * np.pi) ** (-3 / 2)).si.value
+        distFunc = Maxwellian_velocity_3D(vx=self.vx,
                                           vy=self.vy,
+                                          vz=self.vz,
                                           T=self.T,
                                           particle=self.particle,
-                                          vx_drift=self.vx_drift2,
-                                          vy_drift=self.vy_drift2,
+                                          Vx_drift=self.Vx_drift2,
+                                          Vy_drift=self.Vy_drift2,
+                                          Vz_drift=self.Vz_drift2,
                                           units="units")
         errStr = (f"Distribution function should be {testVal} "
                   f"and not {distFunc}.")
@@ -505,141 +516,8 @@ def test_value_drift_units(self):
                           atol=0.0), errStr
 
 
-# test class for Maxwellian_speed_2D function
-class Test_Maxwellian_speed_2D(object):
-    @classmethod
-    def setup_class(self):
-        """initializing parameters for tests """
-        self.T = 1.0 * u.eV
-        self.particle = 'H+'
-        # get thermal velocity and thermal velocity squared
-        self.vTh = thermal_speed(self.T,
-                                 particle=self.particle,
-                                 method="most_probable")
-        self.v = 1e5 * u.m / u.s
-        self.v_drift = 0 * u.m / u.s
-        self.v_drift2 = 1e5 * u.m / u.s
-        self.distFuncTrue = 2.2148166449365907e-26
-
-    def test_norm(self):
-        """
-        Tests whether distribution function is normalized, and integrates to 1.
-        """
-        # setting up integration from 0 to 10*vTh
-        xData1D = np.arange(0, 10.001, 0.001) * self.vTh
-        yData1D = Maxwellian_speed_2D(v=xData1D,
-                                      T=self.T,
-                                      particle=self.particle)
-        # integrating, this should be close to 1
-        integ = spint.trapz(y=yData1D, x=xData1D)
-        exceptStr = "Integral of distribution function should be 1."
-        assert np.isclose(integ.value, 1), exceptStr
-
-    def test_units_no_vTh(self):
-        """
-        Tests distribution function with units, but not passing vTh.
-        """
-        distFunc = Maxwellian_speed_2D(v=self.v,
-                                       T=self.T,
-                                       particle=self.particle,
-                                       units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_units_vTh(self):
-        """
-        Tests distribution function with units and passing vTh.
-        """
-        distFunc = Maxwellian_speed_2D(v=self.v,
-                                       T=self.T,
-                                       vTh=self.vTh,
-                                       particle=self.particle,
-                                       units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_unitless_no_vTh(self):
-        """
-        Tests distribution function without units, and not passing vTh.
-        """
-        # converting T to SI then stripping units
-        T = self.T.to(u.K, equivalencies=u.temperature_energy())
-        T = T.si.value
-        distFunc = Maxwellian_speed_2D(v=self.v.si.value,
-                                       T=T,
-                                       particle=self.particle,
-                                       units="unitless")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_unitless_vTh(self):
-        """
-        Tests distribution function without units, and with passing vTh.
-        """
-        # converting T to SI then stripping units
-        T = self.T.to(u.K, equivalencies=u.temperature_energy())
-        T = T.si.value
-        distFunc = Maxwellian_speed_2D(v=self.v.si.value,
-                                       T=T,
-                                       vTh=self.vTh.si.value,
-                                       particle=self.particle,
-                                       units="unitless")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_zero_drift_units(self):
-        """
-        Testing inputting drift equal to 0 with units. These should just
-        get passed and not have extra units applied to them.
-        """
-        distFunc = Maxwellian_speed_2D(v=self.v,
-                                       T=self.T,
-                                       particle=self.particle,
-                                       v_drift=self.v_drift,
-                                       units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_value_drift_units(self):
-        """
-        Testing vdrifts with values
-        """
-        with pytest.raises(NotImplementedError):
-            distFunc = Maxwellian_speed_2D(v=self.v,
-                                           T=self.T,
-                                           particle=self.particle,
-                                           v_drift=self.v_drift2,
-                                           units="units")
-#        errStr = (f"Distribution function should be 0.0 "
-#                  f"and not {distFunc}.")
-#        assert np.isclose(distFunc.value,
-#                          0.0,
-#                          rtol=1e-8,
-#                          atol=0.0), errStr
-
-
-# test class for Maxwellian_velocity_3D function
-class Test_Maxwellian_velocity_3D(object):
+# test class for Maxwellian_speed_3D function
+class Test_Maxwellian_speed_3D(object):
     @classmethod
     def setup_class(self):
         """initializing parameters for tests """
@@ -652,13 +530,13 @@ def setup_class(self):
         self.vx = 1e5 * u.m / u.s
         self.vy = 1e5 * u.m / u.s
         self.vz = 1e5 * u.m / u.s
-        self.vx_drift = 0 * u.m / u.s
-        self.vy_drift = 0 * u.m / u.s
-        self.vz_drift = 0 * u.m / u.s
-        self.vx_drift2 = 1e5 * u.m / u.s
-        self.vy_drift2 = 1e5 * u.m / u.s
-        self.vz_drift2 = 1e5 * u.m / u.s
-        self.distFuncTrue = 6.465458269306909e-82
+        self.Vx_drift = 0 * u.m / u.s
+        self.Vy_drift = 0 * u.m / u.s
+        self.Vz_drift = 0 * u.m / u.s
+        self.Vx_drift2 = 1e5 * u.m / u.s
+        self.Vy_drift2 = 1e5 * u.m / u.s
+        self.Vz_drift2 = 1e5 * u.m / u.s
+        self.distFuncTrue = 5.888134761477178e-75
 
     def test_norm(self):
         """
@@ -666,15 +544,15 @@ def test_norm(self):
         """
         # converting vTh to unitless
         vTh = self.vTh.si.value
-        # setting up integration from -10*vTh to 10*vTh, which is close to Inf
+        # setting up integration from 0 to 10*vTh, which is close to Inf
         infApprox = (10 * vTh)
-        # integrating, this should be close to 1
-        integ = spint.tplquad(Maxwellian_velocity_3D,
-                              -infApprox,
+        # integral should be close to 1
+        integ = spint.tplquad(Maxwellian_speed_3D,
+                              0,
                               infApprox,
-                              lambda z: -infApprox,
+                              lambda z: 0,
                               lambda z: infApprox,
-                              lambda z, y: -infApprox,
+                              lambda z, y: 0,
                               lambda z, y: infApprox,
                               args=(self.T,
                                     self.particle,
@@ -690,7 +568,7 @@ def test_norm(self):
         # the 1st
         integVal = integ[0]
         exceptStr = ("Integral of distribution function should be 1 "
-                     f"and not {integVal}.")
+                     f"and not {integVal}")
         assert np.isclose(integVal,
                           1,
                           rtol=1e-3,
@@ -700,156 +578,9 @@ def test_units_no_vTh(self):
         """
         Tests distribution function with units, but not passing vTh.
         """
-        distFunc = Maxwellian_velocity_3D(vx=self.vx,
-                                          vy=self.vy,
-                                          vz=self.vz,
-                                          T=self.T,
-                                          particle=self.particle,
-                                          units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_units_vTh(self):
-        """
-        Tests distribution function with units and passing vTh.
-        """
-        distFunc = Maxwellian_velocity_3D(vx=self.vx,
-                                          vy=self.vy,
-                                          vz=self.vz,
-                                          T=self.T,
-                                          vTh=self.vTh,
-                                          particle=self.particle,
-                                          units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_unitless_no_vTh(self):
-        """
-        Tests distribution function without units, and not passing vTh.
-        """
-        # converting T to SI then stripping units
-        T = self.T.to(u.K, equivalencies=u.temperature_energy())
-        T = T.si.value
-        distFunc = Maxwellian_velocity_3D(vx=self.vx.si.value,
-                                          vy=self.vy.si.value,
-                                          vz=self.vz.si.value,
-                                          T=T,
-                                          particle=self.particle,
-                                          units="unitless")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_unitless_vTh(self):
-        """
-        Tests distribution function without units, and with passing vTh.
-        """
-        # converting T to SI then stripping units
-        T = self.T.to(u.K, equivalencies=u.temperature_energy())
-        T = T.si.value
-        distFunc = Maxwellian_velocity_3D(vx=self.vx.si.value,
-                                          vy=self.vy.si.value,
-                                          vz=self.vz.si.value,
-                                          T=T,
-                                          vTh=self.vTh.si.value,
-                                          particle=self.particle,
-                                          units="unitless")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_zero_drift_units(self):
-        """
-        Testing inputting drift equal to 0 with units. These should just
-        get passed and not have extra units applied to them.
-        """
-        distFunc = Maxwellian_velocity_3D(vx=self.vx,
-                                          vy=self.vy,
-                                          vz=self.vz,
-                                          T=self.T,
-                                          particle=self.particle,
-                                          vx_drift=self.vx_drift,
-                                          vy_drift=self.vy_drift,
-                                          vz_drift=self.vz_drift,
-                                          units="units")
-        errStr = (f"Distribution function should be {self.distFuncTrue} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          self.distFuncTrue,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-    def test_value_drift_units(self):
-        """
-        Testing vdrifts with values
-        """
-        testVal = ((self.vTh**2 * np.pi) ** (-3 / 2)).si.value
-        distFunc = Maxwellian_velocity_3D(vx=self.vx,
-                                          vy=self.vy,
-                                          vz=self.vz,
-                                          T=self.T,
-                                          particle=self.particle,
-                                          vx_drift=self.vx_drift2,
-                                          vy_drift=self.vy_drift2,
-                                          vz_drift=self.vz_drift2,
-                                          units="units")
-        errStr = (f"Distribution function should be {testVal} "
-                  f"and not {distFunc}.")
-        assert np.isclose(distFunc.value,
-                          testVal,
-                          rtol=1e-8,
-                          atol=0.0), errStr
-
-
-# test class for Maxwellian_speed_3D function
-class Test_Maxwellian_speed_3D(object):
-    @classmethod
-    def setup_class(self):
-        """initializing parameters for tests """
-        self.T = 1.0 * u.eV
-        self.particle = 'H+'
-        # get thermal velocity and thermal velocity squared
-        self.vTh = thermal_speed(self.T,
-                                 particle=self.particle,
-                                 method="most_probable")
-        self.v = 1e5 * u.m / u.s
-        self.v_drift = 0 * u.m / u.s
-        self.v_drift2 = 1e5 * u.m / u.s
-        self.distFuncTrue = 1.8057567503860518e-25
-
-    def test_norm(self):
-        """
-        Tests whether distribution function is normalized, and integrates to 1.
-        """
-        # setting up integration from 0 to 10*vTh
-        xData1D = np.arange(0, 10.01, 0.01) * self.vTh
-        yData1D = Maxwellian_speed_3D(v=xData1D,
-                                      T=self.T,
-                                      particle=self.particle)
-        # integrating, this should be close to 1
-        integ = spint.trapz(y=yData1D, x=xData1D)
-        exceptStr = "Integral of distribution function should be 1."
-        assert np.isclose(integ.value, 1), exceptStr
-
-    def test_units_no_vTh(self):
-        """
-        Tests distribution function with units, but not passing vTh.
-        """
-        distFunc = Maxwellian_speed_3D(v=self.v,
+        distFunc = Maxwellian_speed_3D(vx=self.vx,
+                                       vy=self.vy,
+                                       vz=self.vz,
                                        T=self.T,
                                        particle=self.particle,
                                        units="units")
@@ -864,7 +595,9 @@ def test_units_vTh(self):
         """
         Tests distribution function with units and passing vTh.
         """
-        distFunc = Maxwellian_speed_3D(v=self.v,
+        distFunc = Maxwellian_speed_3D(vx=self.vx,
+                                       vy=self.vy,
+                                       vz=self.vz,
                                        T=self.T,
                                        vTh=self.vTh,
                                        particle=self.particle,
@@ -883,7 +616,9 @@ def test_unitless_no_vTh(self):
         # converting T to SI then stripping units
         T = self.T.to(u.K, equivalencies=u.temperature_energy())
         T = T.si.value
-        distFunc = Maxwellian_speed_3D(v=self.v.si.value,
+        distFunc = Maxwellian_speed_3D(vx=self.vx.si.value,
+                                       vy=self.vy.si.value,
+                                       vz=self.vz.si.value,
                                        T=T,
                                        particle=self.particle,
                                        units="unitless")
@@ -901,7 +636,9 @@ def test_unitless_vTh(self):
         # converting T to SI then stripping units
         T = self.T.to(u.K, equivalencies=u.temperature_energy())
         T = T.si.value
-        distFunc = Maxwellian_speed_3D(v=self.v.si.value,
+        distFunc = Maxwellian_speed_3D(vx=self.vx.si.value,
+                                       vy=self.vy.si.value,
+                                       vz=self.vz.si.value,
                                        T=T,
                                        vTh=self.vTh.si.value,
                                        particle=self.particle,
@@ -918,10 +655,14 @@ def test_zero_drift_units(self):
         Testing inputting drift equal to 0 with units. These should just
         get passed and not have extra units applied to them.
         """
-        distFunc = Maxwellian_speed_3D(v=self.v,
+        distFunc = Maxwellian_speed_3D(vx=self.vx,
+                                       vy=self.vy,
+                                       vz=self.vz,
                                        T=self.T,
                                        particle=self.particle,
-                                       v_drift=self.v_drift,
+                                       Vx_drift=self.Vx_drift,
+                                       Vy_drift=self.Vy_drift,
+                                       Vz_drift=self.Vz_drift,
                                        units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -934,18 +675,21 @@ def test_value_drift_units(self):
         """
         Testing vdrifts with values
         """
-        with pytest.raises(NotImplementedError):
-            distFunc = Maxwellian_speed_3D(v=self.v,
-                                           T=self.T,
-                                           particle=self.particle,
-                                           v_drift=self.v_drift2,
-                                           units="units")
-#        errStr = (f"Distribution function should be 0.0 "
-#                  f"and not {distFunc}.")
-#        assert np.isclose(distFunc.value,
-#                          0.0,
-#                          rtol=1e-8,
-#                          atol=0.0), errStr
+        distFunc = Maxwellian_speed_3D(vx=self.vx,
+                                       vy=self.vy,
+                                       vz=self.vz,
+                                       T=self.T,
+                                       particle=self.particle,
+                                       Vx_drift=self.Vx_drift2,
+                                       Vy_drift=self.Vy_drift2,
+                                       Vz_drift=self.Vz_drift2,
+                                       units="units")
+        errStr = (f"Distribution function should be 0.0 "
+                  f"and not {distFunc}.")
+        assert np.isclose(distFunc.value,
+                          0.0,
+                          rtol=1e-8,
+                          atol=0.0), errStr
 
 
 # kappa
@@ -959,9 +703,9 @@ def setup_class(self):
         self.kappa = 4
         self.kappaInvalid = 3 / 2
         self.v = 1e5 * u.m / u.s
-        self.v_drift = 1000000 * u.m / u.s
-        self.v_drift2 = 0 * u.m / u.s
-        self.v_drift3 = 1e5 * u.m / u.s
+        self.V_drift = 1000000 * u.m / u.s
+        self.V_drift2 = 0 * u.m / u.s
+        self.V_drift3 = 1e5 * u.m / u.s
         self.start = -5000
         self.stop = -self.start
         self.dv = 10000 * u.m / u.s
@@ -995,7 +739,7 @@ def test_max_noDrift(self):
                                       T=self.T_e,
                                       kappa=self.kappa,
                                       particle=self.particle,
-                                      v_drift=0 * u.m / u.s).argmax()
+                                      V_drift=0 * u.m / u.s).argmax()
         assert np.isclose(self.v_vect[max_index].value, 0.0)
 
     def test_max_drift(self):
@@ -1007,8 +751,8 @@ def test_max_drift(self):
                                       T=self.T_e,
                                       kappa=self.kappa,
                                       particle=self.particle,
-                                      v_drift=self.v_drift).argmax()
-        assert np.isclose(self.v_vect[max_index].value, self.v_drift.value)
+                                      V_drift=self.V_drift).argmax()
+        assert np.isclose(self.v_vect[max_index].value, self.V_drift.value)
 
     def test_maxwellian_limit(self):
         """
@@ -1140,7 +884,7 @@ def test_zero_drift_units(self):
                                      T=self.T_e,
                                      kappa=self.kappa,
                                      particle=self.particle,
-                                     v_drift=self.v_drift2,
+                                     V_drift=self.V_drift2,
                                      units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -1158,7 +902,7 @@ def test_value_drift_units(self):
                                      T=self.T_e,
                                      kappa=self.kappa,
                                      particle=self.particle,
-                                     v_drift=self.v_drift3,
+                                     V_drift=self.V_drift3,
                                      units="units")
         errStr = (f"Distribution function should be {testVal} "
                   f"and not {distFunc}.")
@@ -1184,12 +928,12 @@ def setup_class(self):
         self.vx = 1e5 * u.m / u.s
         self.vy = 1e5 * u.m / u.s
         self.vz = 1e5 * u.m / u.s
-        self.vx_drift = 0 * u.m / u.s
-        self.vy_drift = 0 * u.m / u.s
-        self.vz_drift = 0 * u.m / u.s
-        self.vx_drift2 = 1e5 * u.m / u.s
-        self.vy_drift2 = 1e5 * u.m / u.s
-        self.vz_drift2 = 1e5 * u.m / u.s
+        self.Vx_drift = 0 * u.m / u.s
+        self.Vy_drift = 0 * u.m / u.s
+        self.Vz_drift = 0 * u.m / u.s
+        self.Vx_drift2 = 1e5 * u.m / u.s
+        self.Vy_drift2 = 1e5 * u.m / u.s
+        self.Vz_drift2 = 1e5 * u.m / u.s
         self.distFuncTrue = 1.1847914288918793e-22
 
     def test_invalid_kappa(self):
@@ -1218,9 +962,9 @@ def test_invalid_kappa(self):
 #                                          T=self.T,
 #                                          kappa=kappaLarge,
 #                                          particle=self.particle,
-#                                          vx_drift=self.vx_drift2,
-#                                          vy_drift=self.vy_drift2,
-#                                          vz_drift=self.vz_drift2,
+#                                          Vx_drift=self.Vx_drift2,
+#                                          Vy_drift=self.Vy_drift2,
+#                                          Vz_drift=self.Vz_drift2,
 #                                          units="units")
 #        Teff =  self.T
 #        maxwellDistFunc = Maxwellian_velocity_3D(vx=self.vx,
@@ -1228,9 +972,9 @@ def test_invalid_kappa(self):
 #                                                 vz=self.vz,
 #                                                 T=Teff,
 #                                                 particle=self.particle,
-#                                                 vx_drift=self.vx_drift2,
-#                                                 vy_drift=self.vy_drift2,
-#                                                 vz_drift=self.vz_drift2,
+#                                                 Vx_drift=self.Vx_drift2,
+#                                                 Vy_drift=self.Vy_drift2,
+#                                                 Vz_drift=self.Vz_drift2,
 #                                                 units="units")
 #        errStr = (f"Distribution function should be {maxwellDistFunc} "
 #                  f"and not {kappaDistFunc}.")
@@ -1369,9 +1113,9 @@ def test_zero_drift_units(self):
                                      T=self.T,
                                      kappa=self.kappa,
                                      particle=self.particle,
-                                     vx_drift=self.vx_drift,
-                                     vy_drift=self.vy_drift,
-                                     vz_drift=self.vz_drift,
+                                     Vx_drift=self.Vx_drift,
+                                     Vy_drift=self.Vy_drift,
+                                     Vz_drift=self.Vz_drift,
                                      units="units")
         errStr = (f"Distribution function should be {self.distFuncTrue} "
                   f"and not {distFunc}.")
@@ -1391,9 +1135,9 @@ def test_value_drift_units(self):
                                      T=self.T,
                                      kappa=self.kappa,
                                      particle=self.particle,
-                                     vx_drift=self.vx_drift2,
-                                     vy_drift=self.vy_drift2,
-                                     vz_drift=self.vz_drift2,
+                                     Vx_drift=self.Vx_drift2,
+                                     Vy_drift=self.Vy_drift2,
+                                     Vz_drift=self.Vz_drift2,
                                      units="units")
         errStr = (f"Distribution function should be {testVal} "
                   f"and not {distFunc}.")
diff --git a/plasmapy/utils/tests/test_datatype_factory_base.py b/plasmapy/utils/tests/test_datatype_factory_base.py
new file mode 100644
index 00000000..cabce995
--- /dev/null
+++ b/plasmapy/utils/tests/test_datatype_factory_base.py
@@ -0,0 +1,183 @@
+import pytest
+
+from plasmapy.utils.datatype_factory_base import BasicRegistrationFactory
+from plasmapy.utils.datatype_factory_base import NoMatchError
+from plasmapy.utils.datatype_factory_base import MultipleMatchError
+from plasmapy.utils.datatype_factory_base import ValidationFunctionError
+
+# SunPy is released under a BSD-style open source licence:
+
+# Copyright (c) 2013-2018 The SunPy developers
+# All rights reserved.
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+
+# * Redistributions of source code must retain the above copyright
+#   notice, this list of conditions and the following disclaimer.
+
+# * Redistributions in binary form must reproduce the above copyright
+#   notice, this list of conditions and the following disclaimer in the
+#   documentation and/or other materials provided with the distribution.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+# A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT
+# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+# Project     : https://github.com/sunpy/sunpy
+# File        : sunpy/util/tests/test_datatype_factory_base.py
+# Commit hash : f6330eea602ea796b5b004dee283b8877b24da23
+
+
+class BaseWidget(object):
+    def __init__(self, *args, **kwargs):
+        pass
+
+
+class DefaultWidget(BaseWidget):
+    pass
+
+
+class StandardWidget(BaseWidget):
+    @classmethod
+    def _factory_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'standard'
+
+
+class DuplicateStandardWidget(BaseWidget):
+    @classmethod
+    def _factory_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'standard'
+
+
+class FancyWidget(BaseWidget):
+    @classmethod
+    def _factory_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'fancy' and 'feature' in kwargs
+
+
+class ExternallyValidatedWidget(BaseWidget):
+    pass
+
+
+def external_validation_function(*args, **kwargs):
+    return kwargs.get('style') == 'external'
+
+
+class UnvalidatedWidget(BaseWidget):
+    pass
+
+
+class MissingClassMethodWidget(BaseWidget):
+    def _factory_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'missing'
+
+
+class DifferentValidationWidget(BaseWidget):
+    @classmethod
+    def different_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'different'
+
+
+class MissingClassMethodDifferentValidationWidget(BaseWidget):
+    def different_validation_function(cls, *args, **kwargs):
+        return kwargs.get('style') == 'missing-different'
+
+
+class TestBasicRegistrationFactory(object):
+
+    def test_default_factory(self):
+
+        DefaultFactory = BasicRegistrationFactory()
+
+        DefaultFactory.register(DefaultWidget, is_default=True)
+        assert DefaultFactory.default_widget_type == DefaultWidget
+
+        DefaultFactory.register(StandardWidget)
+        DefaultFactory.register(FancyWidget)
+        DefaultFactory.register(ExternallyValidatedWidget,
+                                validation_function=external_validation_function)
+
+        assert type(DefaultFactory()) is DefaultWidget
+        assert type(DefaultFactory(style='standard')) is StandardWidget
+        assert type(DefaultFactory(style='fancy')) is DefaultWidget
+        assert type(DefaultFactory(style='fancy', feature="present")) is FancyWidget
+        assert type(DefaultFactory(style='external')) is ExternallyValidatedWidget
+
+        with pytest.raises(ValidationFunctionError):
+            DefaultFactory.register(UnvalidatedWidget)
+
+        with pytest.raises(ValidationFunctionError):
+            DefaultFactory.register(MissingClassMethodWidget)
+
+        DefaultFactory.unregister(StandardWidget)
+        assert type(DefaultFactory(style='standard')) is not StandardWidget
+
+    def test_validation_fun_not_callable(self):
+        TestFactory = BasicRegistrationFactory()
+
+        with pytest.raises(AttributeError):
+            TestFactory.register(StandardWidget, validation_function='not_callable')
+
+    def test_no_default_factory(self):
+
+        NoDefaultFactory = BasicRegistrationFactory()
+
+        NoDefaultFactory.register(StandardWidget)
+        NoDefaultFactory.register(FancyWidget)
+
+        with pytest.raises(NoMatchError):
+            NoDefaultFactory()
+
+        # Raises because all requirements are not met for FancyWidget and no
+        # default is present.
+        with pytest.raises(NoMatchError):
+            NoDefaultFactory(style='fancy')
+
+        assert type(NoDefaultFactory(style='standard')) is StandardWidget
+        assert type(NoDefaultFactory(style='fancy', feature='present')) is FancyWidget
+
+    def test_with_external_registry(self):
+        external_registry = {}
+
+        FactoryWithExternalRegistry = \
+            BasicRegistrationFactory(registry=external_registry)
+
+        assert len(external_registry) == 0
+
+        FactoryWithExternalRegistry.register(StandardWidget)
+        assert type(FactoryWithExternalRegistry(style='standard')) is StandardWidget
+
+        # Ensure the 'external_registry' is being populated see #1988
+        assert len(external_registry) == 1
+
+    def test_multiple_match_factory(self):
+
+        MultipleMatchFactory = BasicRegistrationFactory()
+
+        MultipleMatchFactory.register(StandardWidget)
+        MultipleMatchFactory.register(DuplicateStandardWidget)
+
+        with pytest.raises(MultipleMatchError):
+            MultipleMatchFactory(style='standard')
+
+    def test_extra_validation_factory(self):
+        ExtraValidationFactory = \
+            BasicRegistrationFactory(
+                additional_validation_functions=['different_validation_function'])
+
+        ExtraValidationFactory.register(DifferentValidationWidget)
+
+        assert type(ExtraValidationFactory(style='different')) is DifferentValidationWidget
+
+        with pytest.raises(ValidationFunctionError):
+            ExtraValidationFactory.register(MissingClassMethodDifferentValidationWidget)