Repository: PlasmaPy/PlasmaPy
Pull request #932 (2020-10-13 19:31:02)

Create plasmapy.dispersion.dispersionsolver & associated functionality
<!--
Thank you for contributing to PlasmaPy! Here's a bunch of pointers to
make things easier for all of us:

* If this PR will solve an issue tracked by GitHub, then please add
  "Closes #42" so the issue automatically closes once this pull request
  is merged.  In this example, issue #42 would have been closed.
  If your PR will not completely solve the issue, then please still reference
  still reference the issue like "issue #42".  (For more info see [GitHub's primer
  on linking issues](https://help.github.com/en/github/managing-your-work-on-github/linking-a-pull-request-to-an-issue#linking-a-pull-request-to-an-issue-using-a-keyword).)

* Remember to add some description of your changes in this text box.

* If your pull request is not yet ready for review - either because
  you're just looking for feedback on a change or because you're not
  perfectly satisfied with your change - submit it as a draft pull request.
  Remember to change its' status once it's ready.

* If this is your first contribution, please add your name to the author
  list in `docs/about/credits.rst`.

* Feel free to chat with other developers on our Matrix channel at:
   https://riot.im/app/#/room/#plasmapy:openastronomy.org

* We have a developer's guide to help answer some of your questions.
  http://docs.plasmapy.org/en/latest/development/index.html

Many thanks in advance for following these pointers and for being willing to contribute!

When submitting a pull request, please ensure that you can (eventually,
sometime before it is merged) check the following basic requirements:

-->

- [x] I have added a changelog entry for this pull request.

<!--

In short: A changelog entry is a short description of your PR's changes.
Each entry is written in a `<PULL REQUEST>.<TYPE>.rst` file and stored in
the `changelog` directory,  where `<PULL REQUEST>` is a pull request
number and `<TYPE>` is one of:

* `breaking`: A change which requires users to change code and is not backwards compatible. (Not to be used for removal of deprecated features.)
* `feature`: New user facing features and any new behavior.
* `bugfix`: Fixes a reported bug.
* `doc`: Documentation addition or improvement, like rewording an entire session or adding missing docs.
* `removal`: Feature deprecation and/or feature removal.
* `trivial`: A change which has no user facing effect or is tiny change.

A PR number is generated after you successfully submit a new PR, so the changelog
file can only be added after you open the PR.```

For more information, see:
https://github.com/PlasmaPy/PlasmaPy/blob/master/changelog/README.rst

--->

- [x] If adding new functionality, I have added tests and
      docstrings.

<!--
(Tests pop up at the bottom, in the checks box.)
-->

- [x] I have fixed any newly failing tests.

<!--
(If you're unsure why they're failing, ask!)
-->


review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/PlasmaPy/PlasmaPy/pull/932"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932?src=pr&el=h1) Report
> Merging [#932](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/6f277181fccb0c26cc623942892f7b18daae4a66?el=desc) will **increase** coverage by `0.02%`.
> The diff coverage is `100.00%`.

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

```diff
@@            Coverage Diff             @@
##           master     #932      +/-   ##
==========================================
+ Coverage   96.19%   96.21%   +0.02%     
==========================================
  Files          60       61       +1     
  Lines        5413     5443      +30     
==========================================
+ Hits         5207     5237      +30     
  Misses        206      206              
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/ionization.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L2lvbml6YXRpb24ucHk=) | `100.00% <ø> (ø)` | |
| [plasmapy/dispersion/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932/diff?src=pr&el=tree#diff-cGxhc21hcHkvZGlzcGVyc2lvbi9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [plasmapy/dispersion/two\_fluid\_dispersion\_solver.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932/diff?src=pr&el=tree#diff-cGxhc21hcHkvZGlzcGVyc2lvbi90d29fZmx1aWRfZGlzcGVyc2lvbl9zb2x2ZXIucHk=) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932?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/932?src=pr&el=footer). Last update [6f27718...e6aedf2](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/932?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

rocco8773: @ahmadryan Ok, I think I've finished updating the tests and the jupyter notebook.  For the notebook I added more background of the dispersion relation, updated the plots so they work better with RTD, and added a section that duplicates Figure 1 of [Bellan 2012](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856).  If you don't mind giving it a read through to ensure everything makes sense and reads correctly, then that would be much appreciated.  No matter how many time I reread the things I write I always seem to miss wordage that make me go 🤦🏻‍♂️.
rocco8773: @namurphy Thanks for the quick read through!  I'm going to leave the PR open for now to give @ahmadryan a chance to read through and comment.  I'm hoping we'll merge Monday or Tuesday.

Patch:
diff --git a/.gitignore b/.gitignore
index 4ce6c4d8..51f1fef2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -102,6 +102,9 @@ ENV/
 # Rope project settings
 .ropeproject
 
+# VS code project settings
+.vscode
+
 # mkdocs documentation
 /site
 
diff --git a/changelog/932.feature.rst b/changelog/932.feature.rst
new file mode 100644
index 00000000..10513c34
--- /dev/null
+++ b/changelog/932.feature.rst
@@ -0,0 +1,5 @@
+Add the function
+:func:`~plasmapy.dispersion.two_fluid_dispersion.two_fluid_dispersion_solution` to
+`plasmapy.dispersion`, which gives an analytical solution to the dispersion relation as
+derived by P. M. Bellan 2012 (DOI: `10.1029/2012JA017856
+<https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856>`_).
diff --git a/changelog/962.breaking.rst b/changelog/962.breaking.rst
new file mode 100644
index 00000000..58a0f2da
--- /dev/null
+++ b/changelog/962.breaking.rst
@@ -0,0 +1,1 @@
+Changed the names of the supported methods of computing the Coulomb logarithm to improve readability. Expanded the documentation of these methods.
diff --git a/docs/api_static/plasmapy.dispersion.dispersionfunction.rst b/docs/api_static/plasmapy.dispersion.dispersionfunction.rst
index 2a884ec4..7bb428ca 100644
--- a/docs/api_static/plasmapy.dispersion.dispersionfunction.rst
+++ b/docs/api_static/plasmapy.dispersion.dispersionfunction.rst
@@ -15,4 +15,4 @@ Example Notebooks
 
 .. nbgallery::
 
-    /notebooks/dispersion_function
+    /notebooks/dispersion/dispersion_function
diff --git a/docs/api_static/plasmapy.dispersion.two_fluid_dispersion.rst b/docs/api_static/plasmapy.dispersion.two_fluid_dispersion.rst
new file mode 100644
index 00000000..127f56d1
--- /dev/null
+++ b/docs/api_static/plasmapy.dispersion.two_fluid_dispersion.rst
@@ -0,0 +1,18 @@
+:orphan:
+
+`plasmapy.dispersion.two_fluid_dispersion`
+==========================================
+
+.. currentmodule:: plasmapy.dispersion.two_fluid_dispersion
+
+.. automodapi:: plasmapy.dispersion.two_fluid_dispersion
+   :include-all-objects:
+   :no-heading:
+   :headings: "-^"
+
+Example Notebooks
+-----------------
+
+.. nbgallery::
+
+    /notebooks/dispersion/two_fluid_dispersion
diff --git a/docs/api_static/plasmapy.rst b/docs/api_static/plasmapy.rst
index 433716ba..4f728a15 100644
--- a/docs/api_static/plasmapy.rst
+++ b/docs/api_static/plasmapy.rst
@@ -13,6 +13,7 @@ Sub-Packages & Modules
 .. autosummary::
 
    diagnostics
+   dispersion
    formulary
    particles
    plasma
diff --git a/docs/dispersion/index.rst b/docs/dispersion/index.rst
index 0fc80bbc..f3eb314b 100644
--- a/docs/dispersion/index.rst
+++ b/docs/dispersion/index.rst
@@ -2,14 +2,18 @@
 Dispersion (`plasmapy.dispersion`)
 **********************************
 
+.. py:currentmodule:: plasmapy.dispersion
+
 .. automodule:: plasmapy.dispersion
 
 Example Notebooks
 -----------------
 
 .. nbgallery::
+   :glob:
+   :reversed:
 
-   /notebooks/dispersion_function
+   /notebooks/dispersion/*
 
 Sub-Packages & Modules
 ----------------------
@@ -17,6 +21,7 @@ Sub-Packages & Modules
 .. autosummary::
 
    dispersionfunction
+   two_fluid_dispersion
 
 .. automodapi:: plasmapy.dispersion
    :no-main-docstr:
diff --git a/docs/index.rst b/docs/index.rst
index 791d6af5..a2638ae7 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -10,7 +10,7 @@ PlasmaPy Documentation
    :hidden:
 
    notebooks/physics
-   notebooks/dispersion_function
+   notebooks/dispersion/two_fluid_dispersion
    notebooks/thomson
 
 .. image:: _static/graphic-circular.png
diff --git a/docs/notebooks/dispersion_function.ipynb b/docs/notebooks/dispersion/dispersion_function.ipynb
similarity index 99%
rename from docs/notebooks/dispersion_function.ipynb
rename to docs/notebooks/dispersion/dispersion_function.ipynb
index 23af7060..f1a0369b 100644
--- a/docs/notebooks/dispersion_function.ipynb
+++ b/docs/notebooks/dispersion/dispersion_function.ipynb
@@ -66,12 +66,12 @@
    ]
   },
   {
-   "cell_type": "markdown",
-   "metadata": {},
    "source": [
     "We'll now make some sample data to visualize the dispersion function:\n",
     "\n"
-   ]
+   ],
+   "cell_type": "markdown",
+   "metadata": {}
   },
   {
    "cell_type": "code",
@@ -296,7 +296,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.7.4-final"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {
diff --git a/docs/notebooks/dispersion/two_fluid_dispersion.ipynb b/docs/notebooks/dispersion/two_fluid_dispersion.ipynb
new file mode 100644
index 00000000..9358feee
--- /dev/null
+++ b/docs/notebooks/dispersion/two_fluid_dispersion.ipynb
@@ -0,0 +1,531 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Dispersion: A Full Two Fluid Solution\n",
+    "\n",
+    "[tfds]: ../../api/plasmapy.dispersion.two_fluid_dispersion.two_fluid_dispersion_solution.rst\n",
+    "[bellan2012]: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856\n",
+    "[stringer1963]: https://doi.org/10.1088/0368-3281/5/2/304\n",
+    "\n",
+    "This notebook walks through the functionalty of the [two_fluid_dispersion_solution()][tfds] function.  This function computes the wave frequencies for given wavenumbers and plasma parameters based on the analytical solution presented by [Bellan 2012][bellan2012] to the [Stringer 1963][stringer1963] two fluid dispersion relation.  The two fluid dispersion equaiton assumes a uniform magnetic field, a zero D.C. electric field, and low-frequency waves $\\omega / k c \\ll 1$ which equates to\n",
+    "\n",
+    "$$\n",
+    "    \\left( \\cos^2 \\theta - Q \\frac{\\omega^2}{k^2 {v_A}^2} \\right)\n",
+    "    \\left[\n",
+    "        \\left( \\cos^2 \\theta - \\frac{\\omega^2}{k^2 {c_s}^2} \\right)\n",
+    "        - Q \\frac{\\omega^2}{k^2 {v_A}^2} \\left(\n",
+    "            1 - \\frac{\\omega^2}{k^2 {c_s}^2}\n",
+    "        \\right)\n",
+    "    \\right] \\\\\n",
+    "    = \\left(1 - \\frac{\\omega^2}{k^2 {c_s}^2} \\right)\n",
+    "      \\frac{\\omega^2}{{\\omega_{ci}}^2} \\cos^2 \\theta\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$$Q = 1 + k^2 c^2/{\\omega_{pe}}^2$$\n",
+    "\n",
+    "$$\\cos \\theta = \\frac{k_z}{k}$$\n",
+    "\n",
+    "$$\\mathbf{B_o} = B_{o} \\mathbf{\\hat{z}}$$\n",
+    "\n",
+    "$\\omega$ is the wave frequency, $k$ is the wavenumber, $v_A$ is the Alfvén velocity, $c_s$ is the sound speed, $\\omega_{ci}$ is the ion gyrofrequency, and $\\omega_{pe}$ is the electron plasma frequency.\n",
+    "\n",
+    "The approach outlined in Section 5 of [Bellan 2012][bellan2012] produces exact roots to the above dispersion equation for all three modes (fast, acoustic, and Alfvén) without having to make additional approximations.  The following dispersion relation is what the [two_fluid_dispersion_solution()][tfds] function computes.\n",
+    "\n",
+    "$$\n",
+    "    \\frac{\\omega}{\\omega_{ci}} = \\sqrt{\n",
+    "        2 \\Lambda \\sqrt{-\\frac{P}{3}} \\cos\\left(\n",
+    "            \\frac{1}{3} \\cos^{-1}\\left(\n",
+    "                \\frac{3q}{2p} \\sqrt{-\\frac{3}{p}}\n",
+    "            \\right)\n",
+    "            - \\frac{2 \\pi}{3}j\n",
+    "        \\right)\n",
+    "        + \\frac{\\Lambda A}{3}\n",
+    "    }\n",
+    "$$\n",
+    "\n",
+    "where $j = 0$ represents the fast mode, $j = 1$ represents the Alfvén mode, and $j = 2$ represents the acoustic mode.  Additionally,\n",
+    "\n",
+    "$$p = \\frac{3B-A^2}{3} \\; , \\; q = \\frac{9AB-2A^3-27C}{27}$$\n",
+    "\n",
+    "$$A = \\frac{Q + Q^2 \\beta + Q \\alpha + \\alpha \\Lambda}{Q^2} \\;\n",
+    "    , \\; B = \\alpha \\frac{1 + 2 Q \\beta + \\Lambda \\beta}{Q^2} \\;\n",
+    "    , \\; C = \\frac{\\alpha^2 \\beta}{Q^2}$$\n",
+    "\n",
+    "$$\\alpha = \\cos^2 \\theta \\;\n",
+    "    , \\; \\beta = \\left( \\frac{c_s}{v_A}\\right)^2 \\;\n",
+    "    , \\; \\Lambda = \\left( \\frac{k v_{A}}{\\omega_{ci}}\\right)^2$$\n",
+    "\n",
+    "## Contents:\n",
+    "\n",
+    "1. [Wave Propagating at 45 Degrees](#Wave-Propagating-at-45-Degrees)\n",
+    "2. [Wave frequencies on the k-theta plane](#Wave-frequencies-on-the-k-theta-plane)\n",
+    "3. [Reproduce Figure 1 from Bellan 2012](#Reproduce-Figure-1-from-Bellan-2012)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import astropy.units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from astropy.constants.si import c\n",
+    "from matplotlib import colors\n",
+    "from matplotlib.ticker import MultipleLocator\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "\n",
+    "from plasmapy.dispersion.two_fluid_dispersion import two_fluid_dispersion_solution\n",
+    "from plasmapy.formulary import parameters as pfp\n",
+    "from plasmapy.particles import Particle\n",
+    "\n",
+    "plt.rcParams['figure.figsize'] = [10.5, 0.56 * 10.5]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Wave Propagating at 45 Degrees\n",
+    "\n",
+    "Below we define the required parameters to compute the wave frequencies."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define input parameters\n",
+    "inputs = {\n",
+    "    \"k\": np.linspace(10**-7, 10**-2, 10000) * u.rad / u.m,\n",
+    "    \"theta\": 45 * u.deg,\n",
+    "    \"n_i\": 5 * u.cm ** -3,\n",
+    "    \"B\": 8.3E-9 * u.T,\n",
+    "    \"T_e\": 1.6e6 * u.K,\n",
+    "    \"T_i\": 4.e5 * u.K,\n",
+    "    \"ion\": Particle(\"p+\"),\n",
+    "}\n",
+    "\n",
+    "# a few useful plasma parameters\n",
+    "params = {\n",
+    "    \"n_e\": inputs[\"n_i\"] * abs(inputs[\"ion\"].integer_charge),\n",
+    "    \"cs\": pfp.ion_sound_speed(\n",
+    "        inputs[\"T_e\"],\n",
+    "        inputs[\"T_i\"],\n",
+    "        inputs[\"ion\"],\n",
+    "    ),\n",
+    "    \"va\": pfp.Alfven_speed(\n",
+    "        inputs[\"B\"],\n",
+    "        inputs[\"n_i\"],\n",
+    "        ion=inputs[\"ion\"],\n",
+    "    ),\n",
+    "    \"wci\": pfp.gyrofrequency(inputs[\"B\"], inputs[\"ion\"]),\n",
+    "}\n",
+    "params[\"lpe\"] = pfp.inertial_length(params[\"n_e\"], \"e-\")\n",
+    "params[\"wpe\"] = pfp.plasma_frequency(params[\"n_e\"], \"e-\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The computed wave frequencies ($rad/s$) are returned in a dictionary with keys representing the wave modes and the values being an Astropy [Quantity](https://docs.astropy.org/en/stable/api/astropy.units.Quantity.html#astropy.units.Quantity).  Since our inputs had a scalar $\\theta$ and a 1D array of $k$'s, the computed wave frequencies will be a 1D array of size equal to the size of the $k$ array."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# compute\n",
+    "omegas = two_fluid_dispersion_solution(**inputs)\n",
+    "\n",
+    "(list(omegas.keys()), omegas[\"fast_mode\"], omegas[\"fast_mode\"].shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's plot the results of each wave mode."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs = 14  # default font size\n",
+    "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
+    "figheight = 1.6 * figheight\n",
+    "fig = plt.figure(figsize=[figwidth, figheight])\n",
+    "\n",
+    "# normalize data\n",
+    "k_prime = inputs[\"k\"] * params[\"lpe\"]\n",
+    "\n",
+    "# plot\n",
+    "plt.plot(\n",
+    "    k_prime,\n",
+    "    np.real(omegas[\"fast_mode\"] / params[\"wpe\"]),\n",
+    "    'r.',\n",
+    "    ms=1,\n",
+    "    label='Fast',\n",
+    ")\n",
+    "ax = plt.gca()\n",
+    "ax.plot(\n",
+    "    k_prime,\n",
+    "    np.real(omegas[\"alfven_mode\"] / params[\"wpe\"]),\n",
+    "    'b.',\n",
+    "    ms=1,\n",
+    "    label='Alfvèn',\n",
+    ")\n",
+    "ax.plot(\n",
+    "    k_prime,\n",
+    "    np.real(omegas[\"acoustic_mode\"] / params[\"wpe\"]),\n",
+    "    'g.',\n",
+    "    ms=1,\n",
+    "    label='Acoustic',\n",
+    ")\n",
+    "\n",
+    "# adjust axes\n",
+    "ax.set_xlabel(r\"$kc / \\omega_{pe}$\", fontsize=fs)\n",
+    "ax.set_ylabel(r\"$Re(\\omega / \\omega_{pe})$\", fontsize=fs)\n",
+    "ax.set_yscale('log')\n",
+    "ax.set_xscale('log')\n",
+    "ax.set_ylim(1e-6, 2e-2)\n",
+    "ax.tick_params(\n",
+    "    which=\"both\",\n",
+    "    direction=\"in\",\n",
+    "    width=1,\n",
+    "    labelsize=fs,\n",
+    "    right=True,\n",
+    "    length=5,\n",
+    ")\n",
+    "\n",
+    "# annotate\n",
+    "text = (\n",
+    "    f\"$v_A/c_s = {params['va'] / params['cs']:.1f} \\qquad \"\n",
+    "    f\"c/v_A = 10^{np.log10(c / params['va']):.0f} \\qquad \"\n",
+    "    f\"\\\\theta = {inputs['theta'].value:.0f}\"\n",
+    "    \"^{\\\\circ}$\"\n",
+    ")\n",
+    "ax.text(0.25, 0.95, text, transform=ax.transAxes, fontsize=18)\n",
+    "ax.legend(loc=\"upper left\", markerscale=5, fontsize=fs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Wave frequencies on the k-theta plane\n",
+    "\n",
+    "Let us now look at the distribution of $\\omega$ on a $k$-$\\theta$ plane."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define input parameters\n",
+    "inputs = {\n",
+    "    \"k\": np.linspace(10**-7, 10**-2, 10000) * u.rad / u.m,\n",
+    "    \"theta\": np.linspace(5, 85, 100) * u.deg,\n",
+    "    \"n_i\": 5 * u.cm ** -3,\n",
+    "    \"B\": 8.3E-9 * u.T,\n",
+    "    \"T_e\": 1.6e6 * u.K,\n",
+    "    \"T_i\": 4.e5 * u.K,\n",
+    "    \"ion\": Particle(\"p+\"),\n",
+    "}\n",
+    "\n",
+    "# a few useful plasma parameters\n",
+    "params = {\n",
+    "    \"n_e\": inputs[\"n_i\"] * abs(inputs[\"ion\"].integer_charge),\n",
+    "    \"cs\": pfp.ion_sound_speed(\n",
+    "        inputs[\"T_e\"],\n",
+    "        inputs[\"T_i\"],\n",
+    "        inputs[\"ion\"],\n",
+    "    ),\n",
+    "    \"va\": pfp.Alfven_speed(\n",
+    "        inputs[\"B\"],\n",
+    "        inputs[\"n_i\"],\n",
+    "        ion=inputs[\"ion\"],\n",
+    "    ),\n",
+    "    \"wci\": pfp.gyrofrequency(inputs[\"B\"], inputs[\"ion\"]),\n",
+    "}\n",
+    "params[\"lpe\"] = pfp.inertial_length(params[\"n_e\"], \"e-\")\n",
+    "params[\"wpe\"] = pfp.plasma_frequency(params[\"n_e\"], \"e-\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since the $\\theta$ and $k$ values are now 1-D arrays, the returned wave frequencies will be 2-D arrays with the first dimension matching the size of $k$ and the second dimension matching the size of $\\theta$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# compute\n",
+    "omegas = two_fluid_dispersion_solution(**inputs)\n",
+    "\n",
+    "(\n",
+    "    omegas[\"fast_mode\"].shape,\n",
+    "    omegas[\"fast_mode\"].shape[0] == inputs[\"k\"].size,\n",
+    "    omegas[\"fast_mode\"].shape[1] == inputs[\"theta\"].size,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's plot (the fast mode)!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fs = 14  # default font size\n",
+    "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
+    "figheight = 1.6 * figheight\n",
+    "fig = plt.figure(figsize=[figwidth, figheight])\n",
+    "\n",
+    "# normalize data\n",
+    "k_prime = inputs[\"k\"] * params[\"lpe\"]\n",
+    "zdata = (\n",
+    "    np.transpose(np.real(omegas[\"fast_mode\"].value))\n",
+    "    / params[\"wpe\"].value\n",
+    ")\n",
+    "\n",
+    "# plot\n",
+    "im = plt.imshow(\n",
+    "    zdata,\n",
+    "    aspect='auto',\n",
+    "    origin='lower',\n",
+    "    extent=[\n",
+    "        np.min(k_prime.value),\n",
+    "        np.max(k_prime.value),\n",
+    "        np.min(inputs[\"theta\"].value),\n",
+    "        np.max(inputs[\"theta\"].value),\n",
+    "    ],\n",
+    "    interpolation=None,\n",
+    "    cmap=plt.cm.Spectral,\n",
+    ")\n",
+    "ax = plt.gca()\n",
+    "\n",
+    "# # adjust axes\n",
+    "ax.set_xscale('linear')\n",
+    "ax.set_xlabel(r'$kc/\\omega_{pe}$', fontsize=fs)\n",
+    "ax.set_ylabel(r'$\\theta$ [$deg.$]', fontsize=fs)\n",
+    "ax.tick_params(\n",
+    "    which=\"both\",\n",
+    "    direction=\"in\",\n",
+    "    width=2,\n",
+    "    labelsize=fs,\n",
+    "    right=True,\n",
+    "    top=True,\n",
+    "    length=10,\n",
+    ")\n",
+    "\n",
+    "# Add colorbar\n",
+    "divider = make_axes_locatable(ax)\n",
+    "cax = divider.append_axes(\"top\", size=\"5%\", pad=0.07)\n",
+    "cbar = plt.colorbar(\n",
+    "    im,\n",
+    "    cax=cax,\n",
+    "    orientation='horizontal',\n",
+    "    ticks=None,\n",
+    "    fraction=0.05,\n",
+    "    pad=0.0,\n",
+    ")\n",
+    "cbar.ax.tick_params(\n",
+    "    axis='x',\n",
+    "    direction='in',\n",
+    "    width=2,\n",
+    "    length=10,\n",
+    "    top=True,\n",
+    "    bottom=False,\n",
+    "    labelsize=fs,    \n",
+    "    pad=0.0,\n",
+    "    labeltop=True,\n",
+    "    labelbottom=False,\n",
+    ")\n",
+    "cbar.ax.xaxis.set_label_position('top')\n",
+    "cbar.set_label(r'$\\omega/\\omega_{pe}$', fontsize=fs, labelpad=8)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Reproduce Figure 1 from Bellan 2012\n",
+    "\n",
+    "[bellan2012]: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856\n",
+    "\n",
+    "\n",
+    "Figure 1 of [Bellan 2012][bellan2012] chooses parameters such that $\\beta = 0.4$ and $\\Lambda=0.4$.  Below we define parameters to approximate Bellan's assumptions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define input parameters\n",
+    "inputs = {\n",
+    "    \"B\": 400e-4 * u.T,\n",
+    "    \"ion\": Particle(\"He+\"),\n",
+    "    \"n_i\": 6.358e19 * u.m ** -3,\n",
+    "    \"T_e\": 20 * u.eV,\n",
+    "    \"T_i\": 10 * u.eV,\n",
+    "    \"theta\": np.linspace(0, 90) * u.deg,\n",
+    "    \"k\": (2 * np.pi * u.rad) / (0.56547 * u.m),\n",
+    "}\n",
+    "\n",
+    "# a few useful plasma parameters\n",
+    "params = {\n",
+    "    \"n_e\": inputs[\"n_i\"] * abs(inputs[\"ion\"].integer_charge),\n",
+    "    \"cs\": pfp.cs_(inputs[\"T_e\"], inputs[\"T_i\"], inputs[\"ion\"]),\n",
+    "    \"wci\": pfp.wc_(inputs[\"B\"], inputs[\"ion\"]),\n",
+    "    \"va\": pfp.va_(inputs[\"B\"], inputs[\"n_i\"], ion=inputs[\"ion\"]),\n",
+    "}\n",
+    "params[\"beta\"] = (params[\"cs\"] / params[\"va\"]).value ** 2\n",
+    "params[\"wpe\"] = pfp.wp_(params[\"n_e\"], \"e-\")\n",
+    "params[\"Lambda\"] = (inputs[\"k\"] * params[\"va\"] / params[\"wci\"]).value ** 2\n",
+    "\n",
+    "(params[\"beta\"], params[\"Lambda\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# compute\n",
+    "omegas = two_fluid_dispersion_solution(**inputs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# generate data for plots\n",
+    "plt_vals = {}\n",
+    "for mode, arr in omegas.items():\n",
+    "    norm = (np.absolute(arr) / (inputs[\"k\"] * params[\"va\"])).value ** 2\n",
+    "    plt_vals[mode] = {\n",
+    "        \"x\": norm * np.sin(inputs[\"theta\"].to(u.rad).value),\n",
+    "        \"y\": norm * np.cos(inputs[\"theta\"].to(u.rad).value),\n",
+    "    }"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+     "tooltip": "Full Two Fluid Dispersion Solution"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fs = 14  # default font size\n",
+    "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
+    "figheight = 1.6 * figheight\n",
+    "fig = plt.figure(figsize=[figwidth, figheight])\n",
+    "\n",
+    "# Fast mode\n",
+    "plt.plot(\n",
+    "    plt_vals[\"fast_mode\"][\"x\"],\n",
+    "    plt_vals[\"fast_mode\"][\"y\"],\n",
+    "    linewidth=2,\n",
+    "    label=\"Fast\",\n",
+    ")\n",
+    "ax = plt.gca()\n",
+    "\n",
+    "# adjust axes\n",
+    "ax.set_xlabel(r\"$(\\omega / k v_A)^2 \\, \\sin \\theta$\", fontsize=fs)\n",
+    "ax.set_ylabel(r\"$(\\omega / k v_A)^2 \\, \\cos \\theta$\", fontsize=fs)\n",
+    "ax.set_xlim(0.0, 1.5)\n",
+    "ax.set_ylim(0.0, 2.0)\n",
+    "for spine in ax.spines.values():\n",
+    "    spine.set_linewidth(2)\n",
+    "ax.minorticks_on()\n",
+    "ax.tick_params(which=\"both\", labelsize=fs, width=2)\n",
+    "ax.tick_params(which=\"major\", length=10)\n",
+    "ax.tick_params(which=\"minor\", length=5)\n",
+    "ax.xaxis.set_major_locator(MultipleLocator(0.5))\n",
+    "ax.xaxis.set_minor_locator(MultipleLocator(0.1))\n",
+    "ax.yaxis.set_major_locator(MultipleLocator(0.5))\n",
+    "ax.yaxis.set_minor_locator(MultipleLocator(0.1))\n",
+    "\n",
+    "\n",
+    "# Alfven mode\n",
+    "plt.plot(\n",
+    "    plt_vals[\"alfven_mode\"][\"x\"],\n",
+    "    plt_vals[\"alfven_mode\"][\"y\"],\n",
+    "    linewidth=2,\n",
+    "    label=\"Alfv$\\`{e}$n\",\n",
+    ")\n",
+    "\n",
+    "# Acoustic mode\n",
+    "plt.plot(\n",
+    "    plt_vals[\"acoustic_mode\"][\"x\"],\n",
+    "    plt_vals[\"acoustic_mode\"][\"y\"],\n",
+    "    linewidth=2,\n",
+    "    label=\"Acoustic\",\n",
+    ")\n",
+    "\n",
+    "# annotations\n",
+    "plt.legend(fontsize=fs, loc=\"upper right\")"
+   ]
+  }
+ ],
+ "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/dispersion/__init__.py b/plasmapy/dispersion/__init__.py
index 18686fc5..ca680645 100644
--- a/plasmapy/dispersion/__init__.py
+++ b/plasmapy/dispersion/__init__.py
@@ -5,9 +5,15 @@
 __all__ = [
     "plasma_dispersion_func",
     "plasma_dispersion_func_deriv",
+    "tfds_",
+    "two_fluid_dispersion_solution",
 ]
 
 from plasmapy.dispersion.dispersionfunction import (
     plasma_dispersion_func,
     plasma_dispersion_func_deriv,
 )
+from plasmapy.dispersion.two_fluid_dispersion import (
+    tfds_,
+    two_fluid_dispersion_solution,
+)
diff --git a/plasmapy/dispersion/two_fluid_dispersion.py b/plasmapy/dispersion/two_fluid_dispersion.py
new file mode 100644
index 00000000..c228f00a
--- /dev/null
+++ b/plasmapy/dispersion/two_fluid_dispersion.py
@@ -0,0 +1,352 @@
+__all__ = ["two_fluid_dispersion_solution", "tfds_"]
+
+import astropy.units as u
+import numpy as np
+import warnings
+
+from astropy.constants.si import c
+from typing import Union
+
+from plasmapy.formulary import parameters as pfp
+from plasmapy.particles import Particle
+from plasmapy.particles.exceptions import ChargeError
+from plasmapy.utils.decorators import validate_quantities
+from plasmapy.utils.exceptions import PhysicsWarning
+
+
+@validate_quantities(
+    B={"can_be_negative": False},
+    n_i={"can_be_negative": False},
+    T_e={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+    T_i={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+)
+def two_fluid_dispersion_solution(
+    *,
+    B: u.T,
+    ion: Union[str, Particle],
+    k: u.rad / u.m,
+    n_i: u.m ** -3,
+    T_e: u.K,
+    T_i: u.K,
+    theta: u.deg,
+    gamma_e: Union[float, int] = 1,
+    gamma_i: Union[float, int] = 3,
+    z_mean: Union[float, int] = None,
+):
+    r"""
+    Using the solution provided by Bellan 2012, calculate the analytical
+    solution to the two fluid, low-frequency (:math:`\omega/kc \ll 1`) dispersion
+    relation presented by Stringer 1963.  This dispersion relation also
+    assummes a uniform magnetic field :math:`\mathbf{B_o}`, no D.C. electric
+    field :math:`\mathbf{E_o}=0`, and quasi-neutrality.  For more information
+    see the **Notes** section below.
+
+    **Aliases:** `tfds_`
+
+    Parameters
+    ----------
+    B : `~astropy.units.Quantity`
+        The magnetic field magnitude in units convertible to :math:`T`.
+    ion : `str` or `~plasmapy.particles.particle_class.Particle`
+        Representation of the ion species (e.g., ``'p'`` for protons, ``'D+'``
+        for deuterium, ``'He-4 +1'`` for singly ionized helium-4, etc.). If no
+        charge state information is provided, then the ions are assumed to be
+        singly ionized.
+    k : `~astropy.units.Quantity`, single valued or 1-D array
+        Wavenumber in units convertible to :math:`rad / m`.  Either single
+        valued or 1-D array of length :math:`N`.
+    n_i : `~astropy.units.Quantity`
+        Ion number density in units convertible to :math:`m^{-3}`.
+    T_e : `~astropy.units.Quantity`
+        The electron temperature in units of :math:`K` or :math:`eV`.
+    T_i : `~astropy.units.Quantity`
+        The ion temperature in units of :math:`K` or :math:`eV`.
+    theta : `~astropy.units.Quantity`, single valued or 1-D array
+        The angle of propagation of the wave with respect to the magnetic field,
+        :math:`\cos^{-1}(k_z / k)`, in units must be convertible to :math:`deg`.
+        Either single valued or 1-D array of size :math:`M`.
+    gamma_e : `float` or `int`, optional
+        The adiabatic index for electrons, which defaults to 1.  This
+        value assumes that the electrons are able to equalize their
+        temperature rapidly enough that the electrons are effectively
+        isothermal.
+    gamma_i : `float` or `int`, optional
+        The adiabatic index for ions, which defaults to 3. This value
+        assumes that ion motion has only one degree of freedom, namely
+        along magnetic field lines.
+    z_mean : `float` or int, optional
+        The average ionization state (arithmetic mean) of the ``ion`` composing
+        the plasma.  Will override any charge state defined by argument ``ion``.
+
+    Returns
+    -------
+    omega : Dict[str, `~astropy.units.Quantity`]
+        A dictionary of computed wave frequencies in units :math:`rad/s`.  The
+        dictionary contains three keys: ``'fast_mode'`` for the fast mode,
+        ``'alfven_mode'`` for the Alfvén mode, and ``'acoustic_mode'`` for the
+        ion-acoustic mode.  The value for each key will be a :math:`N x M` array.
+
+    Raises
+    ------
+    TypeError
+        If applicable arguments are not instances of `~astropy.units.Quantity` or
+        cannot be converted into one.
+
+    TypeError
+        If ``ion`` is not of type or convertible to `~plasmapy.particles.Particle`.
+
+    TypeError
+        If ``gamma_e``, ``gamma_i``, or``z_mean`` are not of type `int` or `float`.
+
+    ~astropy.units.UnitTypeError
+        If applicable arguments do not have units convertible to the expected
+        units.
+
+    ValueError
+        If any of ``B``, ``k``, ``n_i``, ``T_e``, or ``T_i`` is negative.
+
+    ValueError
+        If ``k`` is negative or zero.
+
+    ValueError
+        If ``ion`` is not of category ion or element.
+
+    ValueError
+        If ``B``, ``n_i``, ``T_e``, or ``T_I`` are not single valued
+        `astropy.units.Quantity` (i.e. an array).
+
+    ValueError
+        If ``k`` or ``theta`` are not single valued or a 1-D array.
+
+    Warns
+    -----
+    : `~plasmapy.utils.exceptions.PhysicsWarning`
+        When the computed wave frequencies violate the low-frequency
+        (:math:`\omega/kc \ll 1`) assumption of the dispersion relation.
+
+    Notes
+    -----
+
+    The complete dispersion equation presented by Springer 1963 [2]_ (equation 1
+    of Bellan 2012 [1]_) is:
+
+    .. math::
+        \left( \cos^2 \theta - Q \frac{\omega^2}{k^2 {v_A}^2} \right) &
+        \left[
+            \left( \cos^2 \theta - \frac{\omega^2}{k^2 {c_s}^2} \right)
+            - Q \frac{\omega^2}{k^2 {v_A}^2} \left(
+                1 - \frac{\omega^2}{k^2 {c_s}^2}
+            \right)
+        \right] \\
+            &= \left(1 - \frac{\omega^2}{k^2 {c_s}^2} \right)
+              \frac{\omega^2}{{\omega_{ci}}^2} \cos^2 \theta
+
+    where
+
+    .. math::
+        Q &= 1 + k^2 c^2/{\omega_{pe}}^2 \\
+        \cos \theta &= \frac{k_z}{k} \\
+        \mathbf{B_o} &= B_{o} \mathbf{\hat{z}}
+
+    :math:`\omega` is the wave frequency, :math:`k` is the wavenumber, :math:`v_A`
+    is the Alfvén velocity, :math:`c_s` is the sound speed, :math:`\omega_{ci}` is
+    the ion gyrofrequency, and :math:`\omega_{pe}` is the electron plasma frequency.
+    This relation does additionally assume low-frequency waves
+    :math:`\omega/kc \ll 1`, no D.C. electric field :math:`\mathbf{E_o}=0` and
+    quasi-neutrality.
+
+    Following section 5 of Bellan 2012 [1]_ the exact roots of the above dispersion
+    equation can be derived and expressed as one analytical solution (equation 38
+    of Bellan 2012 [1]_):
+
+    .. math::
+        \frac{\omega}{\omega_{ci}} = \sqrt{
+            2 \Lambda \sqrt{-\frac{P}{3}} \cos\left(
+                \frac{1}{3} \cos^{-1}\left(
+                    \frac{3q}{2p} \sqrt{-\frac{3}{p}}
+                \right)
+                - \frac{2 \pi}{3}j
+            \right)
+            + \frac{\Lambda A}{3}
+        }
+
+    where :math:`j = 0` represents the fast mode, :math:`j = 1` represents the
+    Alfvén mode, and :math:`j = 2` represents the acoustic mode.  Additionally,
+
+    .. math::
+        p &= \frac{3B-A^2}{3} \; , \; q = \frac{9AB-2A^3-27C}{27} \\
+        A &= \frac{Q + Q^2 \beta + Q \alpha + \alpha \Lambda}{Q^2} \;
+            , \; B = \alpha \frac{1 + 2 Q \beta + \Lambda \beta}{Q^2} \;
+            , \; C = \frac{\alpha^2 \beta}{Q^2} \\
+        \alpha &= \cos^2 \theta \;
+            , \; \beta = \left( \frac{c_s}{v_A}\right)^2 \;
+            , \; \Lambda = \left( \frac{k v_{A}}{\omega_{ci}}\right)^2
+
+    References
+    ----------
+    .. [1] PM Bellan, Improved basis set for low frequency plasma waves, 2012,
+       JGR, 117, A12219, doi: `10.1029/2012JA017856
+       <https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856>`_.
+
+    .. [2] TE Stringer, Low-frequency waves in an unbounded plasma, 1963, JNE,
+       Part C, doi: `10.1088/0368-3281/5/2/304
+       <https://doi.org/10.1088/0368-3281/5/2/304>`_
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from plasmapy.dispersion import two_fluid_dispersion
+    >>> inputs = {
+    ...     "k": 0.01 * u.rad / u.m,
+    ...     "theta": 30 * u.deg,
+    ...     "B": 8.3e-9 * u.T,
+    ...     "n_i": 5e6 * u.m ** -3,
+    ...     "T_e": 1.6e6 * u.K,
+    ...     "T_i": 4.0e5 * u.K,
+    ...     "ion": "p+",
+    ... }
+    >>> omegas = two_fluid_dispersion_solution(**inputs)
+    >>> omegas
+    {'fast_mode': <Quantity 1520.57... rad / s>,
+     'alfven_mode': <Quantity 1261.75... rad / s>,
+     'acoustic_mode': <Quantity 0.688152... rad / s>}
+
+    >>> inputs = {
+    ...     "k": [1e-7, 2e-7] * u.rad / u.m,
+    ...     "theta": [10, 20] * u.deg,
+    ...     "B": 8.3e-9 * u.T,
+    ...     "n_i": 5e6 * u.m ** -3,
+    ...     "T_e": 1.6e6 * u.K,
+    ...     "T_i": 4.0e5 * u.K,
+    ...     "ion": "He+",
+    ... }
+    >>> omegas = two_fluid_dispersion_solution(**inputs)
+    >>> omegas['fast_mode']
+    <Quantity [[0.00767..., 0.00779... ],
+               [0.01534..., 0.01558...]] rad / s>
+    """
+
+    # validate argument ion
+    if not isinstance(ion, Particle):
+        try:
+            ion = Particle(ion)
+        except TypeError:
+            raise TypeError(
+                f"For argument 'ion' expected type {Particle} but got {type(ion)}."
+            )
+    if not (ion.is_ion or ion.is_category("element")):
+        raise ValueError(f"The particle passed for 'ion' must be an ion or element.")
+
+    # validate z_mean
+    if z_mean is None:
+        try:
+            z_mean = abs(ion.integer_charge)
+        except ChargeError:
+            z_mean = 1
+    else:
+        if not isinstance(z_mean, (int, np.integer, float, np.floating)):
+            raise TypeError(
+                f"Expected int or float for argument 'z_mean', but got {type(z_mean)}."
+            )
+        z_mean = abs(z_mean)
+
+    # validate arguments
+    for arg_name in ("B", "n_i", "T_e", "T_i"):
+        val = locals()[arg_name].squeeze()
+        if val.shape != ():
+            raise ValueError(
+                f"Argument '{arg_name}' must a single value and not an array of "
+                f"shape {val.shape}."
+            )
+        locals()[arg_name] = val
+
+    # validate arguments
+    for arg_name in ("gamma_e", "gamma_i"):
+        if not isinstance(locals()[arg_name], (int, np.integer, float, np.floating)):
+            raise TypeError(
+                f"Expected int or float for argument '{arg_name}', but got "
+                f"{type(locals()[arg_name])}."
+            )
+
+    # validate argument k
+    k = k.squeeze()
+    if not (k.ndim == 0 or k.ndim == 1):
+        raise ValueError(
+            f"Argument 'k' needs to be a single valued or 1D array astropy Quantity,"
+            f" got array of shape {k.shape}."
+        )
+    if np.any(k <= 0):
+        raise ValueError(f"Argument 'k' can not be a or have negative values.")
+
+    # validate argument theta
+    theta = theta.squeeze()
+    theta = theta.to(u.radian)
+    if not (theta.ndim == 0 or theta.ndim == 1):
+        raise ValueError(
+            f"Argument 'theta' needs to be a single valued or 1D array astropy "
+            f"Quantity, got array of shape {k.shape}."
+        )
+
+    # Calc needed plasma parameters
+    n_e = z_mean * n_i
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", category=PhysicsWarning)
+        c_s = pfp.ion_sound_speed(
+            T_e=T_e,
+            T_i=T_i,
+            ion=ion,
+            n_e=n_e,
+            gamma_e=gamma_e,
+            gamma_i=gamma_i,
+            z_mean=z_mean,
+        )
+    v_A = pfp.Alfven_speed(B, n_i, ion=ion, z_mean=z_mean)
+    omega_ci = pfp.gyrofrequency(B=B, particle=ion, signed=False, Z=z_mean)
+    omega_pe = pfp.plasma_frequency(n=n_e, particle="e-")
+
+    # Bellan2012JGR params equation 32
+    alpha = np.cos(theta.value) ** 2
+    beta = (c_s / v_A).to(u.dimensionless_unscaled).value ** 2
+    alphav, kv = np.meshgrid(alpha, k.value)  # create grid
+    Lambda = (kv * v_A.value / omega_ci.value) ** 2
+
+    # Bellan2012JGR params equation 2
+    Q = 1 + (kv * c.value / omega_pe.value) ** 2
+
+    # Bellan2012JGR params equation 35
+    A = ((1 + alphav) / Q) + beta + (alphav * Lambda / Q ** 2)
+    B = alphav * (1 + 2 * Q * beta + Lambda * beta) / Q ** 2
+    C = beta * (alphav / Q) ** 2
+
+    # Bellan2012JGR params equation 36
+    p = (3 * B - A ** 2) / 3
+    q = (9 * A * B - 2 * A ** 3 - 27 * C) / 27
+
+    # Bellan2012JGR params equation 38
+    R = 2 * Lambda * np.emath.sqrt(-p / 3)
+    S = 3 * q / (2 * p) * np.emath.sqrt(-3 / p)
+    T = Lambda * A / 3
+    omega = {}
+    for ind, key in enumerate(("fast_mode", "alfven_mode", "acoustic_mode")):
+        # The solution corresponding to equation 38
+        w = omega_ci * np.emath.sqrt(
+            R * np.cos(1 / 3 * np.emath.arccos(S) - 2 * np.pi / 3 * ind) + T
+        )
+        omega[key] = w.squeeze()
+
+        # check for violation of dispersion relation assumptions
+        # (i.e. low-frequency, w/kc << 0.1)
+        wkc_max = np.max(w.value / (kv * c.value))
+        if wkc_max > 0.1:
+            warnings.warn(
+                f"The {key} calculation produced a high-frequency wave (w/kc == "
+                f"{wkc_max:.3f}), which violates the low-frequency (w/kc << 1) "
+                f"assumption of the dispersion relation.",
+                PhysicsWarning,
+            )
+
+    return omega
+
+
+tfds_ = two_fluid_dispersion_solution
+""" Alias to :func:`two_fluid_dispersion_solution`. """
diff --git a/plasmapy/formulary/braginskii.py b/plasmapy/formulary/braginskii.py
index 73c51b24..fd42440f 100644
--- a/plasmapy/formulary/braginskii.py
+++ b/plasmapy/formulary/braginskii.py
@@ -254,6 +254,16 @@ class can be initialized once with all of the variables necessary for
         on the ion transport coefficients. Defaults to T_e / T_i. Only
         has effect if mu is non-zero.
 
+    coulomb_log_method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of
+        `~plasmapy.formulary.collisions.Coulomb_logarithm` for more
+        information about these methods.
+
     Raises
     ------
     ValueError
diff --git a/plasmapy/formulary/collisions.py b/plasmapy/formulary/collisions.py
index 08d1a5c1..eedd76a5 100644
--- a/plasmapy/formulary/collisions.py
+++ b/plasmapy/formulary/collisions.py
@@ -16,9 +16,8 @@
 Coulomb logarithms
 ==================
 
-Please see the documentation for the `Coulomb_logarithm`
-for a review of the many ways in which one can define and calculate
-that quantity.
+Please read the documentation of `Coulomb_logarithm` below for an explanation of the
+7 PlasmaPy-supported methods of computing the Coulomb logarithm.
 
 Collision rates
 ===============
@@ -89,152 +88,288 @@ def Coulomb_logarithm(
     method="classical",
 ):
     r"""
-    Estimates the Coulomb logarithm.
+    Computes the Coulomb logarithm.
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
         the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
         (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the Notes section of this docstring for more
+        information, including about abbreviated aliases of these names.
 
     Returns
     -------
-    lnLambda : float or numpy.ndarray
-        An estimate of the Coulomb logarithm that is accurate to
-        roughly its reciprocal.
+    ln_Lambda : float or numpy.ndarray
+        The dimensionless Coulomb logarithm.
 
     Raises
     ------
-    ValueError
+    `ValueError`
         If the mass or charge of either particle cannot be found, or
         any of the inputs contain incorrect values.
 
-    UnitConversionError
+    `~astropy.units.UnitConversionError`
         If the units on any of the inputs are incorrect.
 
         If the n_e, T, or V are not Quantities.
 
-    PhysicsError
+    `~plasmapy.utils.exceptions.PhysicsError`
         If the result is smaller than 1.
 
-    RelativityError
+    `~plasmapy.utils.exceptions.RelativityError`
         If the input velocity is same or greater than the speed
         of light.
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
-        If units are not provided, SI units are assumed
+    : `~astropy.units.UnitsWarning`
+        If units are not provided, SI units are assumed.
 
-    ~plasmapy.utils.RelativityWarning
+    : `~plasmapy.utils.exceptions.RelativityWarning`
         If the input velocity is greater than 5% of the speed of
         light.
 
     Notes
     -----
-    The classical Coulomb logarithm is given by
+    **Summary of Supported Methods of Computing the Coulomb Logarithm**
+
+    PlasmaPy supports 7 methods of computing the Coulomb logarithm:
+
+    1. ``"classical"`` or ``"ls"``
+    2. ``"ls_min_interp"`` or ``"GMS-1"``
+    3. ``"ls_full_interp"`` or ``"GMS-2"``
+    4. ``"ls_clamp_mininterp"`` or ``"GMS-3"``
+    5. ``"hls_min_interp"`` or ``"GMS-4"``
+    6. ``"hls_max_interp"`` or ``"GMS-5"``
+    7. ``"hls_full_interp"`` or ``"GMS-6"``
+
+    Options 1–4 are straight-line Landau-Spitzer (``"ls..."``) methods in which the trajectory of a
+    Coulomb collision is modeled as a straight line. For the straight-line Landau-Spitzer methods, the Coulomb
+    logarithm (:math:`\ln{\Lambda}`) is defined to be:
 
     .. math::
         \ln{\Lambda} \equiv \ln\left( \frac{b_{max}}{b_{min}} \right)
 
-    where :math:`b_{min}` and :math:`b_{max}` are the inner and outer
-    impact parameters for Coulomb collisions [1]_.
+    Options 5–7 are hyperbolic Landau-Spitzer (``"hls..."``) methods in which the trajectory of a
+    Coulomb collision is modeled as a hyperbola. For the hyperbolic Landau-Spitzer methods, the Coulomb
+    logarithm (:math:`\ln{\Lambda}`) is defined to be:
 
-    The outer impact parameter is given by the Debye length:
-    :math:`b_{min} = \lambda_D` which is a function of electron
-    temperature and electron density.  At distances greater than the
-    Debye length, electric fields from other particles will be
-    screened out due to electrons rearranging themselves.
+    .. math::
+        \ln{\Lambda} \equiv \frac{1}{2} \ln\left(1 + \frac{b_{max}^2}{b_{min}^2} \right)
 
-    The choice of inner impact parameter is either the distance of closest
-    approach for a 90 degree Coulomb collision or the thermal deBroglie
-    wavelength, whichever is larger. This is because Coulomb-style collisions
-    cannot occur for impact parameters shorter than the deBroglie
-    wavelength because quantum effects will change the fundamental
-    nature of the collision [2]_, [3]_.
+    For all 7 methods, :math:`b_{min}` and :math:`b_{max}` are the inner impact parameter and the outer
+    impact parameter, respectively, for Coulomb collisions [1]_;
+    :math:`b_{min}` and :math:`b_{max}` are each computed by `impact_parameter`, another function.
 
-    Errors associated with the classical Coulomb logarithm are of order its
-    inverse. If the Coulomb logarithm is of order unity, then the
-    assumptions made in the standard analysis of Coulomb collisions
-    are invalid.
+    The abbreviations of Options 2-7 (``"GMS-..."``) refer to the first initials of the three authors
+    of Reference [4]_.
 
-    For dense plasmas where the classical Coulomb logarithm breaks
-    down there are various extended methods. These can be found
-    in D.O. Gericke et al's paper, which has a table summarizing
-    the methods [4]_. The GMS-1 through GMS-6 methods correspond
-    to the methods found it that table.
+    .. note::
+        For strongly-coupled plasma, PlasmaPy recommends Option 7, ``"hls_full_interp"`` or ``"GMS-6"``,
+        because of its high accuracy regardless of a plasma's strength of coupling.
 
-    It should be noted that GMS-4 thru GMS-6 modify the Coulomb
-    logarithm to the form:
+    **Explanation of Supported Methods of Computing the Coulomb Logarithm**
 
-    .. math::
-        \ln{\Lambda} \equiv 0.5 \ln\left(1 + \frac{b_{max}^2}{b_{min}^2} \right)
-
-    This means the Coulomb logarithm will not break down for Lambda < 0,
-    which occurs for dense, cold plasmas.
-
-    Methods
-    ---
-    Classical
-        classical Landau-Spitzer approach. Fails for large coupling
-        parameter where Lambda can become less than zero.
-    GMS-1
-        1st method listed in Table 1 of reference [3]
-        Landau-Spitzer, but with interpolated bmin instead of bmin
-        selected between deBroglie wavelength and distance of closest
-        approach. Fails for large coupling
-        parameter where Lambda can become less than zero.
-    GMS-2
-        2nd method listed in Table 1 of reference [3]
-        Another Landau-Spitzer like approach, but now bmax is also
-        being interpolated. The interpolation is between the Debye
-        length and the ion sphere radius, allowing for descriptions
-        of dilute plasmas. Fails for large coupling
-        parameter where Lambda can become less than zero.
-        3rd method listed in Table 1 of reference [3]
-        classical Landau-Spitzer fails for argument of Coulomb logarithm
-        Lambda < 0, therefore a clamp is placed at Lambda_min = 2
-    GMS-4
-        4th method listed in Table 1 of reference [3]
-        Spitzer-like extension to Coulomb logarithm by noting that
-        Coulomb collisions take hyperbolic trajectories. Removes
-        divergence for small bmin issue in classical Landau-Spitzer
-        approach, so bmin can be zero. Also doesn't break down as
-        Lambda < 0 is now impossible, even when coupling parameter is large.
-    GMS-5
-        5th method listed in Table 1 of reference [3]
-        Similar to GMS-4, but setting bmin as distance of closest approach
-        and bmax interpolated between Debye length and ion sphere radius.
-        Lambda < 0 impossible.
-    GMS-6
-        6th method listed in Table 1 of reference [3]
-        Similar to GMS-4 and GMS-5, but using interpolation methods
-        for both bmin and bmax.
+    In this section, further information about each method, such as about
+    interpolation and other special features, is documented. Please refer
+    to Reference [1]_ and Reference [4]_ for additional information about
+    these methods.
+
+    Option 1: ``"classical"`` or ``"ls"`` (Landau-Spitzer)
+        The classical straight-line Landau-Spitzer method in which :math:`b_{min}` is defined to be the
+        higher of the de Broglie wavelength (:math:`\lambda_{de Broglie}`) and the distance of closest
+        approach (:math:`\rho_{\perp}`) if they are not equal (and either of the two if they are equal) and
+        :math:`b_{max}` is defined to be the Debye length (:math:`\lambda_{Debye}`).
+
+        .. math::
+            \ln{\Lambda} \equiv \ln\left( \frac{b_{max}}{b_{min}} \right)
+
+        .. math::
+            b_{min} \equiv
+            \left\{
+                \begin{array}{ll}
+                           \lambda_{de Broglie} & \mbox{if } \lambda_{de Broglie} \geq \rho_{\perp} \\
+                           \rho_{\perp}         & \mbox{if } \rho_{\perp} \geq \lambda_{de Broglie}
+                \end{array}
+            \right.
+
+        .. math::
+            b_{max} \equiv \lambda_{Debye}
+
+        The inner impact parameter (:math:`b_{min}`) is the higher of :math:`\lambda_{de Broglie}`
+        and :math:`\rho_{\perp}` because for impact parameters lower than :math:`\lambda_{de Broglie}`,
+        quantum effects cause the collision to be non-Coulombic [2]_ [3]_.
+
+        The outer impact parameter (:math:`b_{max}`) is defined to be the Debye length
+        (:math:`\lambda_{Debye}`) because at distances higher than the
+        Debye length, the electric fields created by other particles are
+        screened out by the electrons rearranging themselves.
+
+        The uncertainty of the classical straight-line Landau-Spitzer method is on the order
+        of its reciprocal.
+
+        This method is invalid if :math:`\ln{\Lambda} < 2` because of the uncertainty of this method
+        and is invalid if :math:`\ln{\Lambda} < 0`, which may be true if the
+        coupling parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Please refer to Reference [1]_ for additional information about this method.
+
+    Option 2: ``"ls_min_interp"`` or ``"GMS-1"`` (Landau-Spitzer, interpolation of :math:`b_{min}`)
+        A straight-line Landau-Spitzer method in which :math:`b_{min}` is interpolated between the
+        de Broglie wavelength (:math:`\lambda_{de Broglie}`) and the distance of closest approach
+        (:math:`\rho_{\perp}`) and :math:`b_{max}` is defined to be the Debye length (:math:`\lambda_{Debye}`).
+
+        .. math::
+            \ln{\Lambda} \equiv \ln\left( \frac{b_{max}}{b_{min}} \right)
+
+        .. math::
+            b_{min} \equiv \sqrt{\lambda_{de Broglie}^2 + \rho_{\perp}^2}
+
+        .. math::
+            b_{max} \equiv \lambda_{Debye}
+
+        This method is invalid if :math:`\ln{\Lambda} < 0`, which may be true if the coupling
+        parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Note: This is the first method in Table 1 of Reference [4]_.
+
+    Option 3: ``"ls_full_interp"`` or ``"GMS-2"`` (Landau-Spitzer, interpolation of :math:`b_{min}` and :math:`b_{max}`)
+        A straight-line Landau-Spitzer method in which :math:`b_{min}` and :math:`b_{max}`
+        are each interpolated. :math:`b_{min}` is interpolated between the de Broglie wavelength
+        (:math:`\lambda_{de Broglie}`) and the distance of closest approach (:math:`\rho_{\perp}`).
+        :math:`b_{max}` is interpolated between the Debye length (:math:`\lambda_{Debye}`)
+        and the ion sphere radius (:math:`a_i`).
+
+        .. math::
+            \ln{\Lambda} \equiv \ln\left( \frac{b_{max}}{b_{min}} \right)
+
+        .. math::
+            b_{min} \equiv \sqrt{\lambda_{de Broglie}^2 + \rho_{\perp}^2}
+
+        .. math::
+            b_{max} \equiv \sqrt{\lambda_{Debye}^2 + a_i^2}
+
+        This method is invalid if :math:`\ln{\Lambda} < 0`, which may be true if the coupling
+        parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Note: This is the second method in Table 1 of Reference [4]_.
+
+    Option 4: ``"ls_clamp_mininterp"`` or ``"GMS-3"`` (Landau-Spitzer with a clamp, interpolation of :math:`b_{min}`)
+        A straight-line Landau-Spitzer method in which the value of :math:`\ln{\Lambda}` is clamped at
+        a minimum of :math:`2` so that it is impossible for :math:`\ln{\Lambda} < 0` (which is possible by the
+        classical Landau-Spitzer method). :math:`b_{min}` is interpolated between the de Broglie
+        wavelength (:math:`\lambda_{de Broglie}`) and the distance of closest approach
+        (:math:`\rho_{\perp}`). :math:`b_{max}` is defined to be the Debye length (:math:`\lambda_{Debye}`).
+
+        .. math::
+            \ln{\Lambda} \equiv
+            \left\{
+                \begin{array}{ll}
+                           \ln\left( \frac{b_{max}}{b_{min}} \right) & \mbox{if } \ln\left( \frac{b_{max}}{b_{min}} \right) \geq 2 \\
+                           2                                         & \mbox{if } \ln\left( \frac{b_{max}}{b_{min}} \right) \leq 2
+                \end{array}
+            \right.
+
+        .. math::
+            b_{min} \equiv \sqrt{\lambda_{de Broglie}^2 + \rho_{\perp}^2}
+
+        .. math::
+            b_{max} \equiv \lambda_{Debye}
+
+        This method is valid for any plasma because it is impossible for :math:`\ln{\Lambda} < 0` by this
+        method, even if the coupling parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Note: This is the third method in Table 1 of Reference [4]_.
+
+    Option 5: ``"hls_min_interp"`` or ``"GMS-4"`` (Hyperbolic Landau-Spitzer, interpolation of :math:`b_{min}`)
+        A hyperbolic Landau-Spitzer method in which :math:`b_{min}` is interpolated between the
+        de Broglie wavelength (:math:`\lambda_{de Broglie}`) and the distance of closest approach
+        (:math:`\rho_{\perp}`) and :math:`b_{max}` is defined to be the Debye length (:math:`\lambda_{Debye}`).
+
+        .. math::
+            \ln{\Lambda} \equiv \frac{1}{2} \ln\left(1 + \frac{b_{max}^2}{b_{min}^2} \right)
+
+        .. math::
+            b_{min} \equiv \sqrt{\lambda_{de Broglie}^2 + \rho_{\perp}^2}
+
+        .. math::
+            b_{max} \equiv \lambda_{Debye}
+
+        This method is valid for any plasma because it is impossible for :math:`\ln{\Lambda} < 0` by this
+        method, even if the coupling parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Note: This is the fourth method in Table 1 of Reference [4]_.
+
+    Option 6: ``"hls_max_interp"`` or ``"GMS-5"`` (Hyperbolic Landau-Spitzer, interpolation of :math:`b_{max}`)
+        A hyperbolic Landau-Spitzer method in which :math:`b_{max}` is interpolated between
+        the Debye length (:math:`\lambda_{Debye}`) and the ion sphere radius (:math:`a_i`)
+        and :math:`b_{min}` is defined to be the distance of closest approach (:math:`\rho_{\perp}`).
+
+        .. math::
+            \ln{\Lambda} \equiv \frac{1}{2} \ln\left(1 + \frac{b_{max}^2}{b_{min}^2} \right)
+
+        .. math::
+            b_{min} \equiv \rho_{\perp}
+
+        .. math::
+            b_{max} \equiv \sqrt{\lambda_{Debye}^2 + a_i^2}
+
+        This method is valid for any plasma because it is impossible for :math:`\ln{\Lambda} < 0` by this
+        method, even if the coupling parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Caution: This method overestimates :math:`\ln{\Lambda}` at high temperatures.
+
+        Note: This is the fifth method in Table 1 of Reference [4]_.
+
+    Option 7: ``"hls_full_interp"`` or ``"GMS-6"`` (Hyperbolic Landau-Spitzer, interpolation of :math:`b_{min}` and :math:`b_{max}`)
+        A hyperbolic Landau-Spitzer method in which :math:`b_{min}` and :math:`b_{max}`
+        are each interpolated. :math:`b_{min}` is interpolated between the de Broglie wavelength
+        (:math:`\lambda_{de Broglie}`) and the distance of closest approach (:math:`\rho_{\perp}`).
+        :math:`b_{max}` is interpolated between the Debye length (:math:`\lambda_{Debye}`)
+        and the ion sphere radius (:math:`a_i`).
+
+        .. math::
+            \ln{\Lambda} \equiv \frac{1}{2} \ln\left(1 + \frac{b_{max}^2}{b_{min}^2} \right)
+
+        .. math::
+            b_{min} \equiv \sqrt{\lambda_{de Broglie}^2 + \rho_{\perp}^2}
+
+        .. math::
+            b_{max} \equiv \sqrt{\lambda_{Debye}^2 + a_i^2}
+
+        This method is valid for any plasma because it is impossible for :math:`\ln{\Lambda} < 0` by this
+        method, even if the coupling parameter is high (such as for nonideal, dense, cold plasmas).
+
+        Note: This is the sixth method in Table 1 of Reference [4]_.
 
     Examples
     --------
@@ -261,42 +396,70 @@ def Coulomb_logarithm(
     .. [4] Dense plasma temperature equilibration in the binary collision
        approximation. D. O. Gericke et. al. PRE,  65, 036418 (2002).
        DOI: 10.1103/PhysRevE.65.036418
+
+    See Also
+    --------
+    impact_parameter : Computes :math:`b_{min}` and :math:`b_{max}`.
     """
     # fetching impact min and max impact parameters
     bmin, bmax = impact_parameter(
         T=T, n_e=n_e, species=species, z_mean=z_mean, V=V, method=method
     )
-    if method in ("classical", "GMS-1", "GMS-2"):
+
+    if method in (
+        "classical",
+        "ls",
+        "ls_min_interp",
+        "GMS-1",
+        "ls_full_interp",
+        "GMS-2",
+    ):
         ln_Lambda = np.log(bmax / bmin)
-    elif method == "GMS-3":
+    elif method in ("ls_clamp_mininterp", "GMS-3"):
         ln_Lambda = np.log(bmax / bmin)
         if np.any(ln_Lambda < 2):
             if np.isscalar(ln_Lambda.value):
                 ln_Lambda = 2 * u.dimensionless_unscaled
             else:
                 ln_Lambda[ln_Lambda < 2] = 2 * u.dimensionless_unscaled
-    elif method in ("GMS-4", "GMS-5", "GMS-6"):
+    elif method in (
+        "hls_min_interp",
+        "GMS-4",
+        "hls_max_interp",
+        "GMS-5",
+        "hls_full_interp",
+        "GMS-6",
+    ):
         ln_Lambda = 0.5 * np.log(1 + bmax ** 2 / bmin ** 2)
     else:
         raise ValueError(
-            "Unknown method! Choose from 'classical' and 'GMS-N', N from 1 to 6."
+            'Unknown method. Choose from "classical", "ls_min_interp", "ls_full_interp", "ls_clamp_mininterp", "hls_min_interp", "hls_max_interp", "hls_full_interp", and their aliases. Please refer to the documentation of this function for more information.'
         )
+
     # applying dimensionless units
     ln_Lambda = ln_Lambda.to(u.dimensionless_unscaled).value
+
     # Allow NaNs through the < checks without warning
     with np.errstate(invalid="ignore"):
-        if np.any(ln_Lambda < 2) and method in ["classical", "GMS-1", "GMS-2"]:
+        if np.any(ln_Lambda < 2) and method in [
+            "classical",
+            "ls",
+            "ls_min_interp",
+            "GMS-1",
+            "ls_full_interp",
+            "GMS-2",
+        ]:
             warnings.warn(
-                f"Coulomb logarithm is {ln_Lambda} and {method} relies on "
-                "weak coupling.",
+                f'The Coulomb logarithm is {ln_Lambda}, and the specified method, "{method}", depends on weak coupling.',
                 utils.CouplingWarning,
             )
         elif np.any(ln_Lambda < 4):
             warnings.warn(
-                f"Coulomb logarithm is {ln_Lambda}, you might have strong "
-                "coupling effects",
+                f"The Coulomb logarithm is {ln_Lambda}, so strong "
+                "coupling effects may exist for the plasma.",
                 utils.CouplingWarning,
             )
+
     return ln_Lambda
 
 
@@ -330,6 +493,7 @@ def _replaceNanVwithThermalV(V, T, m):
     if np.any(V == 0):
         raise utils.PhysicsError("You cannot have a collision for zero velocity!")
     # getting thermal velocity of system if no velocity is given
+
     if V is None:
         V = parameters.thermal_speed(T, "e-", mass=m)
     elif np.any(np.isnan(V)):
@@ -341,6 +505,7 @@ def _replaceNanVwithThermalV(V, T, m):
         else:
             V = V.copy()
             V[np.isnan(V)] = parameters.thermal_speed(T[np.isnan(V)], "e-", mass=m)
+
     return V
 
 
@@ -357,7 +522,6 @@ def impact_parameter_perp(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
@@ -395,10 +559,10 @@ def impact_parameter_perp(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -409,7 +573,6 @@ def impact_parameter_perp(
     .. math::
         b_{\perp} = \frac{Z_1 Z_2}{4 \pi \epsilon_0 m v^2}
 
-
     Examples
     --------
     >>> from astropy import units as u
@@ -418,12 +581,10 @@ def impact_parameter_perp(
     >>> impact_parameter_perp(T, species)
     <Quantity 8.3550...e-12 m>
 
-
     References
     ----------
     .. [1] Francis, F. Chen. Introduction to plasma physics and controlled
        fusion 3rd edition. Ch 5 (Springer 2015).
-
     """
     # boiler plate checks
     T, masses, charges, reduced_mass, V = _boilerPlate(T=T, species=species, V=V)
@@ -453,34 +614,39 @@ def impact_parameter(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
-        the target particle
+        the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -506,10 +672,10 @@ def impact_parameter(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -526,13 +692,13 @@ def impact_parameter(
         b_{max} = \left(\lambda_{De}^2 + a_i^2\right)^{1/2}
 
     The minimum impact parameter is typically some combination of the
-    thermal deBroglie wavelength and the distance of closest approach
+    thermal de Broglie wavelength and the distance of closest approach
     for a 90 degree Coulomb collision. A quadratic sum is used for
-    all GMS methods, except for GMS-5, where b_min is simply set to
+    all GMS methods, except for GMS-5, where ``b_min`` is simply set to
     the distance of closest approach [1]_.
 
     .. math::
-        b_{min} = \left(\Lambda_{deBroglie}^2 + \rho_{\perp}^2\right)^{1/2}
+        b_{min} = \left(\Lambda_{de Broglie}^2 + \rho_{\perp}^2\right)^{1/2}
 
     Examples
     --------
@@ -555,25 +721,33 @@ def impact_parameter(
     T, masses, charges, reduced_mass, V = _boilerPlate(T=T, species=species, V=V)
     # catching error where mean charge state is not given for non-classical
     # methods that require the ion density
-    if method in ("GMS-2", "GMS-5", "GMS-6"):
+    if method in (
+        "ls_full_interp",
+        "GMS-2",
+        "hls_max_interp",
+        "GMS-5",
+        "hls_full_interp",
+        "GMS-6",
+    ):
         if np.isnan(z_mean):
             raise ValueError(
-                "Must provide a z_mean for GMS-2, GMS-5, and GMS-6 methods."
+                'Must provide a z_mean for "ls_full_interp", "hls_max_interp", and "hls_full_interp" methods.'
             )
     # Debye length
     lambdaDe = parameters.Debye_length(T, n_e)
-    # deBroglie wavelength
+    # de Broglie wavelength
     lambdaBroglie = hbar / (2 * reduced_mass * V)
     # distance of closest approach in 90 degree Coulomb collision
     bPerp = impact_parameter_perp(T=T, species=species, V=V)
+
     # obtaining minimum and maximum impact parameters depending on which
     # method is requested
-    if method == "classical":
+    if method == "classical" or method == "ls":
         bmax = lambdaDe
         # Coulomb-style collisions will not happen for impact parameters
         # shorter than either of these two impact parameters, so we choose
         # the larger of these two possibilities. That is, between the
-        # deBroglie wavelength and the distance of closest approach.
+        # de Broglie wavelength and the distance of closest approach.
         # ARRAY NOTES
         # T and V should be guaranteed to be same size inputs from _boilerplate
         # therefore, lambdaBroglie and bPerp are either both scalar or both array
@@ -587,14 +761,14 @@ def impact_parameter(
         except ValueError:  # both lambdaBroglie and bPerp are arrays
             bmin = lambdaBroglie
             bmin[bPerp > lambdaBroglie] = bPerp[bPerp > lambdaBroglie]
-    elif method == "GMS-1":
+    elif method == "ls_min_interp" or method == "GMS-1":
         # 1st method listed in Table 1 of reference [1]
         # This is just another form of the classical Landau-Spitzer
-        # approach, but bmin is interpolated between the deBroglie
+        # approach, but bmin is interpolated between the de Broglie
         # wavelength and distance of closest approach.
         bmax = lambdaDe
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
-    elif method == "GMS-2":
+    elif method == "ls_full_interp" or method == "GMS-2":
         # 2nd method listed in Table 1 of reference [1]
         # Another Landau-Spitzer like approach, but now bmax is also
         # being interpolated. The interpolation is between the Debye
@@ -606,17 +780,17 @@ def impact_parameter(
         ionRadius = Wigner_Seitz_radius(n_i)
         bmax = (lambdaDe ** 2 + ionRadius ** 2) ** (1 / 2)
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
-    elif method == "GMS-3":
+    elif method == "ls_clamp_mininterp" or method == "GMS-3":
         # 3rd method listed in Table 1 of reference [1]
         # same as GMS-1, but not Lambda has a clamp at Lambda_min = 2
         # where Lambda is the argument to the Coulomb logarithm.
         bmax = lambdaDe
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
-    elif method == "GMS-4":
+    elif method == "hls_min_interp" or method == "GMS-4":
         # 4th method listed in Table 1 of reference [1]
         bmax = lambdaDe
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
-    elif method == "GMS-5":
+    elif method == "hls_max_interp" or method == "GMS-5":
         # 5th method listed in Table 1 of reference [1]
         # Mean ion density.
         n_i = n_e / z_mean
@@ -624,7 +798,7 @@ def impact_parameter(
         ionRadius = Wigner_Seitz_radius(n_i)
         bmax = (lambdaDe ** 2 + ionRadius ** 2) ** (1 / 2)
         bmin = bPerp
-    elif method == "GMS-6":
+    elif method == "hls_full_interp" or method == "GMS-6":
         # 6th method listed in Table 1 of reference [1]
         # Mean ion density.
         n_i = n_e / z_mean
@@ -634,6 +808,7 @@ def impact_parameter(
         bmin = (lambdaBroglie ** 2 + bPerp ** 2) ** (1 / 2)
     else:
         raise ValueError(f"Method {method} not found!")
+
     # ARRAY NOTES
     # it could be that bmin and bmax have different sizes. If Te is a scalar,
     # T and V will be scalar from _boilerplate, so bmin will scalar.  However
@@ -641,6 +816,7 @@ def impact_parameter(
     # do we want to extend the scalar bmin to equal the length of bmax? Sure.
     if np.isscalar(bmin.value) and not np.isscalar(bmax.value):
         bmin = np.repeat(bmin, len(bmax))
+
     return bmin.to(u.m), bmax.to(u.m)
 
 
@@ -661,14 +837,12 @@ def collision_frequency(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature.
         This should be the electron temperature for electron-electron
         and electron-ion collisions, and the ion temperature for
         ion-ion collisions.
 
-
     n : ~astropy.units.Quantity
         The density in units convertible to per cubic meter.
         This should be the electron density for electron-electron collisions,
@@ -676,23 +850,29 @@ def collision_frequency(
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -717,10 +897,10 @@ def collision_frequency(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -736,7 +916,7 @@ def collision_frequency(
     taken as the thermal velocity), and :math:`\ln{\Lambda}` is the Coulomb
     logarithm accounting for small angle collisions.
 
-    See eq (2.14) in [2]_.
+    See Equation (2.14) in [2]_.
 
     Examples
     --------
@@ -752,13 +932,13 @@ def collision_frequency(
     .. [1] Francis, F. Chen. Introduction to plasma physics and controlled
        fusion 3rd edition. Ch 5 (Springer 2015).
     .. [2] http://homepages.cae.wisc.edu/~callen/chap2.pdf
-
     """
     # boiler plate checks
     T, masses, charges, reduced_mass, V_r = _boilerPlate(T=T, species=species, V=V)
     # using a more descriptive name for the thermal velocity using
     # reduced mass
     V_reduced = V_r
+
     if species[0] in ("e", "e-") and species[1] in ("e", "e-"):
         # electron-electron collision
         # if a velocity was passed, we use that instead of the reduced
@@ -791,6 +971,7 @@ def collision_frequency(
         bPerp = impact_parameter_perp(T=T, species=species, V=V)
         # Coulomb logarithm
         cou_log = Coulomb_logarithm(T, n, species, z_mean, V=V, method=method)
+
     # collisional cross section
     sigma = Coulomb_cross_section(bPerp)
     # collision frequency where Coulomb logarithm accounts for
@@ -809,12 +990,10 @@ def Coulomb_cross_section(impact_param: u.m) -> u.m ** 2:
     impact_param : ~astropy.units.Quantity
         Impact parameter for the collision.
 
-    Examples
-    --------
-    >>> Coulomb_cross_section(7e-10*u.m)
-    <Quantity 6.157...e-18 m2>
-    >>> Coulomb_cross_section(0.5*u.m)
-    <Quantity 3.141... m2>
+    Returns
+    -------
+    sigma : ~astropy.units.Quantity
+        The Coulomb collision cross section area.
 
     Notes
     -----
@@ -829,10 +1008,12 @@ def Coulomb_cross_section(impact_param: u.m) -> u.m ** 2:
     calculation. Please note that it is not guaranteed to return the correct
     results for small angle collisions.
 
-    Returns
-    -------
-    ~astropy.units.Quantity
-        The Coulomb collision cross section area.
+    Examples
+    --------
+    >>> Coulomb_cross_section(7e-10*u.m)
+    <Quantity 6.157...e-18 m2>
+    >>> Coulomb_cross_section(0.5*u.m)
+    <Quantity 3.141... m2>
 
     References
     ----------
@@ -871,7 +1052,7 @@ def fundamental_electron_collision_freq(
     n_e : ~astropy.units.Quantity
         The number density of the Maxwellian test electrons
 
-    ion: str
+    ion : str
         String signifying a particle type of the field ions, including charge
         state information.
 
@@ -885,9 +1066,18 @@ def fundamental_electron_collision_freq(
         If not specified, the Coulomb log will is calculated using the
         `~plasmapy.formulary.Coulomb_logarithm` function.
 
-    coulomb_log_method : string, optional
-        Method used for Coulomb logarithm calculation (see that function
-        for more documentation). Choose from "classical" or "GMS-1" to "GMS-6".
+    coulomb_log_method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
+
+    Returns
+    -------
+    nu_e : ~astropy.units.Quantity
 
     Notes
     -----
@@ -991,7 +1181,7 @@ def fundamental_ion_collision_freq(
     n_i : ~astropy.units.Quantity
         The number density of the Maxwellian test ions
 
-    ion: str
+    ion : str
         String signifying a particle type of the test and field ions,
         including charge state information. This function assumes the test
         and field ions are the same species.
@@ -1007,8 +1197,17 @@ def fundamental_ion_collision_freq(
         ~plasmapy.formulary.Coulomb_logarithm function.
 
     coulomb_log_method : str, optional
-        Method used for Coulomb logarithm calculation (see that function
-        for more documentation). Choose from "classical" or "GMS-1" to "GMS-6".
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
+
+    Returns
+    -------
+    nu_i : ~astropy.units.Quantity
 
     Notes
     -----
@@ -1033,7 +1232,6 @@ def fundamental_ion_collision_freq(
     * 1/tau_i from ref [1]_, equation (2.5i) pp. 215,
     * nu_i from ref [2]_ pp. 33,
 
-
     References
     ----------
     .. [1] Braginskii, S. I. "Transport processes in a plasma." Reviews of
@@ -1115,34 +1313,39 @@ def mean_free_path(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
-        the target particle
+        the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -1167,10 +1370,10 @@ def mean_free_path(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -1231,7 +1434,6 @@ def Spitzer_resistivity(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature.
         This should be the electron temperature for electron-electron
@@ -1244,24 +1446,30 @@ def Spitzer_resistivity(
         and the ion density for electron-ion and ion-ion collisions.
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -1286,10 +1494,10 @@ def Spitzer_resistivity(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -1325,7 +1533,6 @@ def Spitzer_resistivity(
     .. [1] Francis, F. Chen. Introduction to plasma physics and controlled
        fusion 3rd edition. Ch 5 (Springer 2015).
     .. [2] http://homepages.cae.wisc.edu/~callen/chap2.pdf
-
     """
     # collisional frequency
     freq = collision_frequency(
@@ -1358,37 +1565,42 @@ def mobility(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
-        the target particle
+        the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
         impact parameters. It is also used the obtain the average mobility
         of a plasma with multiple charge state species. When z_mean
         is not given, the average charge between the two particles is
-        used instead.
+        used instead. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -1413,10 +1625,10 @@ def mobility(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -1485,37 +1697,42 @@ def Knudsen_number(
 
     Parameters
     ----------
-
     characteristic_length : ~astropy.units.Quantity
         Rough order-of-magnitude estimate of the relevant size of the system.
 
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
-        the target particle
+        the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -1540,10 +1757,10 @@ def Knudsen_number(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -1578,7 +1795,6 @@ def Knudsen_number(
     References
     ----------
     .. [1] https://en.wikipedia.org/wiki/Knudsen_number
-
     """
     path_length = mean_free_path(
         T=T, n_e=n_e, species=species, z_mean=z_mean, V=V, method=method
@@ -1607,34 +1823,39 @@ def coupling_parameter(
 
     Parameters
     ----------
-
     T : ~astropy.units.Quantity
         Temperature in units of temperature or energy per particle,
         which is assumed to be equal for both the test particle and
-        the target particle
+        the target particle.
 
     n_e : ~astropy.units.Quantity
-        The electron density in units convertible to per cubic meter.
+        The electron number density in units convertible to per cubic meter.
 
     species : tuple
         A tuple containing string representations of the test particle
-        (listed first) and the target particle (listed second)
+        (listed first) and the target particle (listed second).
 
     z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. This is used to recover the
+        The average ionization (arithmetic mean) of a plasma for which
+        a macroscopic description is valid. This parameter is used to compute the
         average ion density (given the average ionization and electron
         density) for calculating the ion sphere radius for non-classical
-        impact parameters.
+        impact parameters. ``z_mean`` is a required parameter if ``method`` is
+        ``"ls_full_interp"``, ``"hls_max_interp"``, or ``"hls_full_interp"``.
 
     V : ~astropy.units.Quantity, optional
-        The relative velocity between particles.  If not provided,
+        The relative velocity between particles. If not provided,
         thermal velocity is assumed: :math:`\mu V^2 \sim 2 k_B T`
         where `mu` is the reduced mass.
 
-    method: str, optional
-        Selects which theory to use when calculating the Coulomb
-        logarithm. Defaults to classical method.
+    method : str, optional
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of `Coulomb_logarithm` for more
+        information about these methods.
 
     Returns
     -------
@@ -1659,10 +1880,10 @@ def coupling_parameter(
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : ~astropy.units.UnitsWarning
         If units are not provided, SI units are assumed
 
-    ~plasmapy.utils.RelativityWarning
+    : ~plasmapy.utils.RelativityWarning
         If the input velocity is greater than 5% of the speed of
         light.
 
@@ -1703,10 +1924,10 @@ def coupling_parameter(
     The degeneracy parameter is given by
 
     .. math::
-        \chi = n_e \Lambda_{deBroglie} ^ 3
+        \chi = n_e \Lambda_{de Broglie} ^ 3
 
-    where :math:`n_e` is the electron density and :math:`\Lambda_{deBroglie}`
-    is the thermal deBroglie wavelength.
+    where :math:`n_e` is the electron density and :math:`\Lambda_{de Broglie}`
+    is the thermal de Broglie wavelength.
 
     See equations 1.2, 1.3 and footnote 5 in [2]_ for details on the ideal
     chemical potential.
@@ -1728,11 +1949,10 @@ def coupling_parameter(
        approximation. D. O. Gericke et. al. PRE,  65, 036418 (2002).
        DOI: 10.1103/PhysRevE.65.036418
     .. [2] Bonitz, Michael. Quantum kinetic theory. Stuttgart: Teubner, 1998.
-
-
     """
     # boiler plate checks
     T, masses, charges, reduced_mass, V = _boilerPlate(T=T, species=species, V=V)
+
     if np.isnan(z_mean):
         # using mean charge to get average ion density.
         # If you are running this, you should strongly consider giving
@@ -1749,11 +1969,13 @@ def coupling_parameter(
         n_i = n_e / z_mean
         # getting Wigner-Seitz radius based on ion density
         radius = Wigner_Seitz_radius(n_i)
+
     # Coulomb potential energy between particles
     if np.isnan(z_mean):
         coulombEnergy = charges[0] * charges[1] / (4 * np.pi * eps0 * radius)
     else:
         coulombEnergy = (z_mean * e) ** 2 / (4 * np.pi * eps0 * radius)
+
     if method == "classical":
         # classical thermal kinetic energy
         kineticEnergy = k_B * T
diff --git a/plasmapy/formulary/ionization.py b/plasmapy/formulary/ionization.py
index 28b43535..ea795bec 100644
--- a/plasmapy/formulary/ionization.py
+++ b/plasmapy/formulary/ionization.py
@@ -102,13 +102,16 @@ def Saha(g_j, g_k, n_e: u.m ** -3, E_jk: u.J, T_e: u.K) -> u.dimensionless_unsca
     and recombination of ions with electrons are balanced.
 
     .. math::
-        \frac{N_j}{N_k} = \frac{1}{n_e} \frac{g_j}{4 g_k a_0^{3}} \left( \frac{k_B T_e}{\pi E_H} \right)^{\frac{3}{2}}
-         \exp\left( \frac{-E_{jk}}{k_B T_e} \right)
+        \frac{N_j}{N_k} = \frac{1}{n_e} \frac{g_j}{4 g_k a_0^{3}} \left(
+                              \frac{k_B T_e}{\pi E_H}
+                          \right)^{\frac{3}{2}}
+                          \exp\left( \frac{-E_{jk}}{k_B T_e} \right)
 
     Where :math:`k_B` is the Boltzmann constant,
     :math:`a_0` is the Bohr radius,
     :math:`E_H` is the ionization energy of Hydrogen,
-    :math:`N_j` and :math:`N_k` are the population of ions in the j and k states respectively.
+    :math:`N_j` and :math:`N_k` are the population of ions in the j and k states
+    respectively.
     This function is equivalent to Eq. 3.47 in `Drake`_.
 
     Parameters
@@ -164,9 +167,10 @@ def Saha(g_j, g_k, n_e: u.m ** -3, E_jk: u.J, T_e: u.K) -> u.dimensionless_unsca
 
     Notes
     -----
-    For reference to this function and for more information regarding the Saha equation,
-    see chapter 3 of R Paul Drake's book, "High-Energy-Density Physics: Foundation of
-    Inertial Fusion and Experimental Astrophysics" (`DOI: 10.1007/978-3-319-67711-8_3`_).
+    For reference to this function and for more information regarding the Saha
+    equation, see chapter 3 of R Paul Drake's book, "High-Energy-Density
+    Physics: Foundation of Inertial Fusion and Experimental Astrophysics" (
+    `DOI: 10.1007/978-3-319-67711-8_3`_).
 
     .. _`Drake`: https://doi.org/10.1007/978-3-319-67711-8
     .. _`DOI: 10.1007/978-3-319-67711-8_3`: https://doi.org/10.1007/978-3-319-67711-8_3
diff --git a/plasmapy/formulary/parameters.py b/plasmapy/formulary/parameters.py
index 53d09e9a..a408026e 100644
--- a/plasmapy/formulary/parameters.py
+++ b/plasmapy/formulary/parameters.py
@@ -906,8 +906,14 @@ def Hall_parameter(
         then the ``particle`` thermal velocity is assumed
         (`~plasmapy.formulary.parameters.thermal_speed`).
     coulomb_log_method : `str`, optional
-        Method used for Coulomb logarithm calculation. (see
-        `~plasmapy.formulary.collisions.Coulomb_logarithm`)
+        The method by which to compute the Coulomb logarithm.
+        The default method is the classical straight-line Landau-Spitzer
+        method (``"classical"`` or ``"ls"``). The other 6 supported methods
+        are ``"ls_min_interp"``, ``"ls_full_interp"``, ``"ls_clamp_mininterp"``,
+        ``"hls_min_interp"``, ``"hls_max_interp"``, and ``"hls_full_interp"``.
+        Please refer to the docstring of
+        `~plasmapy.formulary.collisions.Coulomb_logarithm` for more
+        information about these methods.
 
     See Also
     --------


Test patch:
diff --git a/plasmapy/dispersion/tests/test_two_fluid_dispersion.py b/plasmapy/dispersion/tests/test_two_fluid_dispersion.py
new file mode 100644
index 00000000..52113dc3
--- /dev/null
+++ b/plasmapy/dispersion/tests/test_two_fluid_dispersion.py
@@ -0,0 +1,212 @@
+"""Tests for the two fluid dispersion solution."""
+
+import numpy as np
+import pytest
+
+from astropy import units as u
+
+from plasmapy.dispersion.two_fluid_dispersion import (
+    tfds_,
+    two_fluid_dispersion_solution,
+)
+from plasmapy.formulary import parameters as pfp
+from plasmapy.particles import Particle
+from plasmapy.utils.exceptions import PhysicsWarning
+
+
+class TestTwoFluidDispersionSolution:
+    _kwargs_single_valued = {
+        "B": 8.3e-9 * u.T,
+        "ion": "p+",
+        "k": 0.0001 * u.rad / u.m,
+        "n_i": 5.0e6 * u.m ** -3,
+        "T_e": 1.6e6 * u.K,
+        "T_i": 4.0e5 * u.K,
+        "theta": 45 * u.deg,
+    }
+    _kwargs_bellan2012 = {
+        "B": 400e-4 * u.T,
+        "ion": Particle("He+"),
+        "n_i": 6.358e19 * u.m ** -3,
+        "T_e": 20 * u.eV,
+        "T_i": 10 * u.eV,
+        "k": (2 * np.pi * u.rad) / (0.56547 * u.m),
+    }
+
+    def test_alias(self):
+        """Test the associated alias is defined correctly."""
+        assert tfds_ is two_fluid_dispersion_solution
+
+    @pytest.mark.parametrize(
+        "kwargs, _error",
+        [
+            ({**_kwargs_single_valued, "B": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "B": [8e-9, 8.5e-9] * u.T}, ValueError),
+            ({**_kwargs_single_valued, "B": -1 * u.T}, ValueError),
+            ({**_kwargs_single_valued, "B": 5 * u.m}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "ion": {"not": "a particle"}}, TypeError),
+            ({**_kwargs_single_valued, "ion": "e-"}, ValueError),
+            ({**_kwargs_single_valued, "ion": "He", "z_mean": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "k": np.ones((3, 2)) * u.rad / u.m}, ValueError),
+            ({**_kwargs_single_valued, "k": 0 * u.rad / u.m}, ValueError),
+            ({**_kwargs_single_valued, "k": -1.0 * u.rad / u.m}, ValueError),
+            ({**_kwargs_single_valued, "k": 5 * u.s}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "n_i": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "n_i": [5e6, 6e6] * u.m ** -3}, ValueError),
+            ({**_kwargs_single_valued, "n_i": -5e6 * u.m ** -3}, ValueError),
+            ({**_kwargs_single_valued, "n_i": 2 * u.s}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "T_e": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "T_e": [1.4e6, 1.7e6] * u.K}, ValueError),
+            ({**_kwargs_single_valued, "T_e": -10 * u.eV}, ValueError),
+            ({**_kwargs_single_valued, "T_e": 2 * u.s}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "T_i": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "T_i": [4e5, 5e5] * u.K}, ValueError),
+            ({**_kwargs_single_valued, "T_i": -1 * u.eV}, ValueError),
+            ({**_kwargs_single_valued, "T_i": 2 * u.s}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "theta": np.ones((3, 2)) * u.deg}, ValueError),
+            ({**_kwargs_single_valued, "theta": 5 * u.eV}, u.UnitTypeError),
+            ({**_kwargs_single_valued, "gamma_e": "wrong type"}, TypeError),
+            ({**_kwargs_single_valued, "gamma_i": "wrong type"}, TypeError),
+        ],
+    )
+    def test_raises(self, kwargs, _error):
+        """Test scenarios that raise an `Exception`."""
+        with pytest.raises(_error):
+            two_fluid_dispersion_solution(**kwargs)
+
+    @pytest.mark.parametrize(
+        "kwargs, _warning",
+        [
+            # violates the low-frequency assumption (w/kc << 1)
+            (
+                {
+                    "B": 8.3e-7 * u.T,
+                    "ion": "p+",
+                    "k": 0.0001 * u.rad / u.m,
+                    "n_i": 3.0e6 * u.m ** -3,
+                    "T_e": 1.6e6 * u.K,
+                    "T_i": 4.0e5 * u.K,
+                    "theta": 5 * u.deg,
+                },
+                PhysicsWarning,
+            ),
+        ],
+    )
+    def test_warns(self, kwargs, _warning):
+        """Test scenarios the issue a `Warning`."""
+        with pytest.warns(_warning):
+            two_fluid_dispersion_solution(**kwargs)
+
+    @pytest.mark.parametrize(
+        "kwargs, expected",
+        [
+            (
+                {**_kwargs_bellan2012, "theta": 0 * u.deg},
+                {
+                    "fast_mode": 1.8631944,
+                    "alfven_mode": 0.5366538,
+                    "acoustic_mode": 0.4000832,
+                },
+            ),
+            (
+                {**_kwargs_bellan2012, "theta": 90 * u.deg},
+                {"fast_mode": 1.4000284, "alfven_mode": 0.0, "acoustic_mode": 0.0},
+            ),
+        ],
+    )
+    def test_on_bellan2012_vals(self, kwargs, expected):
+        """
+        Test calculated values based on Figure 1 of Bellan 2012
+        (DOI: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JA017856).
+        """
+        # theta and k values need to be single valued for this test to function
+        # correctly
+
+        cs = pfp.cs_(kwargs["T_e"], kwargs["T_i"], kwargs["ion"])
+        va = pfp.va_(kwargs["B"], kwargs["n_i"], ion=kwargs["ion"])
+        wci = pfp.wc_(kwargs["B"], kwargs["ion"])
+
+        beta = (cs / va).value ** 2
+        if not np.isclose(beta, 0.4, atol=1e-4):
+            pytest.fail(
+                f"The Bellan 2012 paper requires a 'beta' value of 0.4 and the test "
+                f"parameters yielded {beta:.6f}."
+            )
+
+        Lambda = (kwargs["k"] * va / wci).value ** 2
+        if not np.isclose(Lambda, 0.4, atol=1e-4):
+            pytest.fail(
+                f"The Bellan 2012 paper requires a 'Lambda' value of 0.4 and the test "
+                f"parameters yielded {Lambda:.6f}."
+            )
+
+        ws = two_fluid_dispersion_solution(**kwargs)
+        for mode, val in ws.items():
+            norm = (np.absolute(val) / (kwargs["k"] * va)).value ** 2
+            assert np.isclose(norm, expected[mode])
+
+    @pytest.mark.parametrize(
+        "kwargs, expected",
+        [
+            (
+                {
+                    **_kwargs_bellan2012,
+                    "ion": Particle("He"),
+                    "z_mean": 2.0,
+                    "theta": 0 * u.deg,
+                },
+                {**_kwargs_bellan2012, "ion": Particle("He +2"), "theta": 0 * u.deg},
+            ),
+            #
+            # z_mean defaults to 1
+            (
+                {**_kwargs_bellan2012, "ion": Particle("He"), "theta": 0 * u.deg},
+                {**_kwargs_bellan2012, "ion": Particle("He+"), "theta": 0 * u.deg},
+            ),
+        ],
+    )
+    def test_z_mean_override(self, kwargs, expected):
+        """Test overriding behavior of kw 'z_mean'."""
+        ws = two_fluid_dispersion_solution(**kwargs)
+        ws_expected = two_fluid_dispersion_solution(**expected)
+
+        for mode in ws:
+            assert np.isclose(ws[mode], ws_expected[mode], atol=0, rtol=1.7e-4)
+
+    @pytest.mark.parametrize(
+        "kwargs, expected",
+        [
+            ({**_kwargs_bellan2012, "theta": 0 * u.deg}, {"shape": ()}),
+            (
+                {
+                    **_kwargs_bellan2012,
+                    "theta": 0 * u.deg,
+                    "k": [1, 2, 3] * u.rad / u.m,
+                },
+                {"shape": (3,)},
+            ),
+            (
+                {
+                    **_kwargs_bellan2012,
+                    "theta": [10, 20, 30, 40, 50] * u.deg,
+                    "k": [1, 2, 3] * u.rad / u.m,
+                },
+                {"shape": (3, 5)},
+            ),
+            (
+                {**_kwargs_bellan2012, "theta": [10, 20, 30, 40, 50] * u.deg},
+                {"shape": (5,)},
+            ),
+        ],
+    )
+    def test_return_structure(self, kwargs, expected):
+        """Test the structure of the returned values."""
+        ws = two_fluid_dispersion_solution(**kwargs)
+
+        assert isinstance(ws, dict)
+        assert len({"acoustic_mode", "alfven_mode", "fast_mode"} - set(ws.keys())) == 0
+
+        for mode, val in ws.items():
+            assert isinstance(val, u.Quantity)
+            assert val.unit == u.rad / u.s
+            assert val.shape == expected["shape"]
diff --git a/plasmapy/formulary/tests/test_collisions.py b/plasmapy/formulary/tests/test_collisions.py
index 85495552..da42c46f 100644
--- a/plasmapy/formulary/tests/test_collisions.py
+++ b/plasmapy/formulary/tests/test_collisions.py
@@ -37,18 +37,31 @@ def setup_class(self):
         self.density2 = 1e23 * u.cm ** -3
         self.z_mean = 2.5
         self.particles = ("e", "p")
+        self.ls_min_interp = 3.4014290066940966
         self.gms1 = 3.4014290066940966
+        self.ls_min_interp_negative = -3.4310536971592493
         self.gms1_negative = -3.4310536971592493
+        self.ls_full_interp = 3.6349941014645157
         self.gms2 = 3.6349941014645157
+        self.ls_full_interp_negative = -1.379394033464292
         self.gms2_negative = -1.379394033464292
+        self.ls_clamp_mininterp = 3.4014290066940966
         self.gms3 = 3.4014290066940966
+        self.ls_clamp_mininterp_negative = 2
         self.gms3_negative = 2
+        self.ls_clamp_mininterp_non_scalar = (2, 2)
         self.gms3_non_scalar = (2, 2)
+        self.hls_min_interp = 3.401983996820073
         self.gms4 = 3.401983996820073
+        self.hls_min_interp_negative = 0.0005230791851781715
         self.gms4_negative = 0.0005230791851781715
+        self.hls_max_interp = 3.7196690506837693
         self.gms5 = 3.7196690506837693
+        self.hls_max_interp_negative = 0.03126832674323108
         self.gms5_negative = 0.03126832674323108
+        self.hls_full_interp = 3.635342040477818
         self.gms6 = 3.635342040477818
+        self.hls_full_interp_negative = 0.030720859361047514
         self.gms6_negative = 0.030720859361047514
 
     @pytest.mark.parametrize("insert_some_nans", [[], ["V"]])
@@ -57,11 +70,18 @@ def setup_class(self):
         "kwargs",
         [
             {"method": "classical"},
+            {"method": "ls"},
+            {"method": "ls_min_interp"},
             {"method": "GMS-1"},
+            {"method": "ls_full_interp", "z_mean": 1.0},
             {"method": "GMS-2", "z_mean": 1.0},
+            {"method": "ls_clamp_mininterp"},
             {"method": "GMS-3"},
+            {"method": "hls_min_interp"},
             {"method": "GMS-4"},
+            {"method": "hls_max_interp", "z_mean": 1.0},
             {"method": "GMS-5", "z_mean": 1.0},
+            {"method": "hls_full_interp", "z_mean": 1.0},
             {"method": "GMS-6", "z_mean": 1.0},
         ],
     )
@@ -262,6 +282,27 @@ def test_Chen_laser(self):
         )
         assert testTrue, errStr
 
+    def test_ls_min_interp(self):
+        """
+        Test for first version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=np.nan * u.dimensionless_unscaled,
+                V=np.nan * u.m / u.s,
+                method="ls_min_interp",
+            )
+        testTrue = np.isclose(methodVal, self.ls_min_interp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for ls_min_interp should be "
+            f"{self.ls_min_interp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS1(self):
         """
         Test for first version of Coulomb logarithm from Gericke,
@@ -283,13 +324,29 @@ def test_GMS1(self):
         )
         assert testTrue, errStr
 
+    def test_ls_min_interp_negative(self):
+        """
+        Test for first version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks for when
+        a negative (invalid) Coulomb logarithm is returned.
+        """
+        with pytest.warns(exceptions.CouplingWarning, match="depends on weak coupling"):
+            Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=np.nan * u.dimensionless_unscaled,
+                V=np.nan * u.m / u.s,
+                method="ls_min_interp",
+            )
+
     def test_GMS1_negative(self):
         """
         Test for first version of Coulomb logarithm from Gericke,
         Murillo, and Schlanges PRE (2002). This checks for when
         a negative (invalid) Coulomb logarithm is returned.
         """
-        with pytest.warns(exceptions.CouplingWarning, match="relies on weak coupling"):
+        with pytest.warns(exceptions.CouplingWarning, match="depends on weak coupling"):
             Coulomb_logarithm(
                 self.temperature2,
                 self.density2,
@@ -299,6 +356,27 @@ def test_GMS1_negative(self):
                 method="GMS-1",
             )
 
+    def test_ls_full_interp(self):
+        """
+        Test for second version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="ls_full_interp",
+            )
+        testTrue = np.isclose(methodVal, self.ls_full_interp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for GMS-2 should be "
+            f"{self.ls_full_interp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS2(self):
         """
         Test for second version of Coulomb logarithm from Gericke,
@@ -320,13 +398,29 @@ def test_GMS2(self):
         )
         assert testTrue, errStr
 
+    def test_ls_full_interp_negative(self):
+        """
+        Test for second version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks for when
+        a negative (invalid) Coulomb logarithm is returned.
+        """
+        with pytest.warns(exceptions.CouplingWarning, match="depends on weak coupling"):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="ls_full_interp",
+            )
+
     def test_GMS2_negative(self):
         """
         Test for second version of Coulomb logarithm from Gericke,
         Murillo, and Schlanges PRE (2002). This checks for when
         a negative (invalid) Coulomb logarithm is returned.
         """
-        with pytest.warns(exceptions.CouplingWarning, match="relies on weak coupling"):
+        with pytest.warns(exceptions.CouplingWarning, match="depends on weak coupling"):
             methodVal = Coulomb_logarithm(
                 self.temperature2,
                 self.density2,
@@ -336,6 +430,27 @@ def test_GMS2_negative(self):
                 method="GMS-2",
             )
 
+    def test_ls_clamp_mininterp(self):
+        """
+        Test for third version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="ls_clamp_mininterp",
+            )
+        testTrue = np.isclose(methodVal, self.ls_clamp_mininterp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for ls_clamp_minterp should be "
+            f"{self.ls_clamp_mininterp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS3(self):
         """
         Test for third version of Coulomb logarithm from Gericke,
@@ -357,6 +472,31 @@ def test_GMS3(self):
         )
         assert testTrue, errStr
 
+    def test_ls_clamp_mininterp_negative(self):
+        """
+        Test for third version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks whether
+        a positive value is returned whereas the classical Coulomb
+        logarithm would return a negative value.
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="ls_clamp_mininterp",
+            )
+        testTrue = np.isclose(
+            methodVal, self.ls_clamp_mininterp_negative, rtol=1e-6, atol=0.0
+        )
+        errStr = (
+            f"Coulomb logarithm for GMS-3 should be "
+            f"{self.ls_clamp_mininterp_negative} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS3_negative(self):
         """
         Test for third version of Coulomb logarithm from Gericke,
@@ -380,6 +520,31 @@ def test_GMS3_negative(self):
         )
         assert testTrue, errStr
 
+    def test_ls_clamp_mininterp_non_scalar_density(self):
+        """
+        Test for third version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks whether
+        passing in a collection of density values returns a
+        collection of Coulomb logarithm values.
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                10 * 1160 * u.K,
+                (1e23 * u.cm ** -3, 1e20 * u.cm ** -3),
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="ls_clamp_mininterp",
+            )
+        testTrue = np.isclose(
+            methodVal, self.ls_clamp_mininterp_non_scalar, rtol=1e-6, atol=0.0
+        )
+        errStr = (
+            f"Coulomb logarithm for GMS-3 should be "
+            f"{self.ls_clamp_mininterp_non_scalar} and not {methodVal}."
+        )
+        assert testTrue.all(), errStr
+
     def test_GMS3_non_scalar_density(self):
         """
         Test for third version of Coulomb logarithm from Gericke,
@@ -403,6 +568,27 @@ def test_GMS3_non_scalar_density(self):
         )
         assert testTrue.all(), errStr
 
+    def test_hls_min_interp(self):
+        """
+        Test for fourth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_min_interp",
+            )
+        testTrue = np.isclose(methodVal, self.hls_min_interp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for GMS-4 should be "
+            f"{self.hls_min_interp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS4(self):
         """
         Test for fourth version of Coulomb logarithm from Gericke,
@@ -424,6 +610,31 @@ def test_GMS4(self):
         )
         assert testTrue, errStr
 
+    def test_hls_min_interp_negative(self):
+        """
+        Test for fourth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks whether
+        a positive value is returned whereas the classical Coulomb
+        logarithm would return a negative value.
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_min_interp",
+            )
+        testTrue = np.isclose(
+            methodVal, self.hls_min_interp_negative, rtol=1e-6, atol=0.0
+        )
+        errStr = (
+            f"Coulomb logarithm for GMS-4 should be "
+            f"{self.hls_min_interp_negative} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS4_negative(self):
         """
         Test for fourth version of Coulomb logarithm from Gericke,
@@ -447,6 +658,27 @@ def test_GMS4_negative(self):
         )
         assert testTrue, errStr
 
+    def test_hls_max_interp(self):
+        """
+        Test for fifth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_max_interp",
+            )
+        testTrue = np.isclose(methodVal, self.hls_max_interp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for GMS-5 should be "
+            f"{self.hls_max_interp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS5(self):
         """
         Test for fifth version of Coulomb logarithm from Gericke,
@@ -468,6 +700,31 @@ def test_GMS5(self):
         )
         assert testTrue, errStr
 
+    def test_hls_max_interp_negative(self):
+        """
+        Test for fifth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks whether
+        a positive value is returned whereas the classical Coulomb
+        logarithm would return a negative value.
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_max_interp",
+            )
+        testTrue = np.isclose(
+            methodVal, self.hls_max_interp_negative, rtol=1e-6, atol=0.0
+        )
+        errStr = (
+            f"Coulomb logarithm for GMS-5 should be "
+            f"{self.hls_max_interp_negative} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS5_negative(self):
         """
         Test for fifth version of Coulomb logarithm from Gericke,
@@ -491,6 +748,27 @@ def test_GMS5_negative(self):
         )
         assert testTrue, errStr
 
+    def test_hls_full_interp(self):
+        """
+        Test for sixth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002).
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature1,
+                self.density1,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_full_interp",
+            )
+        testTrue = np.isclose(methodVal, self.hls_full_interp, rtol=1e-6, atol=0.0)
+        errStr = (
+            f"Coulomb logarithm for GMS-6 should be "
+            f"{self.hls_full_interp} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS6(self):
         """
         Test for sixth version of Coulomb logarithm from Gericke,
@@ -512,6 +790,31 @@ def test_GMS6(self):
         )
         assert testTrue, errStr
 
+    def test_hls_full_interp_negative(self):
+        """
+        Test for sixth version of Coulomb logarithm from Gericke,
+        Murillo, and Schlanges PRE (2002). This checks whether
+        a positive value is returned whereas the classical Coulomb
+        logarithm would return a negative value.
+        """
+        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                z_mean=self.z_mean,
+                V=np.nan * u.m / u.s,
+                method="hls_full_interp",
+            )
+        testTrue = np.isclose(
+            methodVal, self.hls_full_interp_negative, rtol=1e-6, atol=0.0
+        )
+        errStr = (
+            f"Coulomb logarithm for GMS-6 should be "
+            f"{self.hls_full_interp_negative} and not {methodVal}."
+        )
+        assert testTrue, errStr
+
     def test_GMS6_negative(self):
         """
         Test for sixth version of Coulomb logarithm from Gericke,
@@ -535,6 +838,19 @@ def test_GMS6_negative(self):
         )
         assert testTrue, errStr
 
+    def test_ls_full_interp_zmean_error(self):
+        """
+        Tests whether ls_full_interp raises z_mean error when a z_mean is not
+        provided.
+        """
+        with pytest.raises(ValueError):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                method="ls_full_interp",
+            )
+
     def test_GMS2_zmean_error(self):
         """
         Tests whether GMS-2 raises z_mean error when a z_mean is not
@@ -545,6 +861,19 @@ def test_GMS2_zmean_error(self):
                 self.temperature2, self.density2, self.particles, method="GMS-2"
             )
 
+    def test_hls_max_interp_zmean_error(self):
+        """
+        Tests whether hls_max_interp raises z_mean error when a z_mean is not
+        provided.
+        """
+        with pytest.raises(ValueError):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                method="hls_max_interp",
+            )
+
     def test_GMS5_zmean_error(self):
         """
         Tests whether GMS-5 raises z_mean error when a z_mean is not
@@ -555,6 +884,19 @@ def test_GMS5_zmean_error(self):
                 self.temperature2, self.density2, self.particles, method="GMS-5"
             )
 
+    def test_hls_full_interp_zmean_error(self):
+        """
+        Tests whether hls_full_interp raises z_mean error when a z_mean is not
+        provided.
+        """
+        with pytest.raises(ValueError):
+            methodVal = Coulomb_logarithm(
+                self.temperature2,
+                self.density2,
+                self.particles,
+                method="hls_full_interp",
+            )
+
     def test_GMS6_zmean_error(self):
         """
         Tests whether GMS-6 raises z_mean error when a z_mean is not