Repository: PlasmaPy/PlasmaPy
Pull request #908 (2020-10-02 19:16:54)

Create plasmapy.analysis sub-pkg & add fit functions as initial functionality
Fit functions are a set of callable classes that wraps together the following functionality to perform curve fitting...

1. A callable version of the fitted function.
1. A curve fitting method to define the fit function parameters and associated errors.
1. Error propagation for the fitted curve
1. Root solver for the fitted function

This functionality was originally developed for the swept Langmuir analysis in PR #889, but soon went beyond the PR's scope so I broke it out into this PR.  This PR has to be merged before PR #889.

----

This PR does the following...

- [x] Creates `plasmapy.analysis`.
- [x] Initialize `plasmapy.analysis` with `fit_functions.py`.
- [x] Updates the documentation to incorporate the "Analysis" documentation with the "Diagnostics" documentation.
- [x] Adds associated tests.
- [x] Add a jupyter notebook example.


codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=h1) Report
> Merging [#908](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/6d5651cb9cfb7bd010c8e81c45a1bffeb9c4e95d?el=desc) will **decrease** coverage by `2.09%`.
> The diff coverage is `47.32%`.

[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908/graphs/tree.svg?width=650&height=150&src=pr&token=LzZ5mnXVad)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #908      +/-   ##
==========================================
- Coverage   96.21%   94.11%   -2.10%     
==========================================
  Files          60       62       +2     
  Lines        5417     5660     +243     
==========================================
+ Hits         5212     5327     +115     
- Misses        205      333     +128     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/analysis/fit\_functions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908/diff?src=pr&el=tree#diff-cGxhc21hcHkvYW5hbHlzaXMvZml0X2Z1bmN0aW9ucy5weQ==) | `46.88% <46.88%> (ø)` | |
| [plasmapy/analysis/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908/diff?src=pr&el=tree#diff-cGxhc21hcHkvYW5hbHlzaXMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=footer). Last update [6d5651c...624aa15](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/908?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

rocco8773: @namurphy I created a little `@modify_docstring` decorator to handle the docstings of `func_err` and made a `_check_func_err_params` method to help streamline the overriding process.  Let me know what you think.  If you like the decorator approach, then I'll do a PR for that separately before we merge this PR.
rocco8773: > Looks pretty much ready to merge, besides pre-commit formatting changes. Thank you for working on this!

The problem is, what makes `isort` happy makes `black` grumpy, visa-versa.  The only way to  make both happy is to make the import statements look like...

```python
from abc import ABC, abstractmethod
from contextlib import ExitStack as does_not_raise

# ^ ExitStack can be replaced with nullcontext when we require >= python 3.7
```

But this disassociates the comment with the import statement it's suppose to be tied with. Thoughts?

The current version looks like...

```python
from abc import ABC, abstractmethod
from contextlib import ExitStack as does_not_raise
# ^ ExitStack can be replaced with nullcontext when we require >= python 3.7
```

and I'm happy leaving it like that until we drop python 3.6.
namurphy: > and I'm happy leaving it like that until we drop python 3.6.

Maybe we could drop Python 3.6 support? :smile_cat: 


Patch:
diff --git a/.coveragerc b/.coveragerc
deleted file mode 100644
index 47a30a6a..00000000
--- a/.coveragerc
+++ /dev/null
@@ -1,11 +0,0 @@
-[report]
-omit =
-    ci-helpers/*
-    */tests/*
-    plasmapy/setup_package.py
-    plasmapy/version.py
-    plasmapy/_base_init.py
-exclude_lines =
-    coverage: ignore
-    ImportError
-    ModuleNotFoundError
diff --git a/.gitignore b/.gitignore
index 878d978f..4ce6c4d8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,6 +77,7 @@ target/
 
 # Jupyter Notebook
 .ipynb_checkpoints
+notebooks_for_proto
 
 # pyenv
 .python-version
diff --git a/changelog/908.feature.rst b/changelog/908.feature.rst
new file mode 100644
index 00000000..0acd9222
--- /dev/null
+++ b/changelog/908.feature.rst
@@ -0,0 +1,6 @@
+Create the `plasmapy.analysis` package as per
+`PLEP-7 <https://github.com/PlasmaPy/PlasmaPy-PLEPs/blob/master/PLEP-0007.rst>`_ and
+initialize the package with the `~plasmapy.analysis.fit_functions` module.  Fit
+functions are designed to wrap together an analytical function, a curve fitter,
+uncertainty propagation, and a root solver to make curve fitting a little less
+painful.
diff --git a/docs/ad/ad_attention.inc b/docs/ad/ad_attention.inc
new file mode 100644
index 00000000..48f7036f
--- /dev/null
+++ b/docs/ad/ad_attention.inc
@@ -0,0 +1,6 @@
+.. important::
+   The Analysis and Diagnostic framework is in active development at the moment.
+   For the foreseeable future, the API will be in continuous flux as functionality
+   is added and modified.  To follow the package development please visit our
+   GitHub Project ( https://github.com/PlasmaPy/PlasmaPy/projects/19 ) and
+   comment on any of the relevant issues and/or pull requests.
diff --git a/docs/ad/analysis/index.rst b/docs/ad/analysis/index.rst
new file mode 100644
index 00000000..3e46e554
--- /dev/null
+++ b/docs/ad/analysis/index.rst
@@ -0,0 +1,41 @@
+.. include:: /ad/ad_attention.inc
+
+.. py:currentmodule:: plasmapy.analysis
+
+.. _analysis:
+
+===================
+`plasmapy.analysis`
+===================
+
+.. toctree::
+   :maxdepth: 1
+
+   Fit Functions </api_static/plasmapy.analysis.fit_functions>
+
+Notebooks
+---------
+
+.. nbgallery::
+   :glob:
+   :reversed:
+
+   /notebooks/analysis/*
+
+API
+---
+
+.. automodule:: plasmapy.analysis
+
+Sub-Packages & Modules
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+
+   fit_functions
+
+.. automodapi:: plasmapy.analysis
+   :no-main-docstr:
+   :no-heading:
+   :include-all-objects:
+   :headings: "^-"
diff --git a/docs/diagnostics/index.rst b/docs/ad/diagnostics/index.rst
similarity index 51%
rename from docs/diagnostics/index.rst
rename to docs/ad/diagnostics/index.rst
index 719a2dd1..40a7420f 100644
--- a/docs/diagnostics/index.rst
+++ b/docs/ad/diagnostics/index.rst
@@ -1,8 +1,10 @@
+.. include:: /ad/ad_attention.inc
+
 .. _diagnostics:
 
-*******************************************
-Plasma diagnostics (`plasmapy.diagnostics`)
-*******************************************
+======================
+`plasmapy.diagnostics`
+======================
 
 .. automodapi:: plasmapy.diagnostics
    :no-heading:
@@ -11,5 +13,5 @@ Plasma diagnostics (`plasmapy.diagnostics`)
    :maxdepth: 1
    :caption: Sub-Packages and Modules
 
-   langmuir
+   Swept Langmuir Analysis <langmuir>
    thomson
diff --git a/docs/diagnostics/langmuir.rst b/docs/ad/diagnostics/langmuir.rst
similarity index 88%
rename from docs/diagnostics/langmuir.rst
rename to docs/ad/diagnostics/langmuir.rst
index 3d71ec16..7f874569 100644
--- a/docs/diagnostics/langmuir.rst
+++ b/docs/ad/diagnostics/langmuir.rst
@@ -1,3 +1,5 @@
+.. include:: /ad/ad_attention.inc
+
 .. _langmuir:
 
 *****************
diff --git a/docs/diagnostics/thomson.rst b/docs/ad/diagnostics/thomson.rst
similarity index 100%
rename from docs/diagnostics/thomson.rst
rename to docs/ad/diagnostics/thomson.rst
diff --git a/docs/ad/index.rst b/docs/ad/index.rst
new file mode 100644
index 00000000..12521251
--- /dev/null
+++ b/docs/ad/index.rst
@@ -0,0 +1,43 @@
+.. include:: ad_attention.inc
+
+.. _ad:
+
+==============================
+Analysis & Diagnostic Toolkits
+==============================
+
+Analyses and diagnostics go hand-in-hand, but have subtle differences.  Thus,
+PlasmaPy gives each their own sub-packages, `plasmapy.analysis` and
+`plasmapy.diagnostics` respectively.
+
+
+Think of the `plasmapy.analysis` as your toolbox.  It has all the tools
+(functionality) you need to analyze your data.  Functionality is built around
+`numpy arrays <https://numpy.org/doc/stable/reference/arrays.html>`_ and each
+function has a well-defined, focused task.
+
+Diagnostics have a much broader scope and leverage the tools
+defined in `plasmapy.analysis` to give a more integrated user experience when
+analyzing data.  Diagnostics try to enhance the analysis workflow by focusing
+on some of following key areas...
+
+#. A more human-friendly way of managing data by building an interface around
+   `xarray` arrays and datasets via custom diagnostic accessors.
+
+   - `xarray` provides labelled multi-dimensional arrays and datasets.
+   - Diagnostics self-manage the computed analysis data within a `xarray`
+     dataset while maintaining the computed data's relation to the original
+     data.
+
+#. Quick viewing of analyzed data with default plotting routines.
+#. Fully defining the physical parameters of a diagnostic with purposely
+   designed `Probe` classes that are integrated into the analysis workflow.
+#. Adding graphical user interfaces (GUIs) to the analysis workflow via notebook widgets.
+
+----
+
+.. toctree::
+   :maxdepth: 2
+
+   Analysis <analysis/index>
+   Diagnostics <diagnostics/index>
diff --git a/docs/api_static/plasmapy.analysis.fit_functions.rst b/docs/api_static/plasmapy.analysis.fit_functions.rst
new file mode 100644
index 00000000..025eca8b
--- /dev/null
+++ b/docs/api_static/plasmapy.analysis.fit_functions.rst
@@ -0,0 +1,18 @@
+:orphan:
+
+`plasmapy.analysis.fit_functions`
+=================================
+
+.. currentmodule:: plasmapy.analysis.fit_functions
+
+.. automodapi:: plasmapy.analysis.fit_functions
+   :include-all-objects:
+   :inherited-members:
+   :no-heading:
+
+Example Notebook
+----------------
+
+.. nbgallery::
+
+   /notebooks/analysis/fit_functions
diff --git a/docs/index.rst b/docs/index.rst
index 1a6319bc..a94d827d 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -40,7 +40,7 @@ package for plasma physics currently under development.
    :maxdepth: 1
    :caption: Package features
 
-   Diagnostics <diagnostics/index>
+   ad/index
    Dispersion <dispersion/index>
    Formulary <formulary/index>
    Particles <particles/index>
diff --git a/docs/notebooks/analysis/fit_functions.ipynb b/docs/notebooks/analysis/fit_functions.ipynb
new file mode 100644
index 00000000..70fc6b52
--- /dev/null
+++ b/docs/notebooks/analysis/fit_functions.ipynb
@@ -0,0 +1,370 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fit Functions\n",
+    "\n",
+    "[Fit functions]: ../../api_static/plasmapy.analysis.fit_functions.rst\n",
+    "[linregress()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html\n",
+    "[curve_fit()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html\n",
+    "[fsolve()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html\n",
+    "\n",
+    "**[Fit functions]** are a set of callable classes designed to aid in fitting analytical functions to data.  A fit function class combines the following functionality:\n",
+    "\n",
+    "1. An analytical function that is callable with given parameters or fitted parameters.\n",
+    "1. Curve fitting functionality (usually SciPy's [curve_fit()] or [linregress()]), which stores the fit statistics and parameters into the class.  This makes the function easily callable with the fitted parameters.\n",
+    "1. Error propagation calculations.\n",
+    "1. A root solver that returns either the known analytical solutions or uses SciPy's [fsolve()] to calculate the roots."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from pathlib import Path\n",
+    "\n",
+    "from plasmapy.analysis import fit_functions as ffuncs\n",
+    "\n",
+    "plt.rcParams['figure.figsize'] = [10.5, 0.56 * 10.5]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Contents:\n",
+    "\n",
+    "1. [Fit function basics](#Fit-function-basics)\n",
+    "1. [Fitting to data](#Fitting-to-data)\n",
+    "    1. [Getting fit results](#Getting-fit-results)\n",
+    "    1. [Fit function is callable](#Fit-function-is-callable)\n",
+    "    1. [Plotting results](#Plotting-results)\n",
+    "    1. [Root solving](#Root-solving)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fit function basics\n",
+    "\n",
+    "[fit functions]: ../../api_static/plasmapy.analysis.fit_functions.rst\n",
+    "[ExponentialPlusLinear]: ../../api/plasmapy.analysis.fit_functions.ExponentialPlusLinear.rst\n",
+    "\n",
+    "There is an ever expanding collection of [fit functions], but this notebook will use [ExponentialPlusLinear] as an example.\n",
+    "\n",
+    "A fit function class has no required arguments at time of instantiation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# basic instantiation\n",
+    "explin = ffuncs.ExponentialPlusLinear()\n",
+    "\n",
+    "# fit parameters are not set yet\n",
+    "(explin.params, explin.param_errors)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each fit parameter is given a name."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "explin.param_names"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "These names are used throughout the [fit function's documentation](../../api/plasmapy.analysis.fit_functions.ExponentialPlusLinear.rst), as well as in its `__repr__`, `__str__`, and `latex_str` methods."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(explin, explin.__str__(), explin.latex_str)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fitting to data\n",
+    "\n",
+    "[curve_fit()]: ../../api/plasmapy.analysis.fit_functions.ExponentialPlusLinear.rst#plasmapy.analysis.fit_functions.ExponentialPlusLinear.curve_fit\n",
+    "[linregress()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html\n",
+    "[Linear]: ../../api/plasmapy.analysis.fit_functions.Linear.rst#plasmapy.analysis.fit_functions.Linear.curve_fit\n",
+    "\n",
+    "Fit functions provide the [curve_fit()] method to fit the analytical function to a set of $(x, y)$ data.  This is typically done with SciPy's [curve_fit()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) function, but fitting is done with SciPy's [linregress()] for the [Linear] fit funciton."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's generate some noisy data to fit to..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = (5.0, 0.1, -0.5, -8.0)  # (a, alpha, m, b)\n",
+    "xdata = np.linspace(-20, 15, num=100)\n",
+    "ydata = (explin.func(xdata, *params)\n",
+    "        + np.random.normal(0.0, 0.6, xdata.size))\n",
+    "\n",
+    "plt.plot(xdata, ydata)\n",
+    "plt.xlabel(\"X\", fontsize=14)\n",
+    "plt.ylabel(\"Y\", fontsize=14)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[curve_fit()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html\n",
+    "\n",
+    "The fit function `curve_fit()` shares the same signature as SciPy's [curve_fit()], so any `**kwargs` will be passed on.  By default, only the $(x, y)$ values are needed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "explin.curve_fit(xdata, ydata)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Getting fit results\n",
+    "\n",
+    "After fitting, the fitted parameters, uncertainties, and [coefficient of determination](https://en.wikipedia.org/wiki/Coefficient_of_determination), or $r^2$, values can be retrieved through their respective properties, `params`, `parame_errors`, and `rsq`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(explin.params, explin.params.a, explin.params.alpha)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "(explin.param_errors, explin.param_errors.a, explin.param_errors.alpha)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "explin.rsq"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Fit function is callable\n",
+    "\n",
+    "Now that parameters are set, the fit function is callable."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "explin(0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Associated errors can also be generated."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y, y_err = explin(np.linspace(-1, 1, num=10), reterr=True)\n",
+    "(y, y_err)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Known uncertainties in $x$ can be specified too."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y, y_err = explin(np.linspace(-1, 1, num=10), reterr=True, x_err=0.1)\n",
+    "(y, y_err)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Plotting results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+     "tooltip": "Fit Function Tutorial"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# plot original data\n",
+    "plt.plot(xdata, ydata, marker=\"o\", linestyle=\" \", label=\"Data\")\n",
+    "ax = plt.gca()\n",
+    "ax.set_xlabel(\"X\", fontsize=14)\n",
+    "ax.set_ylabel(\"Y\", fontsize=14)\n",
+    "\n",
+    "ax.axhline(0.0, color='r', linestyle='--')\n",
+    "\n",
+    "# plot fitted curve + error\n",
+    "yfit, yfit_err = explin(xdata, reterr=True)\n",
+    "ax.plot(xdata, yfit, color=\"orange\", label=\"Fit\")\n",
+    "ax.fill_between(xdata, yfit + yfit_err, yfit - yfit_err,\n",
+    "                 color=\"orange\", alpha=0.12, zorder=0,\n",
+    "                 label=\"Fit Error\")\n",
+    "\n",
+    "# plot annotations\n",
+    "plt.legend(fontsize=14, loc=\"upper left\")\n",
+    "\n",
+    "txt = (\n",
+    "    f\"$f(x) = {explin.latex_str}$\\n\"\n",
+    "    f\"$r^2 = {explin.rsq:.3f}$\\n\"\n",
+    ")\n",
+    "for name, param, err in zip(explin.param_names, explin.params, explin.param_errors):\n",
+    "    txt += f\"{name} = {param:.3f} $\\\\pm$ {err:.3f}\\n\"\n",
+    "txt_loc = [-13.0, ax.get_ylim()[1]]\n",
+    "txt_loc = ax.transAxes.inverted().transform(ax.transData.transform(txt_loc))\n",
+    "txt_loc[0] -= 0.02\n",
+    "txt_loc[1] -= 0.05\n",
+    "ax.text(txt_loc[0], txt_loc[1], txt, \n",
+    "        fontsize='large', transform=ax.transAxes,\n",
+    "        va=\"top\", linespacing=1.5)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Root solving\n",
+    "\n",
+    "[fsolve()]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html\n",
+    "[ExponentialPlusLinear().root_solve()]: ../../api/plasmapy.analysis.fit_functions.ExponentialPlusLinear.rst#plasmapy.analysis.fit_functions.ExponentialPlusLinear.root_solve\n",
+    "[Linear().root_solve()]: ../../api/plasmapy.analysis.fit_functions.Linear.rst#plasmapy.analysis.fit_functions.Linear.root_solve\n",
+    "\n",
+    "An exponential plus a linear offset has no analytical solutions for its roots, except for a few specific cases.  To get around this, [ExponentialPlusLinear().root_solve()] uses SciPy's [fsolve()] to calculate it's roots.  If a fit function has analytical solutions to its roots (e.g. [Linear().root_solve()]), then the method is overriden with the known solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "root, err = explin.root_solve(-15.)\n",
+    "(root, err)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[Linear().root_solve()]: ../../api/plasmapy.analysis.fit_functions.Linear.rst#plasmapy.analysis.fit_functions.Linear.root_solve\n",
+    "\n",
+    "Let's use [Linear().root_solve()] as an example for a known solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lin = ffuncs.Linear(params=(1.0, -5.0), param_errors=(0.1, 0.1))\n",
+    "root, err = lin.root_solve()\n",
+    "(root, err)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/plasmapy/analysis/__init__.py b/plasmapy/analysis/__init__.py
new file mode 100644
index 00000000..92c05eb8
--- /dev/null
+++ b/plasmapy/analysis/__init__.py
@@ -0,0 +1,4 @@
+"""The analysis subpackage for PlasmaPy."""
+__all__ = []
+
+from plasmapy.analysis import fit_functions
diff --git a/plasmapy/analysis/fit_functions.py b/plasmapy/analysis/fit_functions.py
new file mode 100644
index 00000000..89fc7718
--- /dev/null
+++ b/plasmapy/analysis/fit_functions.py
@@ -0,0 +1,1056 @@
+"""
+`FitFunction` classes designed to assist in curve fitting of swept Langmuir
+traces.
+"""
+__all__ = [
+    "AbstractFitFunction",
+    "Exponential",
+    "ExponentialPlusLinear",
+    "ExponentialPlusOffset",
+    "Linear",
+]
+
+import numpy as np
+
+from abc import ABC, abstractmethod
+from collections import namedtuple
+from scipy.optimize import curve_fit, fsolve
+from scipy.stats import linregress
+from typing import Tuple, Union
+
+from plasmapy.utils.decorators import modify_docstring
+
+#: Named tuple for :meth:`AbstractFitFunction.root_solve`.
+_RootResults = namedtuple("RootResults", ("root", "err"))
+
+
+class AbstractFitFunction(ABC):
+    """
+    Abstract class for defining fit functions :math:`f(x)` and the tools for
+    fitting the function to a set of data.
+    """
+
+    _param_names = NotImplemented  # type: Tuple[str, ...]
+
+    def __init__(
+        self, params: Tuple[float, ...] = None, param_errors: Tuple[float, ...] = None,
+    ):
+        """
+        Parameters
+        ----------
+        params: Tuple[float, ...], optional
+            Tuple of values for the function parameters. Equal in size to
+            :attr:`param_names`.
+
+        param_errors: Tuple[float, ...], optional
+            Tuple of values for the errors associated with the function
+            parameters.  Equal in size to :attr:`param_names`.
+
+        """
+
+        self.FitParamTuple = namedtuple("FitParamTuple", self._param_names)
+        """
+        A `~collections.namedtuple` used for attributes :attr:`params` and
+        :attr:`param_errors`.  The attribute :attr:`parameter_names` defines
+        the tuple field names.
+        """
+
+        if params is None:
+            self._params = None
+        else:
+            self.params = params
+
+        if param_errors is None:
+            self._param_errors = None
+        else:
+            self.param_errors = param_errors
+
+        self._curve_fit_results = None
+        self._rsq = None
+
+    def __call__(self, x, x_err=None, reterr=False):
+        """
+        Direct call of the fit function :math:`f(x)``.
+
+        Parameters
+        ----------
+        x: array_like
+            Dependent variables.
+
+        x_err: array_like, optional
+            Errors associated with the independent variables ``x``.  Must be of
+            size one or equal to the size of ``x``.
+
+        reterr: bool, optional
+            (Default: `False`) If `True`, return an array of uncertainties
+            associated with the calculated independent variables
+
+        Returns
+        -------
+        y: `numpy.ndarray`
+            Corresponding dependent variables :math:`y=f(x)` of the independent
+            variables ``x``.
+
+        y_err: `numpy.ndarray`
+            Uncertainties associated with the calculated dependent variables
+            :math:`\\delta y`
+        """
+        if reterr:
+            y_err, y = self.func_err(x, x_err=x_err, rety=True)
+
+            return y, y_err
+
+        y = self.func(x, *self.params)
+
+        return y
+
+    def __repr__(self):
+        return f"{self.__str__()} {self.__class__}"
+
+    @abstractmethod
+    def __str__(self):
+        ...
+
+    @abstractmethod
+    def func(self, x, *args):
+        """
+        The fit function.  This signature of the function must first take the
+        independent variable followed by the parameters to be fitted as
+        separate arguments.
+
+        Parameters
+        ----------
+        x: array_like
+            Independent variables to be passed to the fit function.
+
+        *args: Tuple[Union[float, int],...]
+            The parameters that will be adjusted to make the fit.
+
+        Returns
+        -------
+        `numpy.ndarray`:
+            The calculated dependent variables of the independent variables ``x``.
+
+        Notes
+        -----
+        * When sub-classing the definition should look something like::
+
+            def func(self, x, a, b, c):
+                x = self._check_x(x)
+                self._check_params(a, b, c)
+
+                return a * x ** 2 + b * x + c
+        """
+        ...
+
+    @abstractmethod
+    @modify_docstring(
+        prepend="""
+        Calculate dependent variable uncertainties :math:`\\delta y` for
+        dependent variables :math:`y=f(x)`.
+        """,
+        append="""
+        * When sub-classing the definition should look something like::
+
+            @modify_docstring(append=AbstractFitFunction.func_err.__original_doc__)
+            def func_err(self, x, x_err=None, rety=False):
+                '''
+                A simple docstring giving the equation for error propagation, but
+                excluding the parameter descriptions.  The @modify_docstring
+                decorator will append the docstring from the parent class.
+                '''
+                x, x_err = self._check_func_err_params(x, x_err)
+
+                a, b, c = self.params
+                a_err, b_err, c_err = self.param_errors
+
+                # calculate error
+
+                if rety:
+                    y = self.func(x, a, b, c)
+                    return err, y
+
+                return err
+        """,
+    )
+    def func_err(self, x, x_err=None, rety=False):
+        """
+        Parameters
+        ----------
+        x: array_like
+            Independent variables to be passed to the fit function.
+
+        x_err: array_like, optional
+            Errors associated with the independent variables ``x``.  Must be of
+            size one or equal to the size of ``x``.
+
+        rety: bool
+            Set `True` to also return the associated dependent variables
+            :math:`y = f(x)`.
+
+        Returns
+        -------
+        err: `numpy.ndarray`
+            The calculated uncertainties :math:`\\delta y` of the dependent
+            variables (:math:`y = f(x)`) of the independent variables ``x``.
+
+        y: `numpy.ndarray`, optional
+            (if ``rety == True``) The associated dependent variables
+            :math:`y = f(x)`.
+
+        Notes
+        -----
+        * A good reference for formulating propagation of uncertainty expressions is:
+
+            J. R. Taylor.  *An Introduction to Error Analysis: The Study of
+            Uncertainties in Physical Measurements.* University Science Books,
+            second edition, August 1996 (ISBN: 093570275X)
+
+        """
+        ...
+
+    @property
+    def curve_fit_results(self):
+        """
+        The results returned by the curve fitting routine used by
+        :attr:`curve_fit`.  This is typically from `scipy.stats.linregress` or
+        `scipy.optimize.curve_fit`.
+        """
+        return self._curve_fit_results
+
+    @property
+    def params(self) -> Union[None, tuple]:
+        """The fitted parameters for the fit function."""
+        if self._params is None:
+            return self._params
+        else:
+            return self.FitParamTuple(*self._params)
+
+    @params.setter
+    def params(self, val) -> None:
+        if isinstance(val, self.FitParamTuple):
+            self._params = tuple(val)
+        elif (
+            isinstance(val, (tuple, list))
+            and len(val) == len(self.param_names)
+            and all(isinstance(vv, (int, np.integer, float, np.floating)) for vv in val)
+        ):
+            self._params = tuple(val)
+        else:
+            raise ValueError(
+                f"Got {val} for 'val', expecting tuple of ints and "
+                f"floats of length {len(self.param_names)}."
+            )
+
+    @property
+    def param_errors(self) -> Union[None, tuple]:
+        """The associated errors of the fitted :attr:`params`."""
+        if self._param_errors is None:
+            return self._param_errors
+        else:
+            return self.FitParamTuple(*self._param_errors)
+
+    @param_errors.setter
+    def param_errors(self, val) -> None:
+        if isinstance(val, self.FitParamTuple):
+            self._param_errors = tuple(val)
+        elif (
+            isinstance(val, (tuple, list))
+            and len(val) == len(self.param_names)
+            and all(isinstance(vv, (int, np.integer, float, np.floating)) for vv in val)
+        ):
+            self._param_errors = tuple(val)
+        else:
+            raise ValueError(
+                f"Got {val} for 'val', expecting tuple of ints and "
+                f"floats of length {len(self.param_names)}."
+            )
+
+    @property
+    def param_names(self) -> Tuple[str, ...]:
+        """Names of the fitted parameters."""
+        return self._param_names
+
+    @property
+    @abstractmethod
+    def latex_str(self) -> str:
+        """LaTeX friendly representation of the fit function."""
+        ...
+
+    def _check_func_err_params(self, x, x_err):
+        """Check the ``x`` and ``x_err`` parameters for :meth:`func_err`."""
+        x = self._check_x(x)
+        if x_err is not None:
+            x_err = self._check_x(x_err)
+
+            if x_err.shape == ():
+                pass
+            elif x_err.shape != x.shape:
+                raise ValueError(
+                    f"x_err shape {x_err.shape} must be equal the shape of "
+                    f"x {x.shape}."
+                )
+        return x, x_err
+
+    @staticmethod
+    def _check_params(*args) -> None:
+        """
+        Check fitting parameters so that they are an expected type for the
+        class functionality.
+        """
+        for arg in args:
+            if not isinstance(arg, (int, np.integer, float, np.floating)):
+                raise TypeError(
+                    f"Expected int or float for parameter argument, got "
+                    f"{type(arg)}."
+                )
+
+    @staticmethod
+    def _check_x(x):
+        """
+        Check the independent variable ``x`` so that it is an expected
+        type for the class functionality.
+        """
+        if isinstance(x, (int, float, np.integer, np.floating)):
+            x = np.array(x)
+        else:
+            if not isinstance(x, np.ndarray):
+                x = np.array(x)
+
+            if not (
+                np.issubdtype(x.dtype, np.integer)
+                or np.issubdtype(x.dtype, np.floating)
+            ):
+                raise TypeError(
+                    f"Argument x needs to be an array_like object of integers "
+                    f"or floats."
+                )
+
+            x = x.squeeze()
+            if x.shape == ():
+                # force x to be a scalar
+                x = x[()]
+
+        return x
+
+    def root_solve(self, x0):
+        """
+        Solve for the root of the fit function (i.e. :math:`f(x_r) = 0`).  This
+        mehtod used `scipy.optimize.fsolve` to find the function roots.
+
+        Parameters
+        ----------
+        x0: `~numpy.ndarray`
+            The starting estimate for the roots of :math:`f(x_r) = 0`.
+
+        Returns
+        -------
+        x : `~numpy.ndarray`
+            The solution (or the result of the last iteration for an
+            unsuccessful call).
+
+        x_err: `~numpy.ndarray`
+            The uncertainty associated with the root calculation.  **Currently
+            this returns an array of** `numpy.nan` **values equal in shape to**
+            ``x`` **, since there is no determined way to calculate the
+            uncertainties.**
+
+        Notes
+        -----
+        If the full output of `scipy.optimize.fsolve` is desired then one can do:
+
+            >>> func = Linear()
+            >>> func.params = (1.0, 5.0)
+            >>> func.param_errors = (0.0, 0.0)
+            >>> roots = fsolve(func, -4.0, full_output=True)
+            >>> roots
+            (array([-5.]),
+             {'nfev': 4,
+              'fjac': array([[-1.]]),
+              'r': array([-1.]),
+              'qtf': array([2.18...e-12]),
+              'fvec': 0.0},
+             1,
+             'The solution converged.')
+
+        """
+        results = fsolve(self.func, x0, args=self.params)
+
+        root = np.squeeze(results[0])
+        err = np.tile(np.nan, root.shape)
+        if root.shape == ():
+            # force x to be a scalar
+            root = root[()]
+            err = np.nan
+
+        return _RootResults(root, err)
+
+    @property
+    def rsq(self):
+        """
+        Coefficient of determination (r-squared) value of the fit.
+
+        .. math::
+
+            r^2 &= 1 - \\frac{SS_{res}}{SS_{tot}}
+
+            SS_{res} &= \\sum\\limits_{i} (y_i - f(x_i))^2
+
+            SS_{tot} &= \\sum\\limits_{i} (y_i - \\bar{y})^2
+
+        where :math:`(x_i, y_i)` are the sample data pairs, :math:`f(x_i)` is
+        the fitted dependent variable corresponding to :math:`x_i`, and
+        :math:`\\bar{y}` is the average of the :math:`y_i` values.
+
+        The :math:`r^2` value is an indicator of how close the points
+        :math:`(x_i, y_i)` lie to the model :math:`f(x)`.  :math:`r^2` values
+        range between 0 and 1.  Values close to 0 indicate that the points
+        are uncorrelated and have little tendency to lie close to the model,
+        whereas, values close to 1 indicate a high correlation to the model.
+
+        """
+        return self._rsq
+
+    def curve_fit(self, xdata, ydata, **kwargs) -> None:
+        """
+        Use a non-linear least squares method to fit the fit function to
+        (``xdata``, ``ydata``), using `scipy.optimize.curve_fit`.  This will set
+        the attributes :attr:`parameters`, :attr:`parameters_err`, and
+        :attr:`rsq`.
+
+        The results of `scipy.optimize.curve_fit` can be obtained via
+        :attr:`curve_fit_results`.
+
+        Parameters
+        ----------
+        xdata: array_like
+            The independent variable where data is measured.  Should be 1D of
+            length M.
+
+        ydata: array_like
+            The dependent data associated with ``xdata``.
+
+        **kwargs
+            Any keywords accepted by `scipy.optimize.curve_fit`.
+
+        Raises
+        ------
+        ValueError
+            if either ``ydata`` or ``xdata`` contain `numpy.nan`'s, or if
+            incompatible options are used.
+
+        RuntimeError
+            if the least-squares minimization fails.
+
+        ~scipy.optimize.OptimizeWarning
+            if covariance of the parameters can not be estimated.
+
+        """
+        popt, pcov = curve_fit(self.func, xdata, ydata, **kwargs)
+        self._curve_fit_results = (popt, pcov)
+        self.params = tuple(popt.tolist())
+        self.param_errors = tuple(np.sqrt(np.diag(pcov)).tolist())
+
+        # calc rsq
+        # rsq = 1 - (ss_res / ss_tot)
+        residuals = ydata - self.func(xdata, *self.params)
+        ss_res = np.sum(residuals ** 2)
+        ss_tot = np.sum((ydata - np.mean(ydata)) ** 2)
+        self._rsq = 1 - (ss_res / ss_tot)
+
+
+class Linear(AbstractFitFunction):
+    """
+    A sub-class of `AbstractFitFunction` to represent a linear function.
+
+    .. math::
+
+        y &= f(x) = m \\, x + b
+
+        (\\delta y)^2 &= (x \\, \\delta m)^2 + (m \\, \\delta x)^2 + (\\delta b)^2
+
+    where :math:`m` and :math:`b` are real constants to be fitted and :math:`x` is
+    the independent variable.  :math:`\\delta m`, :math:`\\delta b`, and
+    :math:`\\delta x` are the respective uncertainties for :math:`m`, :math:`b`,
+    and :math:`x`.
+    """
+
+    _param_names = ("m", "b")
+
+    def __str__(self):
+        return f"f(x) = m x + b"
+
+    @property
+    def latex_str(self) -> str:
+        return fr"m x + b"
+
+    def func(self, x, m, b):
+        """
+        The fit function, a linear function.
+
+        .. math::
+
+            f(x) = m \\, x + b
+
+        where :math:`m` and :math:`b` are real constants representing the
+        slope and intercept, respectively, and :math:`x` is the independent
+        variable.
+
+        Parameters
+        ----------
+        x: array_like
+            Independent variable.
+
+        m: float
+            value for slope :math:`m`
+
+        b: float
+            value for intercept :math:`b`
+
+        Returns
+        -------
+        y: array_like
+            dependent variables corresponding to `:math:``x``
+
+        """
+        x = self._check_x(x)
+        self._check_params(m, b)
+
+        return m * x + b
+
+    @modify_docstring(append=AbstractFitFunction.func_err.__original_doc__)
+    def func_err(self, x, x_err=None, rety=False):
+        """
+        Calculate dependent variable uncertainties :math:`\\delta y` for
+        dependent variables :math:`y=f(x)`.
+
+        .. math::
+
+            (\\delta y)^2 = (x \\, \\delta m)^2 + (m \\, \\delta x)^2 + (\\delta b)^2
+
+        """
+        x, x_err = self._check_func_err_params(x, x_err)
+
+        m, b = self.params
+        m_err, b_err = self.param_errors
+
+        m_term = (m_err * x) ** 2
+        b_term = b_err ** 2
+        err = m_term + b_term
+
+        if x_err is not None:
+            x_term = (m * x_err) ** 2
+            err += x_term
+        err = np.sqrt(err)
+
+        if rety:
+            y = self.func(x, m, b)
+            return err, y
+
+        return err
+
+    @property
+    def rsq(self):
+        """
+        Coefficient of determination (r-squared) value of the fit.  Calculated
+        by `scipy.stats.linregress` from the fit.
+        """
+        return self._rsq
+
+    def root_solve(self, *args, **kwargs):
+        """
+        The root :math:`f(x_r) = 0` for the fit function.
+
+        .. math::
+
+            x_r &= \\frac{-b}{m}
+
+            \\delta x_r &= |x_r| \\sqrt{
+                \\left( \\frac{\\delta m}{m} \\right)^2
+                + \\left( \\frac{\\delta b}{b} \\right)^2
+            }
+
+        Parameters
+        ----------
+        *args
+            Not needed.  This is to ensure signature comparability with
+            `AbstractFitFunction`.
+
+        *kwargs
+            Not needed.  This is to ensure signature comparability with
+            `AbstractFitFunction`.
+
+        Returns
+        -------
+        root: float
+            The root value for the given fit :attr:`parameters`.
+
+        err: float
+            The uncertainty in the calculated root for the given fit
+            :attr:`parameters` and :attr:`parameters_err`.
+        """
+        m, b = self.params
+        root = -b / m
+
+        m_err, b_err = self.param_errors
+        err = np.abs(root) * np.sqrt((m_err / m) ** 2 + (b_err / b) ** 2)
+
+        return _RootResults(root, err)
+
+    def curve_fit(self, xdata, ydata, **kwargs) -> None:
+        """
+        Calculate a linear least-squares regression of (``xdata``, ``ydata``)
+        using `scipy.stats.linregress`.  This will set the attributes
+        :attr:`parameters`, :attr:`parameters_err`, and :attr:`rsq`.
+
+        The results of `scipy.stats.linregress` can be obtained via
+        :attr:`curve_fit_results`.
+
+        Parameters
+        ----------
+        xdata: array_like
+            The independent variable where data is measured.  Should be 1D of
+            length M.
+
+        ydata: array_like
+            The dependent data associated with ``xdata``.
+
+        **kwargs
+            Any keywords accepted by `scipy.stats.linregress.curve_fit`.
+
+        """
+        results = linregress(xdata, ydata)
+        self._curve_fit_results = results
+
+        m = results[0]
+        b = results[1]
+        self.params = (m, b)
+
+        m_err = results[4]
+        b_err = np.sum(xdata ** 2) - ((np.sum(xdata) ** 2) / xdata.size)
+        b_err = m_err * np.sqrt(1.0 / b_err)
+        self.param_errors = (m_err, b_err)
+
+        self._rsq = results[2] ** 2
+
+
+class Exponential(AbstractFitFunction):
+    """
+    A sub-class of `AbstractFitFunction` to represent an exponential with an
+    offset.
+
+    .. math::
+
+        y &= f(x) = a \\, e^{\\alpha \\, x}
+
+        \\left( \\frac{\\delta y}{|y|} \\right)^2 &=
+            \\left( \\frac{\\delta a}{a} \\right)^2
+            + (x \\, \\delta \\alpha)^2
+            + (\\alpha \\, \\delta x)^2
+
+    where :math:`a` and :math:`\\alpha` are the real constants to be fitted and
+    :math:`x` is the independent variable.  :math:`\\delta a`,
+    :math:`\\delta \\alpha`, and :math:`\\delta x` are the respective
+    uncertainties for :math:`a`, :math:`\\alpha`, and :math:`x`.
+    """
+
+    _param_names = ("a", "alpha")
+
+    def __str__(self):
+        return f"f(x) = a exp(alpha x)"
+
+    @property
+    def latex_str(self) -> str:
+        return fr"a \, \exp(\alpha x)"
+
+    def func(self, x, a, alpha):
+        """
+        The fit function, a exponential function.
+
+        .. math::
+
+            f(x) = a \\, e^{\\alpha \\, x}
+
+        where :math:`a` and :math:`\\alpha` are real constants and :math:`x`
+        is the independent variable.
+
+        Parameters
+        ----------
+        x: array_like
+            Independent variable.
+
+        a: float
+            value for the exponential "normalization" constant, :math:`a`
+
+        alpha: float
+            value for the growth constant, :math:`\\alpha`
+
+        Returns
+        -------
+        y: array_like
+            dependent variables corresponding to ``x``
+
+        """
+        x = self._check_x(x)
+        self._check_params(a, alpha)
+
+        return a * np.exp(alpha * x)
+
+    @modify_docstring(append=AbstractFitFunction.func_err.__original_doc__)
+    def func_err(self, x, x_err=None, rety=False):
+        """
+        Calculate dependent variable uncertainties :math:`\\delta y` for
+        dependent variables :math:`y=f(x)`.
+
+        .. math::
+
+            \\left( \\frac{\\delta y}{|y|} \\right)^2 =
+                \\left( \\frac{\\delta a}{a} \\right)^2
+                + (x \\, \\delta \\alpha)^2
+                + (\\alpha \\, \\delta x)^2
+
+        """
+        x, x_err = self._check_func_err_params(x, x_err)
+
+        a, alpha = self.params
+        a_err, alpha_err = self.param_errors
+        y = self.func(x, a, alpha)
+
+        a_term = (a_err / a) ** 2
+        alpha_term = (x * alpha_err) ** 2
+
+        err = a_term + alpha_term
+
+        if x_err is not None:
+            x_term = (alpha * x_err) ** 2
+            err += x_term
+
+        err = np.abs(y) * np.sqrt(err)
+
+        if rety:
+            return err, y
+
+        return err
+
+    def root_solve(self, *args, **kwargs):
+        """
+        The root :math:`f(x_r) = 0` for the fit function. **An exponential has no
+        real roots.**
+
+        Parameters
+        ----------
+        *args
+            Not needed.  This is to ensure signature compatibility with
+            `AbstractFitFunction`.
+
+        *kwargs
+            Not needed.  This is to ensure signature compatibility with
+            `AbstractFitFunction`.
+
+        Returns
+        -------
+        root: float
+            The root value for the given fit :attr:`parameters`.
+
+        err: float
+            The uncertainty in the calculated root for the given fit
+            :attr:`parameters` and :attr:`parameters_err`.
+        """
+
+        return _RootResults(np.nan, np.nan)
+
+
+class ExponentialPlusLinear(AbstractFitFunction):
+    """
+    A sub-class of `AbstractFitFunction` to represent an exponential with an
+    linear offset.
+
+    .. math::
+
+        y =& f(x) = a \\, e^{\\alpha \\, x} + m \\, x + b\\\\
+        (\\delta y)^2 =&
+            \\left( a e^{\\alpha x}\\right)^2 \\left[
+                \\left( \\frac{\\delta a}{a} \\right)^2
+                + (x \\, \\delta \\alpha)^2
+                + (\\alpha \\, \\delta x)^2
+            \\right]\\\\
+            & + \\left(2 \\, a \\, \\alpha \\, m \\, e^{\\alpha x}\\right)
+                (\\delta x)^2\\\\
+            & + \\left[(x \\, \\delta m)^2 + (\\delta b)^2 +(m \\, \\delta x)^2\\right]
+
+    where :math:`a`, :math:`\\alpha`, :math:`m`, and :math:`b` are the real
+    constants to be fitted and :math:`x` is the independent variable.
+    :math:`\\delta a`, :math:`\\delta \\alpha`, :math:`\\delta m`, :math:`\\delta b`,
+    and :math:`\\delta x` are the respective uncertainties for :math:`a`,
+    :math:`\\alpha`, :math:`m`, and :math:`b`, and :math:`x`.
+    """
+
+    _param_names = ("a", "alpha", "m", "b")
+
+    def __init__(
+        self, params: Tuple[float, ...] = None, param_errors: Tuple[float, ...] = None,
+    ):
+        self._exponential = Exponential()
+        self._linear = Linear()
+        super().__init__(params=params, param_errors=param_errors)
+
+    def __str__(self):
+        exp_str = self._exponential.__str__().lstrip("f(x) = ")
+        lin_str = self._linear.__str__().lstrip("f(x) = ")
+        return f"f(x) = {exp_str} + {lin_str}"
+
+    @property
+    def latex_str(self) -> str:
+        exp_str = self._exponential.latex_str
+        lin_str = self._linear.latex_str
+        return fr"{exp_str} + {lin_str}"
+
+    @AbstractFitFunction.params.setter
+    def params(self, val) -> None:
+        AbstractFitFunction.params.fset(self, val)
+        self._exponential.params = (self.params.a, self.params.alpha)
+        self._linear.params = (self.params.m, self.params.b)
+
+    @AbstractFitFunction.param_errors.setter
+    def param_errors(self, val) -> None:
+        AbstractFitFunction.param_errors.fset(self, val)
+        self._exponential.param_errors = (
+            self.param_errors.a,
+            self.param_errors.alpha,
+        )
+        self._linear.param_errors = (self.param_errors.m, self.param_errors.b)
+
+    def func(self, x, a, alpha, m, b):
+        """
+        The fit function, an exponential with a linear offset.
+
+        .. math::
+
+            f(x) = a \\, e^{\\alpha \\, x} + m \\, x + b\\\\
+
+        where :math:`a`, :math:`\\alpha`, :math:`m`, and :math:`b` are the real
+        constants and :math:`x` is the independent variable.
+
+        Parameters
+        ----------
+        x: array_like
+            Independent variable.
+
+        a: float
+            value for constant :math:`a`
+
+        alpha: float
+            value for constant :math:`\\alpha`
+
+        m: float
+            value for slope :math:`m`
+
+        b: float
+            value for intercept :math:`b`
+
+        Returns
+        -------
+        y: array_like
+            dependent variables corresponding to ``x``
+
+        """
+        exp_term = self._exponential.func(x, a, alpha)
+        lin_term = self._linear.func(x, m, b)
+        return exp_term + lin_term
+
+    @modify_docstring(append=AbstractFitFunction.func_err.__original_doc__)
+    def func_err(self, x, x_err=None, rety=False):
+        """
+        Calculate dependent variable uncertainties :math:`\\delta y` for
+        dependent variables :math:`y=f(x)`.
+
+        .. math::
+
+            (\\delta y)^2 =&
+                \\left( a e^{\\alpha x}\\right)^2 \\left[
+                    \\left( \\frac{\\delta a}{a} \\right)^2
+                    + (x \\, \\delta \\alpha)^2
+                    + (\\alpha \\, \\delta x)^2
+                \\right]\\\\
+                & + \\left(2 \\, a \\, \\alpha \\, m \\, e^{\\alpha x}\\right)
+                    (\\delta x)^2\\\\
+                & + \\left[(
+                        x \\, \\delta m)^2 + (\\delta b)^2 +(m \\, \\delta x)^2
+                    \\right]
+
+        """
+        x, x_err = self._check_func_err_params(x, x_err)
+
+        a, alpha, m, b = self.params
+
+        exp_y, exp_err = self._exponential(x, x_err=x_err, reterr=True)
+        lin_y, lin_err = self._linear(x, x_err=x_err, reterr=True)
+        err = exp_err ** 2 + lin_err ** 2
+
+        if x_err is not None:
+            blend_err = 2 * a * alpha * m * np.exp(alpha * x) * (x_err ** 2)
+            err += blend_err
+        err = np.sqrt(err)
+
+        if rety:
+            return err, exp_y + lin_y
+
+        return err
+
+
+class ExponentialPlusOffset(AbstractFitFunction):
+    """
+    A sub-class of `AbstractFitFunction` to represent an exponential with a
+    constant offset.
+
+    .. math::
+
+        y =& f(x) = a \\, e^{\\alpha \\, x} + m \\, x + b\\\\
+        (\\delta y)^2 =&
+            \\left( a e^{\\alpha x}\\right)^2 \\left[
+                \\left( \\frac{\\delta a}{a} \\right)^2
+                + (x \\, \\delta \\alpha)^2
+                + (\\alpha \\, \\delta x)^2
+            \\right]
+            + (\\delta b)^2
+
+
+    where :math:`a`, :math:`\\alpha`, and :math:`b` are the real constants to
+    be fitted and :math:`x` is the independent variable.  :math:`\\delta a`,
+    :math:`\\delta \\alpha`, :math:`\\delta b`, and :math:`\\delta x` are the
+    respective uncertainties for :math:`a`, :math:`\\alpha`, and :math:`b`, and
+    :math:`x`.
+
+    """
+
+    _param_names = ("a", "alpha", "b")
+
+    def __init__(
+        self, params: Tuple[float, ...] = None, param_errors: Tuple[float, ...] = None,
+    ):
+        self._explin = ExponentialPlusLinear()
+        super().__init__(params=params, param_errors=param_errors)
+
+    def __str__(self):
+        return f"f(x) = a exp(alpha x) + b"
+
+    @property
+    def latex_str(self) -> str:
+        return fr"a \, \exp(\alpha x) + b"
+
+    @AbstractFitFunction.params.setter
+    def params(self, val) -> None:
+        AbstractFitFunction.params.fset(self, val)
+        self._explin.params = (
+            self.params.a,
+            self.params.alpha,
+            0.0,
+            self.params.b,
+        )
+
+    @AbstractFitFunction.param_errors.setter
+    def param_errors(self, val) -> None:
+        AbstractFitFunction.param_errors.fset(self, val)
+        self._explin.param_errors = (
+            self.param_errors.a,
+            self.param_errors.alpha,
+            0.0,
+            self.param_errors.b,
+        )
+
+    def func(self, x, a, alpha, b):
+        """
+        The fit function, an exponential with a constant offset.
+
+        .. math::
+
+            f(x) = a \\, e^{\\alpha \\, x} + b\\\\
+
+        where :math:`a`, :math:`\\alpha`, and :math:`b` are the real constants
+        and :math:`x` is the independent variable.
+
+        Parameters
+        ----------
+        x: array_like
+            Independent variable.
+
+        a: float
+            value for constant :math:`a`
+
+        alpha: float
+            value for constant :math:`\\alpha`
+
+        b: float
+            value for DC offset :math:`b`
+
+        Returns
+        -------
+        y: array_like
+            dependent variables corresponding to ``x``
+
+        """
+        return self._explin.func(x, a, alpha, 0.0, b)
+
+    @modify_docstring(append=AbstractFitFunction.func_err.__original_doc__)
+    def func_err(self, x, x_err=None, rety=False):
+        """
+        Calculate dependent variable uncertainties :math:`\\delta y` for
+        dependent variables :math:`y=f(x)`.
+
+        .. math::
+
+            (\\delta y)^2 =
+                \\left( a e^{\\alpha x}\\right)^2 \\left[
+                    \\left( \\frac{\\delta a}{a} \\right)^2
+                    + (x \\, \\delta \\alpha)^2
+                    + (\\alpha \\, \\delta x)^2
+                \\right]
+                + (\\delta b)^2
+
+        """
+        return self._explin.func_err(x, x_err=x_err, rety=rety)
+
+    def root_solve(self, *args, **kwargs):
+        """
+        The root :math:`f(x_r) = 0` for the fit function.
+
+        .. math::
+
+            x_r &= \\frac{1}{\\alpha} \\ln \\left( \\frac{-b}{a} \\right)
+
+            \\delta x_r &= \\sqrt{
+                \\left( \\frac{1}{\\alpha} \\frac{\\delta a}{a} \\right)^2
+                + \\left( x_r \\frac{\\delta \\alpha}{\\alpha} \\right)^2
+                + \\left( \\frac{1}{\\alpha} \\frac{\\delta b}{b} \\right)^2
+            }
+
+        Parameters
+        ----------
+        *args
+            Not needed.  This is to ensure signature compatibility with
+            `AbstractFitFunction`.
+
+        *kwargs
+            Not needed.  This is to ensure signature compability with
+            `AbstractFitFunction`.
+
+        Returns
+        -------
+        root: float
+            The root value for the given fit :attr:`parameters`.
+
+        err: float
+            The uncertainty in the calculated root for the given fit
+            :attr:`parameters` and :attr:`parameters_err`.
+        """
+        a, alpha, b = self.params
+        a_err, b_err, c_err = self.param_errors
+
+        root = np.log(-b / a) / alpha
+
+        a_term = a_err / (a * alpha)
+        b_term = b_err * root / alpha
+        c_term = c_err / (alpha * b)
+        err = np.sqrt(a_term ** 2 + b_term ** 2 + c_term ** 2)
+
+        return _RootResults(root, err)


Test patch:
diff --git a/plasmapy/analysis/tests/__init__.py b/plasmapy/analysis/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/plasmapy/analysis/tests/test_fit_functions.py b/plasmapy/analysis/tests/test_fit_functions.py
new file mode 100644
index 00000000..d25890f8
--- /dev/null
+++ b/plasmapy/analysis/tests/test_fit_functions.py
@@ -0,0 +1,587 @@
+"""
+Tests for the fitting function classes defined in `plasmapy.analysis.fit_functions`.
+"""
+import numpy as np
+import pytest
+
+from abc import ABC, abstractmethod
+from contextlib import ExitStack as does_not_raise
+# ^ ExitStack can be replaced with nullcontext when we require >= python 3.7
+
+import plasmapy.analysis.fit_functions as ffuncs
+
+
+class TestAbstractFitFunction:
+    """
+    Tests for fit function class `plasmapy.analysis.fit_functions.AbstractFitFunction`.
+
+    Notes
+    -----
+    Since `AbstractFitFunction` can not be instantiated, its complete
+    functionality can not be directly tested.  To resolve this, `BaseFFTests`,
+    which is the base test class for the fit function classes, will test all
+    the functionality within `AbstractFitFunction`.
+    """
+
+    ff_class = ffuncs.AbstractFitFunction
+
+    def test_is_abc(self):
+        """Test `AbstractFitFunction` is an abstract base class."""
+        assert issubclass(self.ff_class, ABC)
+
+    @pytest.mark.parametrize(
+        "name, isproperty",
+        [
+            ("__call__", False),
+            ("curve_fit", False),
+            ("curve_fit_results", True),
+            ("func", False),
+            ("func_err", False),
+            ("latex_str", True),
+            ("param_errors", True),
+            ("param_names", True),
+            ("params", True),
+            ("rsq", True),
+            ("root_solve", False),
+        ],
+    )
+    def test_methods(self, name, isproperty):
+        """Test for required methods and properties."""
+        assert hasattr(self.ff_class, name)
+
+        if isproperty:
+            assert isinstance(getattr(self.ff_class, name), property)
+
+    @pytest.mark.parametrize(
+        "name", ["__str__", "func", "func_err", "latex_str"],
+    )
+    def test_abstractmethods(self, name):
+        """Test for required abstract methods."""
+        assert name in self.ff_class.__abstractmethods__
+
+
+class BaseFFTests(ABC):
+    abc = ffuncs.AbstractFitFunction
+    _test_params = NotImplemented  # type: tuple
+    _test_param_errors = NotImplemented  # type: tuple
+    _test_param_names = NotImplemented  # type: tuple
+    _test_latex_str = NotImplemented  # type: str
+    _test__str__ = NotImplemented  # type: str
+
+    @property
+    @abstractmethod
+    def ff_class(self):
+        """Fit function class to be tested."""
+        ...
+
+    @staticmethod
+    @abstractmethod
+    def func(x, *args):
+        """
+        Formula/Function that the fit function class is suppose to be modeling.
+        This is used to test the fit function `func` method.
+        """
+        ...
+
+    @abstractmethod
+    def func_err(self, x, params, param_errors, x_err=None):
+        """
+        Function representing the propagation of error associated with fit function
+        model. This is used to test the fit function `func_err` method.
+        """
+        ...
+
+    def test_inheritance(self):
+        """Test inheritance from `AbstractFitFunction`."""
+        assert issubclass(self.ff_class, self.abc)
+
+    def test_iscallable(self):
+        """Test instantiated fit function is callable."""
+        assert callable(self.ff_class())
+
+    def test_repr(self):
+        """Test __repr__."""
+        ff_obj = self.ff_class()
+        assert ff_obj.__repr__() == f"{ff_obj.__str__()} {ff_obj.__class__}"
+
+    @pytest.mark.parametrize(
+        "name, isproperty",
+        [
+            ("__call__", False),
+            ("_param_names", False),
+            ("curve_fit", False),
+            ("curve_fit_results", True),
+            ("func", False),
+            ("func_err", False),
+            ("latex_str", True),
+            ("param_errors", True),
+            ("param_names", True),
+            ("params", True),
+            ("rsq", True),
+            ("root_solve", False),
+        ],
+    )
+    def test_methods(self, name, isproperty):
+        """Test attribute/method/property existence."""
+        assert hasattr(self.ff_class, name)
+
+        if isproperty:
+            assert isinstance(getattr(self.ff_class, name), property)
+
+        if name == "_param_names":
+            if self.ff_class._param_names == NotImplemented:
+                pytest.fail(
+                    f"{self.ff_class} class attribute '_param_names' needs to "
+                    f" be defined as a tuple of strings representing the names of "
+                    f"the fit parameters."
+                )
+
+    @pytest.mark.parametrize(
+        "name, value_ref_name",
+        [
+            ("param_names", "_test_param_names"),
+            ("latex_str", "_test_latex_str"),
+            ("__str__", "_test__str__"),
+        ],
+    )
+    def test_abstractmethod_values(self, name, value_ref_name):
+        """Test value of all abstract methods, except `func` and `func_err`."""
+        ff_obj = self.ff_class()
+
+        value = getattr(ff_obj, name)
+        if callable(value):
+            value = value()
+
+        exp_value = getattr(self, value_ref_name)
+        if exp_value == NotImplemented:
+            pytest.fail(
+                f"The expected value for abstract method {name} is not "
+                f"implemented/defined in the test class attribute {value_ref_name}."
+            )
+
+        assert value == exp_value
+
+    @pytest.mark.parametrize(
+        "params, param_errors, with_condition",
+        [
+            (None, None, does_not_raise()),
+            ("default", "default", does_not_raise()),
+            (5, None, pytest.raises(ValueError)),
+            (None, 5, pytest.raises(ValueError)),
+            (["wrong"], None, pytest.raises(ValueError)),
+            (None, ["wrong"], pytest.raises(ValueError)),
+            ("default+", None, pytest.raises(ValueError)),
+            (None, "default+", pytest.raises(ValueError)),
+        ],
+    )
+    def test_instantiation(self, params, param_errors, with_condition):
+        """Test behavior of fit function class instantiation."""
+        if params == "default":
+            params = self._test_params
+        elif params == "default+":
+            params = self._test_params
+            params = list(params)
+            params.append(5)
+
+        if param_errors == "default":
+            param_errors = self._test_param_errors
+        elif param_errors == "default+":
+            param_errors = self._test_param_errors
+            param_errors = list(param_errors)
+            param_errors.append(5)
+
+        with with_condition:
+            ff_obj = self.ff_class(params=params, param_errors=param_errors)
+
+            assert ff_obj.curve_fit_results is None
+            assert ff_obj.rsq is None
+
+            if params is None:
+                assert ff_obj.params is None
+            else:
+                assert ff_obj.params == ff_obj.FitParamTuple(*params)
+
+            if param_errors is None:
+                assert ff_obj.param_errors is None
+            else:
+                assert ff_obj.param_errors == ff_obj.FitParamTuple(*param_errors)
+
+    def test_param_namedtuple(self):
+        """
+        Test that the namedtuple used for `params` and `param_errors` is
+        constructed correctly.
+        """
+        ff_obj = self.ff_class()
+        assert hasattr(ff_obj, "FitParamTuple")
+        assert issubclass(ff_obj.FitParamTuple, tuple)
+        for name in ff_obj.param_names:
+            assert hasattr(ff_obj.FitParamTuple, name)
+
+    def test_param_names(self):
+        """Test attribute `param_names` is defined correctly."""
+        ff_obj = self.ff_class()
+        assert isinstance(ff_obj.param_names, tuple)
+        assert len(ff_obj.param_names) != 0
+        assert all(isinstance(val, str) for val in ff_obj.param_names)
+
+    @pytest.mark.parametrize(
+        "params, extra, with_condition",
+        [
+            ([2], None, does_not_raise()),
+            (5, None, pytest.raises(ValueError)),
+            (["wrong"], None, pytest.raises(ValueError)),
+            ([3], 10, pytest.raises(ValueError)),
+        ],
+    )
+    def test_params_setting(self, params, extra, with_condition):
+        """Tests for property setting of attribute `params`."""
+        ff_obj = self.ff_class()
+
+        if isinstance(params, list) and len(params) == 1:
+            params = params * len(ff_obj.param_names)
+        if extra is not None:
+            params.append(extra)
+
+        with with_condition:
+            ff_obj.params = params
+            assert ff_obj.params == ff_obj.FitParamTuple(*params)
+
+    @pytest.mark.parametrize(
+        "param_errors, extra, with_condition",
+        [
+            ([2], None, does_not_raise()),
+            (5, None, pytest.raises(ValueError)),
+            (["wrong"], None, pytest.raises(ValueError)),
+            ([3], 10, pytest.raises(ValueError)),
+        ],
+    )
+    def test_param_errors_setting(self, param_errors, extra, with_condition):
+        """Tests for property setting of attribute `param_errors`."""
+        ff_obj = self.ff_class()
+
+        if isinstance(param_errors, list) and len(param_errors) == 1:
+            param_errors = param_errors * len(ff_obj.param_names)
+        if extra is not None:
+            param_errors.append(extra)
+
+        with with_condition:
+            ff_obj.param_errors = param_errors
+            assert ff_obj.param_errors == ff_obj.FitParamTuple(*param_errors)
+
+    @pytest.mark.parametrize(
+        "x, replace_a_param, with_condition",
+        [
+            (0, None, does_not_raise()),
+            (1.0, None, does_not_raise()),
+            (np.linspace(10, 30, num=20), None, does_not_raise()),
+            ([4, 5, 6], None, does_not_raise()),
+            ("hello", None, pytest.raises(TypeError)),
+            (5, "hello", pytest.raises(TypeError)),
+        ],
+    )
+    def test_func(self, x, replace_a_param, with_condition):
+        """Test the `func` method."""
+        ff_obj = self.ff_class()
+
+        params = self._test_params
+        if replace_a_param is not None:
+            params = list(params)
+            params[0] = replace_a_param
+
+        with with_condition:
+            y = ff_obj.func(x, *params)
+
+            if isinstance(x, list):
+                x = np.array(x)
+            y_expected = self.func(x, *params)
+
+            assert np.allclose(y, y_expected)
+
+    @pytest.mark.parametrize(
+        "x, kwargs, with_condition",
+        [
+            (0, {}, does_not_raise()),
+            (1.0, {}, does_not_raise()),
+            (np.linspace(10, 30, num=20), {}, does_not_raise()),
+            ([4, 5, 6], {"x_err": 0.1, "rety": True}, does_not_raise()),
+            ("hello", {}, pytest.raises(TypeError)),
+            (5, {"x_err": "goodbye"}, pytest.raises(TypeError)),
+            (5, {"x_err": [0.1, 0.1]}, pytest.raises(ValueError)),
+        ],
+    )
+    def test_func_err(self, x, kwargs, with_condition):
+        """Test the `func_err` method."""
+        params = self._test_params
+        param_errors = self._test_param_errors
+        ff_obj = self.ff_class(params=params, param_errors=param_errors)
+
+        with with_condition:
+            results = ff_obj.func_err(x, **kwargs)
+            if "rety" in kwargs and kwargs["rety"]:
+                y_err, y = results
+            else:
+                y_err = results
+                y = None
+
+            x_err = None
+            if "x_err" in kwargs:
+                x_err = kwargs["x_err"]
+            if isinstance(x, list):
+                x = np.array(x)
+            y_err_expected = self.func_err(x, params, param_errors, x_err=x_err)
+
+            assert np.allclose(y_err, y_err_expected)
+
+            if y is not None:
+                assert np.allclose(y, self.func(x, *params))
+
+    @pytest.mark.parametrize(
+        "x, kwargs, with_condition",
+        [
+            (0, {}, does_not_raise()),
+            (0, {"reterr": True}, does_not_raise()),
+            (1.0, {}, does_not_raise()),
+            (1.0, {"reterr": True}, does_not_raise()),
+            (np.linspace(10, 30, num=20), {}, does_not_raise()),
+            (np.linspace(10, 30, num=20), {"reterr": True}, does_not_raise()),
+            ([4, 5, 6], {}, does_not_raise()),
+            ([4, 5, 6], {"x_err": 0.05, "reterr": True}, does_not_raise()),
+            ("hello", {}, pytest.raises(TypeError)),
+            (5, {"x_err": [1, 2], "reterr": True}, pytest.raises(ValueError)),
+        ],
+    )
+    def test_call(self, x, kwargs, with_condition):
+        """Test __call__ behavior."""
+        params = self._test_params
+        param_errors = self._test_param_errors
+        ff_obj = self.ff_class(params=params, param_errors=param_errors)
+
+        reterr = False
+        if "reterr" in kwargs:
+            reterr = kwargs["reterr"]
+        x_err = None
+        if "x_err" in kwargs:
+            x_err = kwargs["x_err"]
+
+        with with_condition:
+            results = ff_obj(x, **kwargs)
+            if reterr:
+                y = results[0]
+                y_err = results[1]
+            else:
+                y = results
+
+            if isinstance(x, list):
+                x = np.array(x)
+            y_expected = self.func(x, *params)
+
+            assert np.allclose(y, y_expected)
+
+            if reterr:
+                y_err_expected = self.func_err(x, params, param_errors, x_err=x_err)
+                assert np.allclose(y_err, y_err_expected)
+
+    @abstractmethod
+    def test_root_solve(self):
+        ...
+
+    def test_curve_fit(self):
+        """Test the `curve_fit` method."""
+        ff_obj = self.ff_class()
+
+        xdata = np.linspace(-10, 10)
+        ydata = self.func(xdata, *self._test_params)
+
+        assert ff_obj.params is None
+        assert ff_obj.param_errors is None
+        assert ff_obj.rsq is None
+        assert ff_obj.curve_fit_results is None
+
+        ff_obj.curve_fit(xdata, ydata)
+
+        assert ff_obj.curve_fit_results is not None
+        assert np.isclose(ff_obj.rsq, 1.0)
+        assert np.allclose(
+            ff_obj.param_errors, tuple([0] * len(ff_obj.param_names)), atol=1.5e-8,
+        )
+        assert np.allclose(ff_obj.params, self._test_params)
+
+
+class TestFFExponential(BaseFFTests):
+    """
+    Tests for fit function class `plasmapy.analysis.fit_functions.Exponential`.
+    """
+
+    ff_class = ffuncs.Exponential
+    _test_params = (5.0, 1.0)
+    _test_param_errors = (0.1, 0.1)
+    _test_param_names = ("a", "alpha")
+    _test_latex_str = fr"a \, \exp(\alpha x)"
+    _test__str__ = f"f(x) = a exp(alpha x)"
+
+    @staticmethod
+    def func(x, a, alpha):
+        return a * np.exp(alpha * x)
+
+    def func_err(self, x, params, param_errors, x_err=None):
+        a, alpha = params
+        a_err, alpha_err = param_errors
+        y = self.func(x, *params)
+
+        a_term = (a_err / a) ** 2
+        alpha_term = (x * alpha_err) ** 2
+
+        err = a_term + alpha_term
+
+        if x_err is not None:
+            x_term = (alpha * x_err) ** 2
+            err += x_term
+
+        err = np.abs(y) * np.sqrt(err)
+
+        return err
+
+    def test_root_solve(self):
+        ff_obj = self.ff_class(params=(1, 1), param_errors=(0, 0))
+        root, err = ff_obj.root_solve()
+        assert np.isnan(root)
+        assert np.isnan(err)
+
+
+class TestFFExponentialPlusLinear(BaseFFTests):
+    """
+    Tests for fit function class
+    `plasmapy.analysis.fit_functions.ExponentialPlusLinear`.
+    """
+
+    ff_class = ffuncs.ExponentialPlusLinear
+    _test_params = (2.0, 1.0, 5.0, -10.0)
+    _test_param_errors = (0.1, 0.1, 0.1, 0.1)
+    _test_param_names = ("a", "alpha", "m", "b")
+    _test_latex_str = fr"a \, \exp(\alpha x) + m x + b"
+    _test__str__ = f"f(x) = a exp(alpha x) + m x + b"
+
+    @staticmethod
+    def func(x, a, alpha, m, b):
+        return a * np.exp(alpha * x) + m * x + b
+
+    def func_err(self, x, params, param_errors, x_err=None):
+        a, alpha, m, b = params
+        a_err, alpha_err, m_err, b_err = param_errors
+
+        exp_y = a * np.exp(alpha * x)
+
+        a_term = (exp_y * a_err / a) ** 2
+        alpha_term = (exp_y * x * alpha_err) ** 2
+        m_term = (m_err * x) ** 2
+        b_term = b_err ** 2
+
+        err = a_term + alpha_term + m_term + b_term
+
+        if x_err is not None:
+            x_term = (exp_y * alpha * x_err) ** 2
+            x_term += (m * x_err) ** 2
+            x_term += 2 * a * alpha * m * np.exp(alpha * x) * (x_err ** 2)
+
+            err += x_term
+
+        err = np.sqrt(err)
+
+        return err
+
+    def test_root_solve(self):
+        ff_obj = self.ff_class(params=(5.0, 0.5, 1.0, 5.0), param_errors=(0, 0, 0, 0))
+        root, err = ff_obj.root_solve(-5)
+        assert np.isclose(root, -5.345338)
+        assert np.isnan(err)
+
+
+class TestFFExponentialPlusOffset(BaseFFTests):
+    """
+    Tests for fit function class
+    `plasmapy.analysis.fit_functions.ExponentialPlusOffset`.
+    """
+
+    ff_class = ffuncs.ExponentialPlusOffset
+    _test_params = (2.0, 1.0, -10.0)
+    _test_param_errors = (0.1, 0.1, 0.1)
+    _test_param_names = ("a", "alpha", "b")
+    _test_latex_str = fr"a \, \exp(\alpha x) + b"
+    _test__str__ = f"f(x) = a exp(alpha x) + b"
+
+    @staticmethod
+    def func(x, a, alpha, b):
+        return a * np.exp(alpha * x) + b
+
+    def func_err(self, x, params, param_errors, x_err=None):
+        a, alpha, b = params
+        a_err, alpha_err, b_err = param_errors
+
+        exp_y = a * np.exp(alpha * x)
+
+        a_term = (exp_y * a_err / a) ** 2
+        alpha_term = (exp_y * x * alpha_err) ** 2
+        b_term = b_err ** 2
+
+        err = a_term + alpha_term + b_term
+
+        if x_err is not None:
+            x_term = (exp_y * alpha * x_err) ** 2
+            err += x_term
+
+        err = np.sqrt(err)
+
+        return err
+
+    def test_root_solve(self):
+        ff_obj = self.ff_class(params=(3.0, 0.5, -5.0), param_errors=(0, 0, 0))
+        root, err = ff_obj.root_solve()
+        assert root == np.log(5.0 / 3.0) / 0.5
+        assert err == 0
+
+        ff_obj.params = (3.0, 0.5, 5.0)
+        root, err = ff_obj.root_solve()
+        assert np.isnan(root)
+        assert np.isnan(err)
+
+
+class TestFFLinear(BaseFFTests):
+    """
+    Tests for fit function class `plasmapy.analysis.fit_functions.Linear`.
+    """
+
+    ff_class = ffuncs.Linear
+    _test_params = (5.0, 4.0)
+    _test_param_errors = (0.1, 0.1)
+    _test_param_names = ("m", "b")
+    _test_latex_str = fr"m x + b"
+    _test__str__ = f"f(x) = m x + b"
+
+    @staticmethod
+    def func(x, m, b):
+        return m * x + b
+
+    def func_err(self, x, params, param_errors, x_err=None):
+        m, b = params
+        m_err, b_err = param_errors
+
+        m_term = (m_err * x) ** 2
+        b_term = b_err ** 2
+        err = m_term + b_term
+
+        if x_err is not None:
+            x_term = (m * x_err) ** 2
+            err += x_term
+        err = np.sqrt(err)
+
+        return err
+
+    def test_root_solve(self):
+        ff_obj = self.ff_class(params=(1, 1), param_errors=(0, 0))
+        assert ff_obj.root_solve() == (-1, 0)
+
+        ff_obj.params = (5.0, 1.3)
+        ff_obj.param_errors = (0.1, 0.1)
+        root, err = ff_obj.root_solve()
+        assert root == -1.3 / 5.0
+        assert err == np.abs(root) * np.sqrt((0.1 / 5.0) ** 2 + (0.1 / 1.3) ** 2)