Repository: PlasmaPy/PlasmaPy

Issue:
Describe `plasma_calculator` in narrative documentation
Following up on #1229, one of the next things that we have to do for the `plasma_calculator` is to describe it somewhere in the narrative documentation.  A possibility would be to create a new page to describe how to access and use `plasma_calculator`, and then link to it from the main index page of the docs.  I'd probably add it to the "Package features" section, or another possibility would be the "All the Rest" section.

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) Report
> Merging [#1390](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (f9d028a) into [main](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/6df9583cc47375687a07300c0aa11ba31634d770?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (6df9583) will **increase** coverage by `0.00%`.
> The diff coverage is `n/a`.

> :exclamation: Current head f9d028a differs from pull request most recent head 308def6. Consider uploading reports for the commit 308def6 to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390/graphs/tree.svg?width=650&height=150&src=pr&token=LzZ5mnXVad&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)

```diff
@@           Coverage Diff           @@
##             main    #1390   +/-   ##
=======================================
  Coverage   97.11%   97.11%           
=======================================
  Files          76       76           
  Lines        7279     7289   +10     
=======================================
+ Hits         7069     7079   +10     
  Misses        210      210           
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/parameters.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3BhcmFtZXRlcnMucHk=) | `99.18% <0.00%> (+0.03%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Last update [6df9583...308def6](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1390?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).

RAJAGOPALAN-GANGADHARAN: pre-commit.ci autofix
RAJAGOPALAN-GANGADHARAN: @namurphy requesting for review, please let me know if it needs more details or any corrections
namurphy: Oh, and the Python 3.7 test failure appears completely unrelated.  I'm actually not sure why that test is failing but will look into it elsewhere.

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index eeca7a98..6341f79a 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -7,8 +7,8 @@ repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
   rev: v3.4.0
   hooks:
-  - id: check-yaml
   - id: check-json
+  - id: check-yaml
   - id: check-symlinks
   - id: end-of-file-fixer
   - id: trailing-whitespace
@@ -52,5 +52,3 @@ repos:
     args: [--autofix]
   - id: pretty-format-toml
     args: [--autofix]
-  - id: pretty-format-yaml
-    args: [--autofix, --indent, '2']
diff --git a/README.md b/README.md
index 67c843f2..89f57f46 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ # PlasmaPy
 [![YouTube](https://img.shields.io/badge/YouTube%20-subscribe-red?style=flat&logo=youtube)](https://www.youtube.com/channel/UCSH6qzslhqIZKTAJmHPxIxw)
 
 [![GitHub Actions — CI](https://github.com/PlasmaPy/PlasmaPy/workflows/CI/badge.svg)](https://github.com/PlasmaPy/PlasmaPy/actions?query=workflow%3ACI+branch%3Amain)
-[![weekly tests](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/weekly.yml/badge.svg?branch=main)](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/weekly.yml)
+[![fortnightly tests](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/fortnightly.yml/badge.svg?branch=main)](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/fortnightly.yml)
 [![Style linters](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/linters.yml/badge.svg)](https://github.com/PlasmaPy/PlasmaPy/actions/workflows/linters.yml)
 [![pre-commit.ci status](https://results.pre-commit.ci/badge/github/PlasmaPy/PlasmaPy/main.svg)](https://results.pre-commit.ci/latest/github/PlasmaPy/PlasmaPy/main)
 [![codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/branch/main/graph/badge.svg)](https://codecov.io/gh/PlasmaPy/PlasmaPy)
diff --git a/changelog/1189.feature.rst b/changelog/1189.feature.rst
new file mode 100644
index 00000000..f9210631
--- /dev/null
+++ b/changelog/1189.feature.rst
@@ -0,0 +1,3 @@
+Addition of the :file:`hollweg_.py` module to the `~plasmapy.dispersion.numerical`
+subpackage to numerically solve the dispersion relation using Hollweg's method,
+:cite:p:`hollweg:1999,bellan:2012`.
diff --git a/changelog/1292.trivial.rst b/changelog/1292.trivial.rst
new file mode 100644
index 00000000..1de67377
--- /dev/null
+++ b/changelog/1292.trivial.rst
@@ -0,0 +1,1 @@
+Added Python 3.10 to the GitHub Actions test suite.
diff --git a/changelog/1390.doc.rst b/changelog/1390.doc.rst
new file mode 100644
index 00000000..4ed69a82
--- /dev/null
+++ b/changelog/1390.doc.rst
@@ -0,0 +1,1 @@
+Described the |plasma-calculator| in the narrative documentation.
diff --git a/docs/api_static/plasmapy.dispersion.numerical.hollweg_.rst b/docs/api_static/plasmapy.dispersion.numerical.hollweg_.rst
new file mode 100644
index 00000000..970edfe9
--- /dev/null
+++ b/docs/api_static/plasmapy.dispersion.numerical.hollweg_.rst
@@ -0,0 +1,8 @@
+:orphan:
+
+`plasmapy.dispersion.numerical.hollweg_`
+========================================
+
+.. currentmodule:: plasmapy.dispersion.numerical.hollweg_
+
+.. automodapi:: plasmapy.dispersion.numerical.hollweg_
diff --git a/docs/api_static/plasmapy.dispersion.numerical.rst b/docs/api_static/plasmapy.dispersion.numerical.rst
new file mode 100644
index 00000000..79af74be
--- /dev/null
+++ b/docs/api_static/plasmapy.dispersion.numerical.rst
@@ -0,0 +1,8 @@
+:orphan:
+
+`plasmapy.dispersion.numerical`
+===============================
+
+.. currentmodule:: plasmapy.dispersion.numerical
+
+.. automodapi:: plasmapy.dispersion.numerical
diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 2fdaff23..a3c2a111 100644
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -148,6 +148,15 @@ @article{ji:2013
    number = {042114},
    doi = {10.1063/1.4801022}
 }
+@article{hollweg:1999,
+   author = {Hollweg, Joseph. V.},
+   title = {Kinetic Alfvén wave revisited},
+   journal = {Journal of Geophysical Research},
+   volume = {104},
+   number = {A7},
+   doi = {10.1029/1998JA900132},
+   year = {1999}
+}
 @book{khorikov:2020,
    author = {V. Khorikov},
    title = "{Unit Testing Principles, Practices, and Patterns}",
diff --git a/docs/common_links.rst b/docs/common_links.rst
index cc1ef450..0c587b51 100644
--- a/docs/common_links.rst
+++ b/docs/common_links.rst
@@ -169,6 +169,8 @@
 .. _nbsphinx: https://nbsphinx.readthedocs.io
 .. |nbsphinx| replace:: `nbsphinx`
 
+.. |plasma-calculator| replace:: :ref:`plasmapy-calculator`
+
 .. _`setup.cfg`: https://github.com/PlasmaPy/PlasmaPy/blob/main/setup.cfg
 .. |setup.cfg| replace:: :file:`setup.cfg`
 
diff --git a/docs/index.rst b/docs/index.rst
index 5720b945..79fe1af1 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -49,6 +49,7 @@ Example highlights
    Formulary <formulary/index>
    Particles <particles/index>
    Simulation <simulation/index>
+   Plasma Calculator <plasma_calculator/index>
    Plasma Objects <plasma/index>
    Package Utilities <utils/index>
 
diff --git a/docs/plasma_calculator/index.rst b/docs/plasma_calculator/index.rst
new file mode 100644
index 00000000..8140474e
--- /dev/null
+++ b/docs/plasma_calculator/index.rst
@@ -0,0 +1,26 @@
+.. _plasmapy-calculator:
+
+=================
+Plasma Calculator
+=================
+
+.. currentmodule:: plasmapy.utils.calculator
+
+
+Overview
+--------
+
+Plasma calculator is an interactive notebook installed alongside PlasmaPy that
+allows users to input a set of plasma properties and get immediate output for
+various functions at once.
+
+Using Plasma Calculator
+-----------------------
+
+To invoke the app use ``plasma-calculator`` in the command line.
+By default this opens the app in a browser, bootstrapping
+the notebook in a light theme and using port 8866 by default.
+
+``plasma-calculator`` takes optional arguments such as ``--dark``, ``--port``,
+and ``--no-browser``. Pass flag ``-h`` or ``--help`` to get full list of supported
+arguments.
diff --git a/plasmapy/dispersion/__init__.py b/plasmapy/dispersion/__init__.py
index 40fa63c9..58f2aaa2 100644
--- a/plasmapy/dispersion/__init__.py
+++ b/plasmapy/dispersion/__init__.py
@@ -7,9 +7,8 @@
     "plasma_dispersion_func",
     "plasma_dispersion_func_deriv",
 ]
-# TODO: import plasmapy.dispersion.numerical once it is populated with
-#       functionality (i.e. expose to namespace)
-from plasmapy.dispersion import analytical
+
+from plasmapy.dispersion import analytical, numerical
 from plasmapy.dispersion.dispersionfunction import (
     plasma_dispersion_func,
     plasma_dispersion_func_deriv,
diff --git a/plasmapy/dispersion/analytical/tests/__init__.py b/plasmapy/dispersion/analytical/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/plasmapy/dispersion/numerical/__init__.py b/plasmapy/dispersion/numerical/__init__.py
index 3a0b405e..974e87c2 100644
--- a/plasmapy/dispersion/numerical/__init__.py
+++ b/plasmapy/dispersion/numerical/__init__.py
@@ -2,4 +2,6 @@
 The `~plasmapy.dispersion.numerical` subpackage contains functionality
 associated with numerical dispersion solvers.
 """
-__all__ = []
+__all__ = ["hollweg"]
+
+from plasmapy.dispersion.numerical.hollweg_ import hollweg
diff --git a/plasmapy/dispersion/numerical/hollweg_.py b/plasmapy/dispersion/numerical/hollweg_.py
new file mode 100644
index 00000000..c5aa07ea
--- /dev/null
+++ b/plasmapy/dispersion/numerical/hollweg_.py
@@ -0,0 +1,360 @@
+"""
+This module contains functionality for calculating various numerical
+solutions to Hollweg's two fluid dispersion relation
+"""
+
+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
+
+c_si_unitless = c.value
+
+
+@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 hollweg(
+    *,
+    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.rad,
+    gamma_e: Union[float, int] = 1,
+    gamma_i: Union[float, int] = 3,
+    z_mean: Union[float, int] = None,
+):
+    r"""
+    Calculate the two fluid dispersion relation presented by
+    :cite:t:`hollweg:1999`, and discussed by :cite:t:`bellan:2012`.
+    This is a numberical solver of equation 3 in :cite:t:`bellan:2012`.
+    See the **Notes** section below for additional details.
+
+    Parameters
+    ----------
+    B : `~astropy.units.Quantity`
+        The magnetic field magnitude in units convertible to 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 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 m\ :sup:`-3`.
+    T_e : `~astropy.units.Quantity`
+        The electron temperature in units of K or eV.
+    T_i : `~astropy.units.Quantity`
+        The ion temperature in units of K or 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 convertible
+        to radians.  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 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 :math:`\omega / \omega_{\rm ci} > 0.1`, violation of the
+        low-frequency assumption.
+
+    : `~plasmapy.utils.exceptions.PhysicsWarning`
+        When :math:`c_{\rm s} / v_{\rm A} > 0.1`, violation of low-β.
+
+    : `~plasmapy.utils.exceptions.PhysicsWarning`
+        When :math:`|θ - π/2| > 0.1`, violation of quasi-perpendicular
+        propagation.
+
+    Notes
+    -----
+
+    The dispersion relation presented in :cite:t:`hollweg:1999`
+    (equation 3 in :cite:t:`bellan:2012`) is:
+
+    .. math::
+        \left( \frac{\omega^2}{k_{\rm z}^2 v_{\rm A}^2} - 1 \right) &
+        \left[
+            \omega^2 \left( \omega^2 - k^2 v_{\rm A}^2 \right)
+            - \beta k^2 v_{\rm A}^2 \left(
+                \omega^2 - k_{\rm z}^2 v_{\rm A}^2
+            \right)
+        \right] \\
+        &= \omega^2 \left(\omega^2 - k^2 v_{\rm A}^2 \right) k_{\rm x}^2
+        \left(
+            \frac{c_{\rm s}^2}{\omega_{\rm ci}^2}
+            - \frac{c^2}{\omega_{\rm pe}^2} \frac{\omega^2}{k_{\rm z}^2v_{\rm A}^2}
+        \right)
+
+    where
+
+    .. math::
+        \mathbf{B_o} &= B_{o} \mathbf{\hat{z}} \\
+        \cos \theta &= \frac{k_z}{k} \\
+        \mathbf{k} &= k_{\rm x} \hat{x} + k_{\rm z} \hat{z}
+
+    :math:`\omega` is the wave frequency, :math:`k` is the wavenumber,
+    :math:`v_{\rm A}` is the Alfvén velocity, :math:`c_{\rm s}` is the
+    sound speed, :math:`\omega_{\rm ci}` is the ion gyrofrequency, and
+    :math:`\omega_{\rm pe}` is the electron plasma frequency. In the
+    derivation of this relation Hollweg assumed low-frequency waves
+    :math:`\omega / \omega_{\rm ci} \ll 1`, no D.C. electric field
+    :math:`\mathbf{E_o}=0`, and quasi-neutrality.
+
+    :cite:t:`hollweg:1999` asserts this expression is valid for
+    arbitrary :math:`c_{\rm s} / v_{\rm A}` (β) and
+    :math:`k_{\rm z} / k` (θ).  Contrarily, :cite:t:`bellan:2012`
+    states in §1.7 that due to the inconsistent retention of the
+    :math:`\omega / \omega_{\rm ci} \ll 1` terms the expression can
+    only be valid if both :math:`c_{\rm s} \ll v_{\rm A}` (low-β) and
+    the wave propgation is nearly perpendicular to the magnetic field.
+
+    This routine solves for ω for given :math:`k` values by numerically
+    solving for the roots of the above expression.
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from plasmapy.dispersion.numerical import hollweg_
+    >>> inputs = {
+    ...    "k": np.logspace(-7, -2, 2) * u.rad / u.m,
+    ...    "theta": 88 * u.deg,
+    ...    "n_i": 5 * u.cm ** -3,
+    ...    "B": 2.2e-8 * u.T,
+    ...    "T_e": 1.6e6 * u.K,
+    ...    "T_i": 4.0e5 * u.K,
+    ...    "ion": Particle("p+"),
+    ... }
+    >>> omegas = hollweg(**inputs)
+    >>> omegas
+    {'fast_mode': <Quantity [2.62911663e-02+0.j, 2.27876968e+03+0.j] rad / s>,
+     'alfven_mode': <Quantity [7.48765909e-04+0.j, 2.13800404e+03+0.j] rad / s>,
+     'acoustic_mode': <Quantity [0.00043295+0.j, 0.07358991+0.j] 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("The particle passed for 'ion' must be an ion or element.")
+
+    # validate z_mean
+    if z_mean is None:
+        try:
+            z_mean = abs(ion.charge_number)
+        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 be single valued 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 single valued or a 1D array astropy Quantity,"
+            f" got array of shape {k.shape}."
+        )
+    if np.any(k <= 0):
+        raise ValueError("Argument 'k' can not be a or have negative values.")
+
+    # validate argument theta
+    theta = theta.squeeze()
+    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}."
+        )
+
+    # Single k value case
+    if np.isscalar(k.value):
+        k = np.array([k.value]) * u.rad / u.m
+
+    # Calc needed plasma parameters
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", category=PhysicsWarning)
+        n_e = z_mean * n_i
+        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,
+        ).value
+        v_A = pfp.Alfven_speed(B, n_i, ion=ion, z_mean=z_mean).value
+        omega_ci = pfp.gyrofrequency(B=B, particle=ion, signed=False, Z=z_mean).value
+        omega_pe = pfp.plasma_frequency(n=n_e, particle="e-").value
+
+    # strip units from select input args
+    k = k.value
+    theta = theta.value
+    cs_vA = c_s / v_A
+
+    # Parameters kx and kz
+    kz = np.cos(theta) * k
+    kx = np.sin(theta) * k
+
+    # Define helpful parameters
+    beta = (c_s / v_A) ** 2
+    alpha_A = (k * v_A) ** 2
+    alpha_s = (k * c_s) ** 2  # == alpha_A * beta
+    sigma = (kz * v_A) ** 2
+    D = (c_s / omega_ci) ** 2
+    F = (c_si_unitless / omega_pe) ** 2
+
+    # Polynomial coefficients: c3*x^3 + c2*x^2 + c1*x + c0 = 0
+    c3 = F * kx ** 2 + 1
+    c2 = -alpha_A * (1 + beta + F * kx ** 2) - sigma * (1 + D * kx ** 2)
+    c1 = sigma * alpha_A * (1 + 2 * beta + D * kx ** 2)
+    c0 = -alpha_s * sigma ** 2
+
+    # Find roots to polynomial
+    coefficients = np.array([c3, c2, c1, c0], ndmin=2)
+    nroots = coefficients.shape[0] - 1  # 3
+    nks = coefficients.shape[1]
+    roots = np.empty((nroots, nks), dtype=np.complex128)
+    for ii in range(nks):
+        roots[:, ii] = np.roots(coefficients[:, ii])[:]
+
+    roots = np.sqrt(roots)
+    roots = np.sort(roots, axis=0)
+
+    # Warn about NOT low-beta
+    if c_s / v_A > 0.1:
+        warnings.warn(
+            f"This solver is valid in the low-beta regime, "
+            f"c_s/v_A << 1 according to Bellan, 2012, Sec. 1.7 "
+            f"(see documentation for DOI). A c_s/v_A value of {cs_vA:.2f} "
+            f"was calculated which may affect the validity of the solution.",
+            PhysicsWarning,
+        )
+
+    # Warn about theta not nearly perpendicular
+    if np.abs(theta - np.pi / 2) > 0.1:
+        warnings.warn(
+            f"This solver is valid in the regime where propogation is "
+            f"nearly perpendicular to B according to Bellan, 2012, Sec. 1.7 "
+            f"(see documentation for DOI). A theta value of {theta:.2f} was "
+            f"entered which may affect the validity of the solution.",
+            PhysicsWarning,
+        )
+
+    # dispersion relation is only valid in the regime w << w_ci
+    w_max = np.max(roots)
+    w_wci_max = w_max / omega_ci
+    if w_wci_max > 0.1:
+        warnings.warn(
+            f"This solver is valid in the regime w/w_ci << 1.  A w "
+            f"value of {w_max:.2f} and a w/w_ci value of "
+            f"{w_wci_max:.2f} were calculated which may affect the "
+            f"validity of the solution.",
+            PhysicsWarning,
+        )
+
+    omegas = {
+        "fast_mode": roots[2, :].squeeze() * u.rad / u.s,
+        "alfven_mode": roots[1, :].squeeze() * u.rad / u.s,
+        "acoustic_mode": roots[0, :].squeeze() * u.rad / u.s,
+    }
+
+    return omegas
diff --git a/plasmapy/dispersion/tests/__init__.py b/plasmapy/dispersion/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/setup.cfg b/setup.cfg
index 12c32dbb..ca51d660 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -28,6 +28,7 @@ classifiers=
   Programming Language :: Python :: 3.8
   Programming Language :: Python :: 3.8
   Programming Language :: Python :: 3.9
+  Programming Language :: Python :: 3.10
   Programming Language :: Python :: Implementation :: CPython
   Topic :: Scientific/Engineering :: Physics
   Topic :: Scientific/Engineering :: Astronomy


Test patch:
diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index 8163f43f..56f07500 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -20,10 +20,10 @@ jobs:
       matrix:
         include:
 
-        - name: Python 3.9
+        - name: Python 3.10
           os: ubuntu-latest
-          python: 3.9
-          toxenv: py39
+          python: "3.10"
+          toxenv: py310
 
         - name: Linters
           os: ubuntu-latest
@@ -78,6 +78,12 @@ jobs:
           python: 3.8
           toxenv: py38-all
 
+        - name: Python 3.9
+          os: ubuntu-latest
+          python: 3.9
+          toxenv: py39
+
+
     steps:
     - name: Checkout code
       uses: actions/checkout@v2.4.0
diff --git a/plasmapy/dispersion/numerical/tests/__init__.py b/plasmapy/dispersion/numerical/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/plasmapy/dispersion/numerical/tests/test_hollweg_.py b/plasmapy/dispersion/numerical/tests/test_hollweg_.py
new file mode 100644
index 00000000..f27f052d
--- /dev/null
+++ b/plasmapy/dispersion/numerical/tests/test_hollweg_.py
@@ -0,0 +1,271 @@
+"""Tests for the hollweg dispersion solution."""
+
+import numpy as np
+import pytest
+
+from astropy import units as u
+
+from plasmapy.dispersion.numerical.hollweg_ import hollweg
+from plasmapy.formulary import parameters as pfp
+from plasmapy.particles import Particle
+from plasmapy.utils.exceptions import PhysicsWarning
+
+
+class TestHollweg:
+    _kwargs_single_valued = {
+        "k": 0.01 * u.rad / u.m,
+        "theta": 88 * u.deg,
+        "n_i": 5 * u.cm ** -3,
+        "B": 2.2e-8 * u.T,
+        "T_e": 1.6e6 * u.K,
+        "T_i": 4.0e5 * u.K,
+        "ion": Particle("p+"),
+    }
+
+    _kwargs_hollweg1999 = {
+        "theta": 90 * u.deg,
+        "n_i": 5 * u.cm ** -3,
+        "T_e": 1.6e6 * u.K,
+        "T_i": 4.0e5 * u.K,
+        "ion": Particle("p+"),
+    }
+
+    @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):
+            hollweg(**kwargs)
+
+    @pytest.mark.parametrize(
+        "kwargs, _warning",
+        [
+            # w/w_ci << 1 PhysicsWarning
+            (
+                {
+                    "k": 0.01 * u.rad / u.m,
+                    "theta": 88 * u.deg,
+                    "n_i": 0.05 * u.cm ** -3,
+                    "B": 2.2e-8 * u.T,
+                    "T_e": 1.6e6 * u.K,
+                    "T_i": 4.0e5 * u.K,
+                    "ion": Particle("p+"),
+                },
+                PhysicsWarning,
+            ),
+            # c_s/v_A << 1 PhysicsWarning
+            (
+                {
+                    "k": 10e-8 * u.rad / u.m,
+                    "theta": 88 * u.deg,
+                    "n_i": 5 * u.cm ** -3,
+                    "B": 6.98e-8 * u.T,
+                    "T_e": 1.6e6 * u.K,
+                    "T_i": 4.0e5 * u.K,
+                    "ion": Particle("p+"),
+                },
+                PhysicsWarning,
+            ),
+            # theta nearly perpendicular PhysicsWarning
+            (
+                {
+                    "k": 10e-8 * u.rad / u.m,
+                    "theta": 84 * u.deg,
+                    "n_i": 1 * u.cm ** -3,
+                    "B": 6.98e-8 * u.T,
+                    "T_e": 1.6e6 * u.K,
+                    "T_i": 4.0e5 * u.K,
+                    "ion": Particle("p+"),
+                },
+                PhysicsWarning,
+            ),
+        ],
+    )
+    def test_warning(self, kwargs, _warning):
+        """Test scenarios that raise a `Warning`."""
+        with pytest.warns(_warning):
+            hollweg(**kwargs)
+
+    @pytest.mark.parametrize(
+        "kwargs, expected, desired_beta",
+        [
+            (  # beta = 1/20 for kx*L = 0
+                {**_kwargs_hollweg1999, "k": 1e-14 * u.rad / u.m, "B": 6.971e-8 * u.T},
+                1 + 0j,
+                1 / 20,
+            ),
+            (  # beta = 1/20 for kx*L = 1
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 0.0000439223874624874 * u.rad / u.m,
+                    "B": 6.971e-8 * u.T,
+                },
+                1.4018 + 0j,
+                1 / 20,
+            ),
+            (  # beta = 1/2 for kx*L = 0
+                {**_kwargs_hollweg1999, "k": 1e-14 * u.rad / u.m, "B": 2.205e-8 * u.T},
+                1 + 0j,
+                0.5,
+            ),
+            (  # beta = 1/2 for kx*L = 1
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 0.000013893109303101 * u.rad / u.m,
+                    "B": 2.205e-8 * u.T,
+                },
+                1.3536 + 0j,
+                0.5,
+            ),
+            (  # beta = 2 for kx*L = 0
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 1e-14 * u.rad / u.m,
+                    "B": 1.10232e-8 * u.T,
+                },
+                1 + 0j,
+                2,
+            ),
+            (  # beta = 2 for kx*L = 1
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 0.00000691190063354451 * u.rad / u.m,
+                    "B": 1.10232e-8 * u.T,
+                },
+                1.2607 + 0j,
+                2,
+            ),
+            (  # beta = 1/2000 for kx*L = 0
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 1e-14 * u.rad / u.m,
+                    "B": 6.97178e-7 * u.T,
+                },
+                1 + 0j,
+                1 / 2000,
+            ),
+            (  # beta = 1/2000 for kx*L = 1
+                {
+                    **_kwargs_hollweg1999,
+                    "k": 0.000439273010336778 * u.rad / u.m,
+                    "B": 6.97178e-7 * u.T,
+                },
+                0.98750 + 0j,
+                1 / 2000,
+            ),
+        ],
+    )
+    def test_hollweg1999_vals(self, kwargs, expected, desired_beta):
+        """
+        Test calculated values based on Figure 2 of Hollweg1999
+        (DOI: https://doi.org/10.1029/1998JA900132) using eqn 3 of
+        Bellan 2012
+        (DOI: https://doi.org/10.1029/2012JA017856).
+
+        The WebPlotDigitizer software was used to determine the test
+        parameters for k, B, and expected omega from Figure 2 of
+        Hollweg1999.
+
+        - GitHub: https://github.com/ankitrohatgi/WebPlotDigitizer
+        - Web Version: https://automeris.io/WebPlotDigitizer/
+        """
+        # k values need to be single valued for this test to function correctly
+
+        cs = pfp.cs_(kwargs["T_e"], kwargs["T_i"], kwargs["ion"]).value
+        va = pfp.va_(kwargs["B"], kwargs["n_i"], ion=kwargs["ion"]).value
+
+        beta = (cs / va) ** 2
+        if not np.isclose(beta, desired_beta, atol=2e-4):
+            pytest.fail(
+                f"The Holweg 1999 paper requires a 'beta' value of {desired_beta:0.5f} "
+                f"and the test parameters yielded {beta:.6f}."
+            )
+
+        kz = (np.cos(kwargs["theta"]) * kwargs["k"]).value
+
+        w_alfven = (hollweg(**kwargs)["alfven_mode"]).value
+        big_omega = np.abs(w_alfven / (kz * va))
+
+        assert np.allclose(big_omega, expected, atol=1e-2)
+
+    @pytest.mark.parametrize(
+        "kwargs, expected",
+        [
+            (
+                {
+                    **_kwargs_single_valued,
+                    "ion": Particle("He"),
+                    "z_mean": 2.0,
+                },
+                {**_kwargs_single_valued, "ion": Particle("He +2")},
+            ),
+            #
+            # z_mean defaults to 1
+            (
+                {**_kwargs_single_valued, "ion": Particle("He")},
+                {**_kwargs_single_valued, "ion": Particle("He+")},
+            ),
+        ],
+    )
+    def test_z_mean_override(self, kwargs, expected):
+        """Test overriding behavior of kw 'z_mean'."""
+        ws = hollweg(**kwargs)
+        ws_expected = hollweg(**expected)
+
+        for mode in ws:
+            assert np.isclose(ws[mode], ws_expected[mode], atol=1e-5, rtol=1.7e-4)
+
+    @pytest.mark.parametrize(
+        "kwargs, expected",
+        [
+            ({**_kwargs_single_valued}, {"shape": ()}),
+            (
+                {
+                    **_kwargs_single_valued,
+                    "k": [1, 2, 3] * u.rad / u.m,
+                },
+                {"shape": (3,)},
+            ),
+        ],
+    )
+    def test_return_structure(self, kwargs, expected):
+        """Test the structure of the returned values."""
+        ws = hollweg(**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"]