Repository: PlasmaPy/PlasmaPy
Pull request #911 (2020-10-05 20:41:47)

Starting an attempt to remove default particle assumptions
This is my attempt to remove default values from functions that take a particle as input in `formulary/parameters.py`.

This is described in issue #453.

- [x] I have added a changelog entry for this pull request.

<!--
* `breaking`: the `ion` argument in the `ion_sound_speed` function is now a keyword-only argument. I'm not sure how to remove the default value without making it keyword-only.

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. (Well, mostly I changed the docstrings to reflect the fact that some arguments that used to be optional now are required.

<!--
(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!)
-->


rocco8773: @arichar6 Thanks for taking this on!  Issue #453 was created before I started working on `plasmapy`, but I absolutely agree with it.  With that said,  we shouldn't even give the `ion` and `particle` keywords are default value of `None`.  Instead we should make them a required argument like...

```python
def Alfven_speed(B, density, ion, z_mean=None):
    ...
```

or


```python
def Alfven_speed(B, density, *, ion, z_mean=None):
    ...
```

The latter keeps `ion` as a keyword-only argument, whereas, the former makes it a keyword and positional argument.

Additionally, we should add a proper annotation `Union[str, Particle]`, which would look something like...

```python
def Alfven_speed(B: u.T, density: [u.m ** -3, u.kg / u.m ** 3], *, ion: Union[str, Particle], z_mean=None):
    ...
```

arichar6: Thanks for the feedback!

I'll remove the `None` defaults, and add the type annotations.
arichar6: I still have to fix the tests that now fail. But the defaults should be fixed now, and type annotations added.
arichar6: As I was working on the tests, I realized that one change might introduce a backwards-incompatibility. 
I made the `ion` argument in the `ion_sound_speed` function a keyword-only argument. I'm not sure how to remove the default value without making it keyword-only, since arguments that come before it have default values. 

An alternative would be to change the order of the arguments to the function, but that seems worse somehow.

Thoughts?
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/911?src=pr&el=h1) Report
> :exclamation: No coverage uploaded for pull request base (`master@6d5651c`). [Click here to learn what that means](https://docs.codecov.io/docs/error-reference#section-missing-base-commit).
> The diff coverage is `100.00%`.

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

```diff
@@            Coverage Diff            @@
##             master     #911   +/-   ##
=========================================
  Coverage          ?   96.21%           
=========================================
  Files             ?       60           
  Lines             ?     5418           
  Branches          ?        0           
=========================================
  Hits              ?     5213           
  Misses            ?      205           
  Partials          ?        0           
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/911?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/parameters.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/911/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3BhcmFtZXRlcnMucHk=) | `99.01% <100.00%> (ø)` | |

------

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

arichar6: One more question. Now that we are using `Particle` as the type annotation, the docstrings should probably be updated from saying something like `ion : str` to `ion : Particle`, correct?
namurphy: > One more question. Now that we are using `Particle` as the type annotation, the docstrings should probably be updated from saying something like `ion : str` to `ion : Particle`, correct?

Hm...for the moment, I'd suggest using `ion:  Particle`, but we're having a bit of discussion in #912 about how to best do this.  
rocco8773: @arichar6 It appears the documentation build errors are due to the jupyter notebooks not being updated, specifically `/docs/notebooks/physics.ipynb`.  Do you want to give fixing this a try?  Just make sure before committing your changes that you clear all cell outputs.
arichar6: Any insights into why the codecov test is failing now? I'm don't really quite know what it is check for, exactly.
rocco8773: > Any insights into why the codecov test is failing now? I'm don't really quite know what it is check for, exactly.

Looks like line 339 in `collisions.py` is not getting hit by tests.  Without digging more, I'm not sure why this is the case.
arichar6: > 
> Looks like line 339 in `collisions.py` is not getting hit by tests. Without digging more, I'm not sure why this is the case.

I might be wrong, but it looks like line 339 in `collisions.py` is not accessible. Here are the precedding lines:
```python
elif np.any(np.isnan(V)):
    if np.isscalar(V.value) and np.isscalar(T.value):
        V = parameters.thermal_speed(T, "e-", mass=m)
    elif np.isscalar(V.value):
        V = parameters.thermal_speed(T, "e-", mass=m)
```
I just tried this, and if `np.any(np.isnan(V))` returns `True` for a numpy array (with astropy units), then `np.isscalar(V.value)` returns `False`. Was there some other case that this was trying to test for?

namurphy: > I might be wrong, but it looks like line 339 in `collisions.py` is not accessible. Here are the preceding lines:

> I just tried this, and if `np.any(np.isnan(V))` returns `True` for a numpy array (with astropy units), then `np.isscalar(V.value)` returns `False`. Was there some other case that this was trying to test for?

Hm...that's a good question.  That's a bit strange so it would be worth looking into it further.  It does seem outside of the intended scope of this pull request, so my suggestion would be to create a brief issue on it so that we can address it in a different pull request.  I'm not a fan of long extensive `if`/`elif` statements so I'm wondering if there's a way to [replace these `if`/`elif` blocks with polymorphism](https://www.codementor.io/@uditrastogi/replace-conditional-statements-if-else-or-switch-with-polymorphism-ryl8mx4ns)...but that is definitely for another day. 
namurphy: To follow up: there are a few other things in `collisions.py` that should probably be cleaned up, so I can create the other issue.
arichar6: Ok. I think this PR is finally done. While this does not complete the changes required for #453, it does complete the changes needed in `parameters.py`.

Patch:
diff --git a/changelog/911.breaking.rst b/changelog/911.breaking.rst
new file mode 100644
index 00000000..253e806b
--- /dev/null
+++ b/changelog/911.breaking.rst
@@ -0,0 +1,1 @@
+Removed default values for the `ion` and `particle` arguments of functions contained in `plasmapy.formulary.parameters`, in accordance with issue [#453](https://github.com/PlasmaPy/PlasmaPy/issues/453), and updated all relevant calls to modified functionality.
diff --git a/docs/about/credits.rst b/docs/about/credits.rst
index 060b79da..abdc2d55 100644
--- a/docs/about/credits.rst
+++ b/docs/about/credits.rst
@@ -85,6 +85,7 @@ in parentheses are `ORCID author identifiers <https://orcid.org>`__.
 * `Sixue Xu <https://github.com/hzxusx>`__
 * `Carol Zhang <https://github.com/carolyz>`__
 * `Cody Skinner <https://github.com/wskinner74>`__
+* `Steve Richardson <https://github.com/arichar6>`__ (`0000-0002-3056-6334 <https://orcid.org/0000-0002-3056-6334>`__)
 
 This list contains contributors to PlasmaPy's core package and vision
 statement, including a few people who do not show up as `PlasmaPy
diff --git a/docs/notebooks/physics.ipynb b/docs/notebooks/physics.ipynb
index 2292d27a..ba9b6582 100644
--- a/docs/notebooks/physics.ipynb
+++ b/docs/notebooks/physics.ipynb
@@ -187,7 +187,7 @@
    },
    "outputs": [],
    "source": [
-    "print(formulary.plasma_frequency(electron_concentration))"
+    "print(formulary.plasma_frequency(electron_concentration, particle=\"e-\"))"
    ]
   },
   {
@@ -214,10 +214,10 @@
    "outputs": [],
    "source": [
     "n_e = np.logspace(4, 30, 100) / u.m**3\n",
-    "plt.plot(n_e, formulary.plasma_frequency(n_e))\n",
+    "plt.plot(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"))\n",
     "plt.scatter(\n",
     "    electron_concentration,\n",
-    "    formulary.plasma_frequency(electron_concentration),\n",
+    "    formulary.plasma_frequency(electron_concentration, particle=\"e-\"),\n",
     "    label=\"Our Data\")\n",
     "\n",
     "# IRT1 Tokamak Data\n",
@@ -225,14 +225,14 @@
     "n_e = 1.2e19 / u.m**3\n",
     "T_e = 136.8323 * u.eV\n",
     "B = 0.82 * u.T\n",
-    "plt.scatter(n_e, formulary.plasma_frequency(n_e), label=\"IRT1 Tokamak\")\n",
+    "plt.scatter(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"), label=\"IRT1 Tokamak\")\n",
     "\n",
     "# Wendelstein 7-X Stellerator Data\n",
     "# https://nucleus.iaea.org/sites/fusionportal/Shared%20Documents/FEC%202016/fec2016-preprints/preprint0541.pdf\n",
     "n_e = 3e19 / u.m**3\n",
     "T_e = 6 * u.keV\n",
     "B = 3 * u.T\n",
-    "plt.scatter(n_e, formulary.plasma_frequency(n_e), label=\"W7-X Stellerator\")\n",
+    "plt.scatter(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"), label=\"W7-X Stellerator\")\n",
     "\n",
     "# Solar Corona\n",
     "# Estimated by Nick Murphy\n",
@@ -240,7 +240,7 @@
     "T_e = 1e6 * u.K\n",
     "B = 0.005 * u.T\n",
     "T_e.to(u.eV, equivalencies=u.temperature_energy())\n",
-    "plt.scatter(n_e, formulary.plasma_frequency(n_e), label=\"Solar Corona\")\n",
+    "plt.scatter(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"), label=\"Solar Corona\")\n",
     "\n",
     "# Interstellar (warm neutral) Medium\n",
     "# Estimated by Nick Murphy\n",
@@ -248,7 +248,7 @@
     "T_e = 5e3 * u.K\n",
     "B = 0.005 * u.T\n",
     "T_e.to(u.eV, equivalencies=u.temperature_energy())\n",
-    "plt.scatter(n_e, formulary.plasma_frequency(n_e), label=\"Interstellar Medium\")\n",
+    "plt.scatter(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"), label=\"Interstellar Medium\")\n",
     "\n",
     "# Solar Wind at 1 AU\n",
     "# Estimated by Nick Murphy\n",
@@ -256,7 +256,7 @@
     "T_e = 1e5 * u.K\n",
     "B = 5e-9 * u.T\n",
     "T_e.to(u.eV, equivalencies=u.temperature_energy())\n",
-    "plt.scatter(n_e, formulary.plasma_frequency(n_e), label=\"Solar Wind (1AU)\")\n",
+    "plt.scatter(n_e, formulary.plasma_frequency(n_e, particle=\"e-\"), label=\"Solar Wind (1AU)\")\n",
     "\n",
     "\n",
     "plt.xlabel(\"Electron Concentration (m^-3)\")\n",
diff --git a/plasmapy/formulary/collisions.py b/plasmapy/formulary/collisions.py
index f4668d06..08d1a5c1 100644
--- a/plasmapy/formulary/collisions.py
+++ b/plasmapy/formulary/collisions.py
@@ -331,18 +331,16 @@ def _replaceNanVwithThermalV(V, T, m):
         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, mass=m)
+        V = parameters.thermal_speed(T, "e-", mass=m)
     elif np.any(np.isnan(V)):
-        if np.isscalar(V.value) and np.isscalar(T.value):
-            V = parameters.thermal_speed(T, mass=m)
-        elif np.isscalar(V.value):
-            V = parameters.thermal_speed(T, mass=m)
-        elif np.isscalar(T.value):
-            V = V.copy()
-            V[np.isnan(V)] = parameters.thermal_speed(T, mass=m)
+        if np.isscalar(V.value) or np.isscalar(T.value):
+            if np.isscalar(V.value):
+                V = parameters.thermal_speed(T, "e-", mass=m)
+            if np.isscalar(T.value):
+                V[np.isnan(V)] = parameters.thermal_speed(T, "e-", mass=m)
         else:
             V = V.copy()
-            V[np.isnan(V)] = parameters.thermal_speed(T[np.isnan(V)], mass=m)
+            V[np.isnan(V)] = parameters.thermal_speed(T[np.isnan(V)], "e-", mass=m)
     return V
 
 
diff --git a/plasmapy/formulary/parameters.py b/plasmapy/formulary/parameters.py
index fd594b7a..b517b720 100644
--- a/plasmapy/formulary/parameters.py
+++ b/plasmapy/formulary/parameters.py
@@ -52,6 +52,7 @@
 from typing import Optional
 
 from plasmapy import particles
+from plasmapy.particles import Particle
 from plasmapy.utils import PhysicsError
 from plasmapy.utils.decorators import (
     angular_freq_to_hz,
@@ -61,12 +62,12 @@
 from plasmapy.utils.exceptions import PhysicsWarning
 
 
-def _grab_charge(ion, z_mean=None):
+def _grab_charge(ion: Particle, z_mean=None):
     """Utility function to merge two possible inputs for particle charge.
 
     Parameters
     ----------
-    ion : str or `plasmapy.particles.Particle`
+    ion : ~plasmapy.particles.Particle
         a string representing a charged particle, or a Particle object.
 
     z_mean : float
@@ -95,7 +96,7 @@ def _grab_charge(ion, z_mean=None):
 )
 def mass_density(
     density: [u.m ** -3, u.kg / (u.m ** 3)],
-    particle: Optional[str] = None,
+    particle: Optional[Particle] = None,
     z_mean: Optional[numbers.Real] = None,
 ) -> u.kg / u.m ** 3:
     """Utility function to merge two possible inputs for particle charge.
@@ -108,11 +109,12 @@ def mass_density(
         Either a particle density (number of particles per unit volume, in units
         of 1/m^3) or a mass density (in units of kg/m^3 or equivalent).
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle, optional
         Representation of the particle species (e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or `'He-4 +1'` for singly ionized helium-4).  If no
+        charge state information is provided, then the particles are assumed
+        to be singly charged. A particle is required if using a number
+        density (1/m^3) for the density parameter.
 
     z_mean : float
         An optional float describing the average ionization of a particle
@@ -161,7 +163,10 @@ def mass_density(
 @check_relativistic
 @validate_quantities(density={"can_be_negative": False})
 def Alfven_speed(
-    B: u.T, density: [u.m ** -3, u.kg / u.m ** 3], ion="p+", z_mean=None
+    B: u.T,
+    density: [u.m ** -3, u.kg / u.m ** 3],
+    ion: Optional[Particle] = None,
+    z_mean=None,
 ) -> u.m / u.s:
     r"""
     Return the Alfvén speed.
@@ -177,12 +182,13 @@ def Alfven_speed(
         Either the ion number density in units convertible to 1 / m**3,
         or the mass density in units convertible to kg / m**3.
 
-    ion : str, optional
+    ion : ~plasmapy.particles.Particle, optional
         Representation of the ion species (e.g., `'p'` for protons,
         `'D+'` for deuterium, or `'He-4 +1'` for singly ionized
-        helium-4), which defaults to protons.  If no charge state
-        information is provided, then the ions are assumed to be
-        singly charged.
+        helium-4). If no charge state information is provided, then the
+        ions are assumed to be singly charged. If the density is an ion
+        number density, then this paramter is required in order to convert
+        to mass density.
 
     z_mean : ~astropy.units.Quantity, optional
         The average ionization (arithmetic mean) for a plasma where the
@@ -271,11 +277,11 @@ def Alfven_speed(
 def ion_sound_speed(
     T_e: u.K,
     T_i: u.K,
+    ion: Particle,
     n_e: u.m ** -3 = None,
     k: u.m ** -1 = None,
     gamma_e=1,
     gamma_i=3,
-    ion="p+",
     z_mean=None,
 ) -> u.m / u.s:
     r"""
@@ -295,6 +301,12 @@ def ion_sound_speed(
         particle.  If this is not given, then the ion temperature is
         assumed to be zero.
 
+    ion : ~plasmapy.particles.Particle
+        Representation of the ion species (e.g., `'p'` for protons,
+        `'D+'` for deuterium, or 'He-4 +1' for singly ionized
+        helium-4). If no charge state information is provided, then the
+        ions are assumed to be singly charged.
+
     n_e : ~astropy.units.Quantity
         Electron number density. If this is not given, then ion_sound_speed
         will be approximated in the non-dispersive limit
@@ -318,13 +330,6 @@ def ion_sound_speed(
         assumes that ion motion has only one degree of freedom, namely
         along magnetic field lines.
 
-    ion : str, optional
-        Representation of the ion species (e.g., `'p'` for protons,
-        `'D+'` for deuterium, or 'He-4 +1' for singly ionized
-        helium-4), which defaults to protons.  If no charge state
-        information is provided, then the ions are assumed to be
-        singly charged.
-
     z_mean : ~astropy.units.Quantity, optional
         The average ionization (arithmetic mean) for a plasma where the
         a macroscopic description is valid. If this quantity is not
@@ -401,7 +406,7 @@ def ion_sound_speed(
     <Quantity 203155... m / s>
     >>> ion_sound_speed(T_e=5e6*u.K, T_i=0*u.K, n_e=n, k=k_2, ion='p', gamma_e=1, gamma_i=3)
     <Quantity 310.31... m / s>
-    >>> ion_sound_speed(T_e=5e6*u.K, T_i=0*u.K, n_e=n, k=k_1)
+    >>> ion_sound_speed(T_e=5e6*u.K, T_i=0*u.K, n_e=n, k=k_1, ion='p')
     <Quantity 203155... m / s>
     >>> ion_sound_speed(T_e=500*u.eV, T_i=200*u.eV, n_e=n, k=k_1, ion='D+')
     <Quantity 229585... m / s>
@@ -468,7 +473,7 @@ def ion_sound_speed(
 @particles.particle_input
 def thermal_speed(
     T: u.K,
-    particle: particles.Particle = "e-",
+    particle: Particle,
     method="most_probable",
     mass: u.kg = np.nan * u.kg,
     ndim=3,
@@ -484,11 +489,11 @@ def thermal_speed(
     T : ~astropy.units.Quantity
         The particle temperature in either kelvin or energy per particle
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or `'He-4 +1'` for singly ionized helium-4). If no
+        charge state information is provided, then the particles are
+        assumed to be singly charged.
 
     method : str, optional
         Method to be used for calculating the thermal speed. Options are
@@ -579,9 +584,9 @@ def thermal_speed(
     <Quantity 132620... m / s>
     >>> thermal_speed(1e6*u.K, particle='e-')
     <Quantity 550569... m / s>
-    >>> thermal_speed(1e6*u.K, method="rms")
+    >>> thermal_speed(1e6*u.K, "e-", method="rms")
     <Quantity 674307... m / s>
-    >>> thermal_speed(1e6*u.K, method="mean_magnitude")
+    >>> thermal_speed(1e6*u.K, "e-", method="mean_magnitude")
     <Quantity 621251... m / s>
 
     """
@@ -591,7 +596,7 @@ def thermal_speed(
     try:
         coef = _coefficients[ndim]
     except KeyError:
-        raise ValueError("{ndim} is not a supported value for ndim in " "thermal_speed")
+        raise ValueError("{ndim} is not a supported value for ndim in thermal_speed")
     try:
         coef = coef[method]
     except KeyError:
@@ -664,7 +669,7 @@ def thermal_pressure(T: u.K, n: u.m ** -3) -> u.Pa:
     T={"can_be_negative": False, "equivalencies": u.temperature_energy()}
 )
 def kappa_thermal_speed(
-    T: u.K, kappa, particle="e-", method="most_probable"
+    T: u.K, kappa, particle: Particle, method="most_probable"
 ) -> u.m / u.s:
     r"""Return the most probable speed for a particle within a Kappa
     distribution.
@@ -682,11 +687,11 @@ def kappa_thermal_speed(
         of the Kappa velocity distribution function. Kappa must be greater
         than 3/2.
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., 'p' for protons, 'D+'
-        for deuterium, or 'He-4 +1' for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or 'He-4 +1' for singly ionized helium-4). If no
+        charge state information is provided, then the particles are
+        assumed to be singly charged.
 
     method : str, optional
         Method to be used for calculating the thermal speed. Options are
@@ -780,8 +785,8 @@ def Hall_parameter(
     n: u.m ** -3,
     T: u.K,
     B: u.T,
-    ion,
-    particle="e-",
+    ion: Particle,
+    particle: Particle,
     coulomb_log=None,
     V=None,
     coulomb_log_method="classical",
@@ -801,10 +806,10 @@ def Hall_parameter(
         The temperature of particles
     B : ~astropy.units.quantity.Quantity
         The magnetic field
-    ion : str
-        String signifying the type of ion.
-    particle : str, optional
-        String signifying the type of particles. Defaults to electrons.
+    ion : ~plasmapy.particles.Particle
+        The type of ion ``particle`` is colliding with.
+    particle : ~plasmapy.particles.Particle
+        The particle of interest.
     coulomb_log : float, optional
         Preset value for the Coulomb logarithm. Used mostly for testing purposes.
     V : ~astropy.units.quantity.Quantity
@@ -825,9 +830,9 @@ def Hall_parameter(
     Examples
     --------
     >>> from astropy import units as u
-    >>> Hall_parameter(1e10 * u.m**-3, 2.8e3 * u.eV, 2.3 * u.T, 'He-4 +1')
+    >>> Hall_parameter(1e10 * u.m**-3, 2.8e3 * u.eV, 2.3 * u.T, 'He-4 +1', 'e-')
     <Quantity 7.26446...e+16>
-    >>> Hall_parameter(1e10 * u.m**-3, 5.8e3 * u.eV, 2.3 * u.T, 'He-4 +1')
+    >>> Hall_parameter(1e10 * u.m**-3, 5.8e3 * u.eV, 2.3 * u.T, 'He-4 +1', 'e-')
     <Quantity 2.11158...e+17>
 
     """
@@ -858,7 +863,7 @@ def Hall_parameter(
     }
 )
 @angular_freq_to_hz
-def gyrofrequency(B: u.T, particle="e-", signed=False, Z=None) -> u.rad / u.s:
+def gyrofrequency(B: u.T, particle: Particle, signed=False, Z=None) -> u.rad / u.s:
     r"""Calculate the particle gyrofrequency in units of radians per second.
 
     **Aliases:** `oc_`, `wc_`
@@ -868,11 +873,11 @@ def gyrofrequency(B: u.T, particle="e-", signed=False, Z=None) -> u.rad / u.s:
     B : ~astropy.units.Quantity
         The magnetic field magnitude in units convertible to tesla.
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., 'p' for protons, 'D+'
-        for deuterium, or 'He-4 +1' for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or 'He-4 +1' for singly ionized helium-4). If no
+        charge state information is provided, then the particles are assumed
+        to be singly charged.
 
     signed : bool, optional
         The gyrofrequency can be defined as signed (negative for electron,
@@ -929,11 +934,11 @@ def gyrofrequency(B: u.T, particle="e-", signed=False, Z=None) -> u.rad / u.s:
     Examples
     --------
     >>> from astropy import units as u
-    >>> gyrofrequency(0.1*u.T)
+    >>> gyrofrequency(0.1*u.T, 'e-')
     <Quantity 1.7588...e+10 rad / s>
-    >>> gyrofrequency(0.1*u.T, to_hz=True)
+    >>> gyrofrequency(0.1*u.T, 'e-', to_hz=True)
     <Quantity 2.79924...e+09 Hz>
-    >>> gyrofrequency(0.1*u.T, signed=True)
+    >>> gyrofrequency(0.1*u.T, 'e-', signed=True)
     <Quantity -1.75882...e+10 rad / s>
     >>> gyrofrequency(0.01*u.T, 'p')
     <Quantity 957883.32... rad / s>
@@ -943,7 +948,7 @@ def gyrofrequency(B: u.T, particle="e-", signed=False, Z=None) -> u.rad / u.s:
     <Quantity 319964.5... rad / s>
     >>> gyrofrequency(0.01*u.T, particle='T+', to_hz=True)
     <Quantity 50923.9... Hz>
-    >>> omega_ce = gyrofrequency(0.1*u.T)
+    >>> omega_ce = gyrofrequency(0.1*u.T, 'e-')
     >>> print(omega_ce)
     1758820... rad / s
     >>> f_ce = omega_ce.to(u.Hz, equivalencies=[(u.cy/u.s, u.Hz)])
@@ -975,7 +980,7 @@ def gyrofrequency(B: u.T, particle="e-", signed=False, Z=None) -> u.rad / u.s:
 )
 def gyroradius(
     B: u.T,
-    particle="e-",
+    particle: Particle,
     *,
     Vperp: u.m / u.s = np.nan * u.m / u.s,
     T_i: u.K = np.nan * u.K,
@@ -989,11 +994,11 @@ def gyroradius(
     B : ~astropy.units.Quantity
         The magnetic field magnitude in units convertible to tesla.
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or `'He-4 +1'` for singly ionized helium-4).  If no
+        charge state information is provided, then the particles are assumed
+        to be singly charged.
 
     Vperp : ~astropy.units.Quantity, optional
         The component of particle velocity that is perpendicular to the
@@ -1052,23 +1057,23 @@ def gyroradius(
     Examples
     --------
     >>> from astropy import units as u
-    >>> gyroradius(0.2*u.T,particle='p+',T_i=1e5*u.K)
+    >>> gyroradius(0.2*u.T, particle='p+', T_i=1e5*u.K)
     <Quantity 0.002120... m>
-    >>> gyroradius(0.2*u.T,particle='p+',T_i=1e5*u.K)
+    >>> gyroradius(0.2*u.T, particle='p+', T_i=1e5*u.K)
     <Quantity 0.002120... m>
-    >>> gyroradius(5*u.uG,particle='alpha',T_i=1*u.eV)
+    >>> gyroradius(5*u.uG, particle='alpha', T_i=1*u.eV)
     <Quantity 288002.38... m>
-    >>> gyroradius(400*u.G,particle='Fe+++',Vperp=1e7*u.m/u.s)
+    >>> gyroradius(400*u.G, particle='Fe+++', Vperp=1e7*u.m/u.s)
     <Quantity 48.23129... m>
-    >>> gyroradius(B=0.01*u.T,T_i=1e6*u.K)
+    >>> gyroradius(B=0.01*u.T, particle='e-', T_i=1e6*u.K)
     <Quantity 0.003130... m>
-    >>> gyroradius(B=0.01*u.T,Vperp=1e6*u.m/u.s)
+    >>> gyroradius(0.01*u.T, 'e-', Vperp=1e6*u.m/u.s)
     <Quantity 0.000568... m>
-    >>> gyroradius(0.2*u.T,T_i=1e5*u.K)
+    >>> gyroradius(0.2*u.T, 'e-', T_i=1e5*u.K)
     <Quantity 4.94949...e-05 m>
-    >>> gyroradius(5*u.uG,T_i=1*u.eV)
+    >>> gyroradius(5*u.uG, 'e-', T_i=1*u.eV)
     <Quantity 6744.25... m>
-    >>> gyroradius(400*u.G,Vperp=1e7*u.m/u.s)
+    >>> gyroradius(400*u.G, 'e-', Vperp=1e7*u.m/u.s)
     <Quantity 0.001421... m>
 
     """
@@ -1139,7 +1144,7 @@ def gyroradius(
     },
 )
 @angular_freq_to_hz
-def plasma_frequency(n: u.m ** -3, particle="e-", z_mean=None) -> u.rad / u.s:
+def plasma_frequency(n: u.m ** -3, particle: Particle, z_mean=None) -> u.rad / u.s:
     r"""Calculate the particle plasma frequency.
 
     **Aliases:** `wp_`
@@ -1149,11 +1154,11 @@ def plasma_frequency(n: u.m ** -3, particle="e-", z_mean=None) -> u.rad / u.s:
     n : ~astropy.units.Quantity
         Particle number density in units convertible to per cubic meter
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., 'p' for protons, 'D+'
-        for deuterium, or 'He-4 +1' for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        for deuterium, or 'He-4 +1' for singly ionized helium-4). If no
+        charge state information is provided, then the particles are assumed
+        to be singly charged.
 
     z_mean : ~astropy.units.Quantity, optional
         The average ionization (arithmetic mean) for a plasma where the
@@ -1208,9 +1213,9 @@ def plasma_frequency(n: u.m ** -3, particle="e-", z_mean=None) -> u.rad / u.s:
     <Quantity 6.62608...e+08 Hz>
     >>> plasma_frequency(1e19*u.m**-3, particle='D+')
     <Quantity 2.94462...e+09 rad / s>
-    >>> plasma_frequency(1e19*u.m**-3)
+    >>> plasma_frequency(1e19*u.m**-3, 'e-')
     <Quantity 1.78398...e+11 rad / s>
-    >>> plasma_frequency(1e19*u.m**-3, to_hz=True)
+    >>> plasma_frequency(1e19*u.m**-3, 'e-', to_hz=True)
     <Quantity 2.83930...e+10 Hz>
 
     """
@@ -1393,7 +1398,7 @@ def Debye_number(T_e: u.K, n_e: u.m ** -3) -> u.dimensionless_unscaled:
     validations_on_return={"equivalencies": u.dimensionless_angles()},
 )
 @particles.particle_input(require="charged")
-def inertial_length(n: u.m ** -3, particle: particles.Particle) -> u.m:
+def inertial_length(n: u.m ** -3, particle: Particle) -> u.m:
     r"""
     Calculate a charged particle's inertial length.
 
@@ -1404,7 +1409,7 @@ def inertial_length(n: u.m ** -3, particle: particles.Particle) -> u.m:
     n : ~astropy.units.Quantity
         Particle number density in units convertible to m ** -3.
 
-    particle : str, optional
+    particle : ~plasmapy.particles.Particle
         Representation of the particle species (e.g., 'p+' for protons,
         'D+' for deuterium, or 'He-4 +1' for singly ionized helium-4).
 
@@ -1657,8 +1662,8 @@ def upper_hybrid_frequency(B: u.T, n_e: u.m ** -3) -> u.rad / u.s:
     <Quantity 6.37350...e+10 Hz>
 
     """
-    omega_pe = plasma_frequency(n=n_e)
-    omega_ce = gyrofrequency(B)
+    omega_pe = plasma_frequency(n=n_e, particle="e-")
+    omega_ce = gyrofrequency(B, "e-")
     omega_uh = np.sqrt(omega_pe ** 2 + omega_ce ** 2)
 
     return omega_uh
@@ -1676,7 +1681,7 @@ def upper_hybrid_frequency(B: u.T, n_e: u.m ** -3) -> u.rad / u.s:
     },
 )
 @angular_freq_to_hz
-def lower_hybrid_frequency(B: u.T, n_i: u.m ** -3, ion="p+") -> u.rad / u.s:
+def lower_hybrid_frequency(B: u.T, n_i: u.m ** -3, ion: Particle) -> u.rad / u.s:
     r"""
     Return the lower hybrid frequency.
 
@@ -1690,11 +1695,11 @@ def lower_hybrid_frequency(B: u.T, n_i: u.m ** -3, ion="p+") -> u.rad / u.s:
     n_i : ~astropy.units.Quantity
         Ion number density.
 
-    ion : str, optional
+    ion : ~plasmapy.particles.Particle
         Representation of the ion species (e.g., 'p' for protons, 'D+'
-        for deuterium, or 'He-4 +1' for singly ionized helium-4),
-        which defaults to protons.  If no charge state information is
-        provided, then the ions are assumed to be singly charged.
+        for deuterium, or 'He-4 +1' for singly ionized helium-4). If no
+        charge state information is provided, then the ions are assumed to
+        be singly charged.
 
     Returns
     -------
@@ -1751,7 +1756,7 @@ def lower_hybrid_frequency(B: u.T, n_i: u.m ** -3, ion="p+") -> u.rad / u.s:
 
     omega_ci = gyrofrequency(B, particle=ion)
     omega_pi = plasma_frequency(n_i, particle=ion)
-    omega_ce = gyrofrequency(B)
+    omega_ce = gyrofrequency(B, particle="e-")
     omega_lh = ((omega_ci * omega_ce) ** -1 + omega_pi ** -2) ** -0.5
     # TODO possibly optimize the above line via np.sqrt
     omega_lh = omega_lh


Test patch:
diff --git a/plasmapy/formulary/tests/test_parameters.py b/plasmapy/formulary/tests/test_parameters.py
index b9790bac..9182de07 100644
--- a/plasmapy/formulary/tests/test_parameters.py
+++ b/plasmapy/formulary/tests/test_parameters.py
@@ -116,23 +116,23 @@ def test_Alfven_speed():
     # TODO: break this test up until multiple tests
 
     assert np.isclose(
-        Alfven_speed(1 * u.T, 1e-8 * u.kg * u.m ** -3).value,
+        Alfven_speed(1 * u.T, 1e-8 * u.kg * u.m ** -3, "p").value,
         8920620.580763856,
         rtol=1e-6,
     )
 
-    V_A = Alfven_speed(B, n_i)
+    V_A = Alfven_speed(B, n_i, "p")
     assert np.isclose(V_A.value, (B / np.sqrt(mu0 * n_i * (m_p + m_e))).si.value)
 
-    assert Alfven_speed(B, rho) == Alfven_speed(B, n_i)
+    assert Alfven_speed(B, rho, "p") == Alfven_speed(B, n_i, "p")
 
-    assert Alfven_speed(B, rho).unit.is_equivalent(u.m / u.s)
+    assert Alfven_speed(B, rho, "p").unit.is_equivalent(u.m / u.s)
 
-    assert Alfven_speed(B, rho) == Alfven_speed(-B, rho)
+    assert Alfven_speed(B, rho, "p") == Alfven_speed(-B, rho, "p")
 
-    assert Alfven_speed(B, 4 * rho) == 0.5 * Alfven_speed(B, rho)
+    assert Alfven_speed(B, 4 * rho, "p") == 0.5 * Alfven_speed(B, rho, "p")
 
-    assert Alfven_speed(2 * B, rho) == 2 * Alfven_speed(B, rho)
+    assert Alfven_speed(2 * B, rho, "p") == 2 * Alfven_speed(B, rho, "p")
 
     # Case when Z=1 is assumed
     with pytest.warns(RelativityWarning):
@@ -144,26 +144,26 @@ def test_Alfven_speed():
         )
 
     # Case where magnetic field and density are Quantity arrays
-    V_A_arr = Alfven_speed(B_arr, rho_arr)
-    V_A_arr0 = Alfven_speed(B_arr[0], rho_arr[0])
-    V_A_arr1 = Alfven_speed(B_arr[1], rho_arr[1])
+    V_A_arr = Alfven_speed(B_arr, rho_arr, "p")
+    V_A_arr0 = Alfven_speed(B_arr[0], rho_arr[0], "p")
+    V_A_arr1 = Alfven_speed(B_arr[1], rho_arr[1], "p")
     assert np.isclose(V_A_arr0.value, V_A_arr[0].value)
     assert np.isclose(V_A_arr1.value, V_A_arr[1].value)
 
     # Case where magnetic field is an array but density is a scalar Quantity
-    V_A_arr = Alfven_speed(B_arr, rho)
-    V_A_arr0 = Alfven_speed(B_arr[0], rho)
-    V_A_arr1 = Alfven_speed(B_arr[1], rho)
+    V_A_arr = Alfven_speed(B_arr, rho, "p")
+    V_A_arr0 = Alfven_speed(B_arr[0], rho, "p")
+    V_A_arr1 = Alfven_speed(B_arr[1], rho, "p")
     assert np.isclose(V_A_arr0.value, V_A_arr[0].value)
     assert np.isclose(V_A_arr1.value, V_A_arr[1].value)
 
     with pytest.raises(ValueError):
-        Alfven_speed(np.array([5, 6, 7]) * u.T, np.array([5, 6]) * u.m ** -3)
+        Alfven_speed(np.array([5, 6, 7]) * u.T, np.array([5, 6]) * u.m ** -3, "p")
 
-    assert np.isnan(Alfven_speed(B_nanarr, rho_arr)[1])
+    assert np.isnan(Alfven_speed(B_nanarr, rho_arr, "p")[1])
 
     with pytest.raises(ValueError):
-        Alfven_speed(B_arr, rho_negarr)
+        Alfven_speed(B_arr, rho_negarr, "p")
 
     with pytest.raises(u.UnitTypeError):
         Alfven_speed(5 * u.A, n_i, ion="p")
@@ -197,7 +197,7 @@ def test_Alfven_speed():
         assert Alfven_speed(-np.inf * u.T, 1 * u.m ** -3, ion="p") == np.inf * u.m / u.s
 
     with pytest.warns(u.UnitsWarning):
-        assert Alfven_speed(1.0, n_i) == Alfven_speed(1.0 * u.T, n_i)
+        assert Alfven_speed(1.0, n_i, "p") == Alfven_speed(1.0 * u.T, n_i, "p")
 
     Alfven_speed(1 * u.T, 5e19 * u.m ** -3, ion="p")
     # testing for user input z_mean
@@ -219,6 +219,11 @@ def test_ion_sound_speed():
         218816.06086407552,
     )
 
+    # Test that function call without keyword argument works correctly
+    assert np.isclose(
+        ion_sound_speed(1.831 * u.MK, 1.3232 * u.MK, "p").value, 218816.06086407552,
+    )
+
     assert np.isclose(
         ion_sound_speed(
             T_i=1.3232 * u.MK,
@@ -267,13 +272,13 @@ def test_ion_sound_speed():
     ).unit.is_equivalent(u.m / u.s)
 
     with pytest.raises(RelativityError):
-        ion_sound_speed(T_i=T_i, T_e=T_e, n_e=n_e, k=k_1, gamma_i=np.inf)
+        ion_sound_speed(T_i=T_i, T_e=T_e, n_e=n_e, k=k_1, ion="p", gamma_i=np.inf)
 
     with pytest.warns(PhysicsWarning):
-        ion_sound_speed(T_i=T_i, T_e=T_e, n_e=n_e)
+        ion_sound_speed(T_i=T_i, T_e=T_e, n_e=n_e, ion="p")
 
     with pytest.warns(PhysicsWarning):
-        ion_sound_speed(T_i=T_i, T_e=T_e, k=k_1)
+        ion_sound_speed(T_i=T_i, T_e=T_e, k=k_1, ion="p")
 
     with pytest.raises(u.UnitTypeError):
         ion_sound_speed(
@@ -281,71 +286,76 @@ def test_ion_sound_speed():
             T_e=np.array([3, 4]) * u.K,
             n_e=np.array([5, 6, 5]) * u.m ** -3,
             k=np.array([3, 4]) * u.m ** -3,
+            ion="p",
         )
 
     with pytest.raises(TypeError):  # Is this test right??????
-        ion_sound_speed(5 * u.T)
+        ion_sound_speed(5 * u.T, ion="p")
 
     with pytest.raises(TypeError):
-        ion_sound_speed("p")
+        ion_sound_speed(ion="p")
 
     with pytest.raises(PhysicsError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_i=0.9999)
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_i=0.9999, ion="p")
 
     with pytest.raises(PhysicsError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_e=0.9999)
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_e=0.9999, ion="p")
 
     with pytest.raises(TypeError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_e="sdjklsf")
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_e="sdjklsf", ion="p")
 
     with pytest.raises(TypeError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_i="fsdfas")
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, gamma_i="fsdfas", ion="p")
 
     with pytest.raises(InvalidParticleError):
         ion_sound_speed(T_i=T_i, T_e=0 * u.K, ion="cupcakes")
 
     with pytest.raises(ValueError):
-        ion_sound_speed(T_i=-np.abs(T_i), T_e=0 * u.K)
+        ion_sound_speed(T_i=-np.abs(T_i), T_e=0 * u.K, ion="p")
 
     with pytest.raises(ValueError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, n_e=-np.abs(n_e), k=k_1)
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, n_e=-np.abs(n_e), k=k_1, ion="p")
 
     with pytest.raises(ValueError):
-        ion_sound_speed(T_i=T_i, T_e=0 * u.K, n_e=n_e, k=-np.abs(k_1))
+        ion_sound_speed(T_i=T_i, T_e=0 * u.K, n_e=n_e, k=-np.abs(k_1), ion="p")
 
     with pytest.warns(RelativityWarning):
-        ion_sound_speed(T_i=5e11 * u.K, T_e=0 * u.K)
+        ion_sound_speed(T_i=5e11 * u.K, T_e=0 * u.K, ion="p")
 
     with pytest.raises(RelativityError):
-        ion_sound_speed(T_i=5e19 * u.K, T_e=0 * u.K)
+        ion_sound_speed(T_i=5e19 * u.K, T_e=0 * u.K, ion="p")
 
     with pytest.raises(u.UnitTypeError):
-        ion_sound_speed(T_i=5 * u.A, T_e=0 * u.K, n_e=n_e, k=k_1)
+        ion_sound_speed(T_i=5 * u.A, T_e=0 * u.K, n_e=n_e, k=k_1, ion="p")
 
-    assert np.isnan(ion_sound_speed(T_i=T_nanarr, T_e=0 * u.K, n_e=n_e, k=k_1)[1])
+    assert np.isnan(
+        ion_sound_speed(T_i=T_nanarr, T_e=0 * u.K, n_e=n_e, k=k_1, ion="p")[1]
+    )
 
-    assert np.isnan(ion_sound_speed(T_e=T_nanarr, T_i=0 * u.K, n_e=n_e, k=k_1)[1])
+    assert np.isnan(
+        ion_sound_speed(T_e=T_nanarr, T_i=0 * u.K, n_e=n_e, k=k_1, ion="p")[1]
+    )
 
     with pytest.raises(ValueError):
-        ion_sound_speed(T_i=T_negarr, T_e=0 * u.K, n_e=n_e, k=k_1)
+        ion_sound_speed(T_i=T_negarr, T_e=0 * u.K, n_e=n_e, k=k_1, ion="p")
 
     with pytest.raises(ValueError):
-        ion_sound_speed(T_e=T_negarr, T_i=0 * u.K, n_e=n_e, k=k_1)
+        ion_sound_speed(T_e=T_negarr, T_i=0 * u.K, n_e=n_e, k=k_1, ion="p")
 
     with pytest.warns(u.UnitsWarning):
         assert ion_sound_speed(
-            T_e=1.2e6, T_i=0 * u.K, n_e=n_e, k=k_1
-        ) == ion_sound_speed(T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=k_1)
+            T_e=1.2e6, T_i=0 * u.K, n_e=n_e, k=k_1, ion="p"
+        ) == ion_sound_speed(T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=k_1, ion="p")
 
     with pytest.warns(u.UnitsWarning):
         assert ion_sound_speed(
-            T_i=1.3e6, T_e=0 * u.K, n_e=n_e, k=k_1
-        ) == ion_sound_speed(T_i=1.3e6 * u.K, T_e=0 * u.K, n_e=n_e, k=k_1)
+            T_i=1.3e6, T_e=0 * u.K, n_e=n_e, k=k_1, ion="p"
+        ) == ion_sound_speed(T_i=1.3e6 * u.K, T_e=0 * u.K, n_e=n_e, k=k_1, ion="p")
 
-    ion_sound_speed(T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=k_1)
+    ion_sound_speed(T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=k_1, ion="p")
     # testing for user input z_mean
     testMeth1 = ion_sound_speed(
-        T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=0 * u.m ** -1, z_mean=0.8
+        T_e=1.2e6 * u.K, T_i=0 * u.K, n_e=n_e, k=0 * u.m ** -1, z_mean=0.8, ion="p"
     ).si.value
     testTrue1 = 89018.09
     errStr = f"ion_sound_speed() gave {testMeth1}, should be {testTrue1}."
@@ -356,28 +366,28 @@ def test_ion_sound_speed():
 
 def test_thermal_speed():
     r"""Test the thermal_speed function in parameters.py"""
-    assert thermal_speed(T_e).unit.is_equivalent(u.m / u.s)
+    assert thermal_speed(T_e, "e-").unit.is_equivalent(u.m / u.s)
 
-    assert thermal_speed(T_e) > thermal_speed(T_e, "p")
+    assert thermal_speed(T_e, "e-") > thermal_speed(T_e, "p")
 
     # The NRL Plasma Formulary uses a definition of the electron
     # thermal speed that differs by a factor of sqrt(2).
-    assert np.isclose(thermal_speed(1 * u.MK).value, 5505694.743141063)
+    assert np.isclose(thermal_speed(1 * u.MK, "e-").value, 5505694.743141063)
 
     with pytest.raises(u.UnitTypeError):
-        thermal_speed(5 * u.m)
+        thermal_speed(5 * u.m, "e-")
 
     with pytest.raises(ValueError):
-        thermal_speed(-T_e)
+        thermal_speed(-T_e, "e-")
 
     with pytest.warns(RelativityWarning):
-        thermal_speed(1e9 * u.K)
+        thermal_speed(1e9 * u.K, "e-")
 
     with pytest.raises(RelativityError):
-        thermal_speed(5e19 * u.K)
+        thermal_speed(5e19 * u.K, "e-")
 
     with pytest.warns(u.UnitsWarning):
-        assert thermal_speed(1e5) == thermal_speed(1e5 * u.K)
+        assert thermal_speed(1e5, "e-") == thermal_speed(1e5 * u.K, "e-")
 
     assert thermal_speed(T_i, particle="p").unit.is_equivalent(u.m / u.s)
 
@@ -389,47 +399,54 @@ def test_thermal_speed():
 
     # Explicitly check all three modes and dimensionalities
     # ndim = 1
-    assert np.isclose(thermal_speed(T_e, method="most_probable", ndim=1).si.value, 0.0)
+    assert np.isclose(
+        thermal_speed(T_e, "e-", method="most_probable", ndim=1).si.value, 0.0
+    )
 
     # Regression tests start here!
     assert np.isclose(
-        thermal_speed(T_e, method="rms", ndim=1).si.value, 3893114.2008620175
+        thermal_speed(T_e, "e-", method="rms", ndim=1).si.value, 3893114.2008620175
     )
 
     assert np.isclose(
-        thermal_speed(T_e, method="mean_magnitude", ndim=1).si.value, 3106255.714310189
+        thermal_speed(T_e, "e-", method="mean_magnitude", ndim=1).si.value,
+        3106255.714310189,
     )
 
     # ndim = 2
     assert np.isclose(
-        thermal_speed(T_e, method="most_probable", ndim=2).si.value, 3893114.2008620175
+        thermal_speed(T_e, "e-", method="most_probable", ndim=2).si.value,
+        3893114.2008620175,
     )
 
     assert np.isclose(
-        thermal_speed(T_e, method="rms", ndim=2).si.value, 5505694.902726359
+        thermal_speed(T_e, "e-", method="rms", ndim=2).si.value, 5505694.902726359
     )
 
     assert np.isclose(
-        thermal_speed(T_e, method="mean_magnitude", ndim=2).si.value, 4879295.066124102
+        thermal_speed(T_e, "e-", method="mean_magnitude", ndim=2).si.value,
+        4879295.066124102,
     )
 
     # ndim = 3
     assert np.isclose(
-        thermal_speed(T_e, method="most_probable", ndim=3).si.value, 5505694.902726359
+        thermal_speed(T_e, "e-", method="most_probable", ndim=3).si.value,
+        5505694.902726359,
     )
 
     assert np.isclose(
-        thermal_speed(T_e, method="rms", ndim=3).si.value, 6743071.595560921
+        thermal_speed(T_e, "e-", method="rms", ndim=3).si.value, 6743071.595560921
     )
 
     assert np.isclose(
-        thermal_speed(T_e, method="mean_magnitude", ndim=3).si.value, 6212511.428620378
+        thermal_speed(T_e, "e-", method="mean_magnitude", ndim=3).si.value,
+        6212511.428620378,
     )
 
     # Case when Z=1 is assumed
     assert thermal_speed(T_i, particle="p") == thermal_speed(T_i, particle="H-1+")
 
-    assert thermal_speed(1 * u.MK, particle="e+") == thermal_speed(1 * u.MK)
+    assert thermal_speed(1 * u.MK, particle="e+") == thermal_speed(1 * u.MK, "e-")
 
     with pytest.raises(u.UnitTypeError):
         thermal_speed(5 * u.m, particle="p")
@@ -452,20 +469,21 @@ def test_thermal_speed():
         )
 
     assert np.isclose(
-        thermal_speed(1e6 * u.K, method="mean_magnitude").si.value, 6212510.3969422
+        thermal_speed(1e6 * u.K, "e-", method="mean_magnitude").si.value,
+        6212510.3969422,
     )
 
     assert np.isclose(
-        thermal_speed(1e6 * u.K, method="rms").si.value, 6743070.475775486
+        thermal_speed(1e6 * u.K, "e-", method="rms").si.value, 6743070.475775486
     )
 
     # Test invalid method
     with pytest.raises(ValueError):
-        thermal_speed(T_i, method="sadks")
+        thermal_speed(T_i, "e-", method="sadks")
 
     # Test invalid ndim
     with pytest.raises(ValueError):
-        thermal_speed(T_i, ndim=4)
+        thermal_speed(T_i, "e-", ndim=4)
 
     assert_can_handle_nparray(thermal_speed)
 
@@ -564,37 +582,41 @@ def test_handle_nparrays(self, kwargs=None):
 def test_gyrofrequency():
     r"""Test the gyrofrequency function in parameters.py."""
 
-    assert gyrofrequency(B).unit.is_equivalent(u.rad / u.s)
+    assert gyrofrequency(B, "e-").unit.is_equivalent(u.rad / u.s)
 
-    assert gyrofrequency(B, to_hz=True).unit.is_equivalent(u.Hz)
+    assert gyrofrequency(B, "e-", to_hz=True).unit.is_equivalent(u.Hz)
 
-    assert np.isclose(gyrofrequency(1 * u.T).value, 175882008784.72018)
+    assert np.isclose(gyrofrequency(1 * u.T, "e-").value, 175882008784.72018)
 
-    assert np.isclose(gyrofrequency(2.4 * u.T).value, 422116821083.3284)
+    assert np.isclose(gyrofrequency(2.4 * u.T, "e-").value, 422116821083.3284)
 
-    assert np.isclose(gyrofrequency(1 * u.T, to_hz=True).value, 27992490076.528206)
+    assert np.isclose(
+        gyrofrequency(1 * u.T, "e-", to_hz=True).value, 27992490076.528206
+    )
 
-    assert np.isclose(gyrofrequency(2.4 * u.T, signed=True).value, -422116821083.3284)
+    assert np.isclose(
+        gyrofrequency(2.4 * u.T, "e-", signed=True).value, -422116821083.3284
+    )
 
-    assert np.isclose(gyrofrequency(1 * u.G).cgs.value, 1.76e7, rtol=1e-3)
+    assert np.isclose(gyrofrequency(1 * u.G, "e-").cgs.value, 1.76e7, rtol=1e-3)
 
     with pytest.raises(TypeError):
-        gyrofrequency(u.m)
+        gyrofrequency(u.m, "e-")
 
     with pytest.raises(u.UnitTypeError):
-        gyrofrequency(u.m * 1)
+        gyrofrequency(u.m * 1, "e-")
 
-    assert np.isnan(gyrofrequency(B_nanarr)[-1])
+    assert np.isnan(gyrofrequency(B_nanarr, "e-")[-1])
 
     # The following is a test to check that equivalencies from astropy
     # are working.
-    omega_ce = gyrofrequency(2.2 * u.T)
+    omega_ce = gyrofrequency(2.2 * u.T, "e-")
     f_ce = (omega_ce / (2 * np.pi)) / u.rad
     f_ce_use_equiv = omega_ce.to(u.Hz, equivalencies=[(u.cy / u.s, u.Hz)])
     assert np.isclose(f_ce.value, f_ce_use_equiv.value)
 
     with pytest.warns(u.UnitsWarning):
-        assert gyrofrequency(5.0) == gyrofrequency(5.0 * u.T)
+        assert gyrofrequency(5.0, "e-") == gyrofrequency(5.0 * u.T, "e-")
 
     assert gyrofrequency(B, particle=ion).unit.is_equivalent(u.rad / u.s)
 
@@ -609,7 +631,7 @@ def test_gyrofrequency():
     # Case when Z=1 is assumed
     # assert gyrofrequency(B, particle='p+') == gyrofrequency(B, particle='H-1')
 
-    assert gyrofrequency(B, particle="e+") == gyrofrequency(B)
+    assert gyrofrequency(B, particle="e+") == gyrofrequency(B, "e-")
 
     with pytest.warns(u.UnitsWarning):
         gyrofrequency(8, "p")
@@ -639,48 +661,50 @@ def test_gyrofrequency():
 def test_gyroradius():
     r"""Test the gyroradius function in parameters.py."""
 
-    assert gyroradius(B, T_i=T_e).unit.is_equivalent(u.m)
+    assert gyroradius(B, "e-", T_i=T_e).unit.is_equivalent(u.m)
 
-    assert gyroradius(B, Vperp=25 * u.m / u.s).unit.is_equivalent(u.m)
+    assert gyroradius(B, "e-", Vperp=25 * u.m / u.s).unit.is_equivalent(u.m)
 
     Vperp = 1e6 * u.m / u.s
     Bmag = 1 * u.T
-    omega_ce = gyrofrequency(Bmag)
+    omega_ce = gyrofrequency(Bmag, "e-")
     analytical_result = (Vperp / omega_ce).to(
         u.m, equivalencies=u.dimensionless_angles()
     )
-    assert gyroradius(Bmag, Vperp=Vperp) == analytical_result
+    assert gyroradius(Bmag, "e-", Vperp=Vperp) == analytical_result
 
     with pytest.raises(TypeError):
-        gyroradius(u.T)
+        gyroradius(u.T, "e-")
 
     with pytest.raises(u.UnitTypeError):
-        gyroradius(5 * u.A, Vperp=8 * u.m / u.s)
+        gyroradius(5 * u.A, "e-", Vperp=8 * u.m / u.s)
 
     with pytest.raises(u.UnitTypeError):
-        gyroradius(5 * u.T, Vperp=8 * u.m)
+        gyroradius(5 * u.T, "e-", Vperp=8 * u.m)
 
     with pytest.raises(ValueError):
-        gyroradius(np.array([5, 6]) * u.T, Vperp=np.array([5, 6, 7]) * u.m / u.s)
+        gyroradius(np.array([5, 6]) * u.T, "e-", Vperp=np.array([5, 6, 7]) * u.m / u.s)
 
-    assert np.isnan(gyroradius(np.nan * u.T, Vperp=1 * u.m / u.s))
+    assert np.isnan(gyroradius(np.nan * u.T, "e-", Vperp=1 * u.m / u.s))
 
     with pytest.raises(ValueError):
-        gyroradius(3.14159 * u.T, T_i=-1 * u.K)
+        gyroradius(3.14159 * u.T, "e-", T_i=-1 * u.K)
 
     with pytest.warns(u.UnitsWarning):
-        assert gyroradius(1.0, Vperp=1.0) == gyroradius(
-            1.0 * u.T, Vperp=1.0 * u.m / u.s
+        assert gyroradius(1.0, "e-", Vperp=1.0) == gyroradius(
+            1.0 * u.T, "e-", Vperp=1.0 * u.m / u.s
         )
 
     with pytest.warns(u.UnitsWarning):
-        assert gyroradius(1.1, T_i=1.2) == gyroradius(1.1 * u.T, T_i=1.2 * u.K)
+        assert gyroradius(1.1, "e-", T_i=1.2) == gyroradius(
+            1.1 * u.T, "e-", T_i=1.2 * u.K
+        )
 
     with pytest.raises(ValueError):
-        gyroradius(1.1 * u.T, Vperp=1 * u.m / u.s, T_i=1.2 * u.K)
+        gyroradius(1.1 * u.T, "e-", Vperp=1 * u.m / u.s, T_i=1.2 * u.K)
 
     with pytest.raises(u.UnitTypeError):
-        gyroradius(1.1 * u.T, Vperp=1.1 * u.m, T_i=1.2 * u.K)
+        gyroradius(1.1 * u.T, "e-", Vperp=1.1 * u.m, T_i=1.2 * u.K)
 
     assert gyroradius(B, particle="p", T_i=T_i).unit.is_equivalent(u.m)
 
@@ -713,7 +737,7 @@ def test_gyroradius():
     assert gyro_by_vperp == gyroradius(B2, particle="alpha", T_i=T2)
 
     explicit_positron_gyro = gyroradius(1 * u.T, particle="positron", T_i=1 * u.MK)
-    assert explicit_positron_gyro == gyroradius(1 * u.T, T_i=1 * u.MK)
+    assert explicit_positron_gyro == gyroradius(1 * u.T, "e-", T_i=1 * u.MK)
 
     with pytest.raises(TypeError):
         gyroradius(u.T, particle="p", Vperp=8 * u.m / u.s)
@@ -746,53 +770,57 @@ class Test_gyroradius:
     # some custom numpy array tests here, because of the T_i / Vperp situation
     def test_handle_numpy_array(self):
         # Tests to verify that can handle Quantities with numpy array as the value:
-        assert gyroradius(B_arr, Vperp=V_arr)[0] == gyroradius(B_arr[0], Vperp=V_arr[0])
-        assert gyroradius(B_arr, T_i=T_arr)[0] == gyroradius(B_arr[0], T_i=T_arr[0])
+        assert gyroradius(B_arr, "e-", Vperp=V_arr)[0] == gyroradius(
+            B_arr[0], "e-", Vperp=V_arr[0]
+        )
+        assert gyroradius(B_arr, "e-", T_i=T_arr)[0] == gyroradius(
+            B_arr[0], "e-", T_i=T_arr[0]
+        )
 
     def test_handle_mixed_Qarrays(self):
         # If both Vperp or Ti are input as Qarrays, but only one of the two is valid
         # at each element, then that's fine, the function should work:
-        assert gyroradius(B_arr, Vperp=V_nanarr, T_i=T_nanarr2)[0] == gyroradius(
-            B_arr[0], Vperp=V_nanarr[0], T_i=T_nanarr2[0]
+        assert gyroradius(B_arr, "e-", Vperp=V_nanarr, T_i=T_nanarr2)[0] == gyroradius(
+            B_arr[0], "e-", Vperp=V_nanarr[0], T_i=T_nanarr2[0]
         )
 
     def test_raise_two_valid_inputs(self):
         # If both Vperp or Ti are nan-less, Qarrays or not, should raise ValueError:
         with pytest.raises(ValueError):
-            gyroradius(B_arr, Vperp=V, T_i=T_arr)
+            gyroradius(B_arr, "e-", Vperp=V, T_i=T_arr)
         with pytest.raises(ValueError):
-            gyroradius(B_arr, Vperp=V_arr, T_i=T_i)
+            gyroradius(B_arr, "e-", Vperp=V_arr, T_i=T_i)
 
     def test_all_valid_and_one_valid(self):
         # If one of (Vperp, Ti) is a valid and one is Qarray with at least one valid, ValueError:
         with pytest.raises(ValueError):
-            gyroradius(B_arr, Vperp=V, T_i=T_nanarr)
+            gyroradius(B_arr, "e-", Vperp=V, T_i=T_nanarr)
         with pytest.raises(ValueError):
-            gyroradius(B_arr, Vperp=V_nanarr, T_i=T_i)
+            gyroradius(B_arr, "e-", Vperp=V_nanarr, T_i=T_i)
 
     def test_scalar_and_nan_qarray(self):
         # If either Vperp or Ti is a valid scalar and the other is a Qarray of all nans,
         # should do something valid and not raise a ValueError
-        assert np.all(np.isfinite(gyroradius(B_arr, Vperp=V, T_i=T_allnanarr)))
-        assert np.all(np.isfinite(gyroradius(B_arr, Vperp=V_allnanarr, T_i=T_i)))
+        assert np.all(np.isfinite(gyroradius(B_arr, "e-", Vperp=V, T_i=T_allnanarr)))
+        assert np.all(np.isfinite(gyroradius(B_arr, "e-", Vperp=V_allnanarr, T_i=T_i)))
 
     def test_keeps_arguments_unchanged(self):
         Vperp1 = u.Quantity([np.nan, 1], unit=u.m / u.s)
         Vperp2 = u.Quantity([np.nan, 1], unit=u.m / u.s)  # an exact copy
         T_i = u.Quantity([1, np.nan], unit=u.K)
 
-        gyroradius(B_arr, Vperp=Vperp1, T_i=T_i)
+        gyroradius(B_arr, "e-", Vperp=Vperp1, T_i=T_i)
         assert_quantity_allclose(Vperp1, Vperp2)
 
 
 def test_plasma_frequency():
     r"""Test the plasma_frequency function in parameters.py."""
 
-    assert plasma_frequency(n_e).unit.is_equivalent(u.rad / u.s)
+    assert plasma_frequency(n_e, "e-").unit.is_equivalent(u.rad / u.s)
 
-    assert plasma_frequency(n_e, to_hz=True).unit.is_equivalent(u.Hz)
+    assert plasma_frequency(n_e, "e-", to_hz=True).unit.is_equivalent(u.Hz)
 
-    assert np.isclose(plasma_frequency(1 * u.cm ** -3).value, 5.64e4, rtol=1e-2)
+    assert np.isclose(plasma_frequency(1 * u.cm ** -3, "e-").value, 5.64e4, rtol=1e-2)
 
     assert np.isclose(
         plasma_frequency(1 * u.cm ** -3, particle="N").value, 3.53e2, rtol=1e-1
@@ -804,15 +832,15 @@ def test_plasma_frequency():
     )
 
     with pytest.raises(TypeError):
-        plasma_frequency(u.m ** -3)
+        plasma_frequency(u.m ** -3, "e-")
 
     with pytest.raises(u.UnitTypeError):
-        plasma_frequency(5 * u.m ** -2)
+        plasma_frequency(5 * u.m ** -2, "e-")
 
-    assert np.isnan(plasma_frequency(np.nan * u.m ** -3))
+    assert np.isnan(plasma_frequency(np.nan * u.m ** -3, "e-"))
 
     with pytest.warns(u.UnitsWarning):
-        assert plasma_frequency(1e19) == plasma_frequency(1e19 * u.m ** -3)
+        assert plasma_frequency(1e19, "e-") == plasma_frequency(1e19 * u.m ** -3, "e-")
 
         assert plasma_frequency(n_i, particle="p").unit.is_equivalent(u.rad / u.s)
 
@@ -1046,8 +1074,8 @@ def test_upper_hybrid_frequency():
 
     omega_uh = upper_hybrid_frequency(B, n_e=n_e)
     omega_uh_hz = upper_hybrid_frequency(B, n_e=n_e, to_hz=True)
-    omega_ce = gyrofrequency(B)
-    omega_pe = plasma_frequency(n=n_e)
+    omega_ce = gyrofrequency(B, "e-")
+    omega_pe = plasma_frequency(n=n_e, particle="e-")
     assert omega_ce.unit.is_equivalent(u.rad / u.s)
     assert omega_pe.unit.is_equivalent(u.rad / u.s)
     assert omega_uh.unit.is_equivalent(u.rad / u.s)
@@ -1080,7 +1108,7 @@ def test_lower_hybrid_frequency():
     ion = "He-4 1+"
     omega_ci = gyrofrequency(B, particle=ion)
     omega_pi = plasma_frequency(n=n_i, particle=ion)
-    omega_ce = gyrofrequency(B)
+    omega_ce = gyrofrequency(B, "e-")
     omega_lh = lower_hybrid_frequency(B, n_i=n_i, ion=ion)
     omega_lh_hz = lower_hybrid_frequency(B, n_i=n_i, ion=ion, to_hz=True)
     assert omega_ci.unit.is_equivalent(u.rad / u.s)
@@ -1105,8 +1133,8 @@ def test_lower_hybrid_frequency():
         lower_hybrid_frequency(np.nan * u.T, n_i=-5e19 * u.m ** -3, ion="asdfasd")
 
     with pytest.warns(u.UnitsWarning):
-        assert lower_hybrid_frequency(1.3, 1e19) == lower_hybrid_frequency(
-            1.3 * u.T, 1e19 * u.m ** -3
+        assert lower_hybrid_frequency(1.3, 1e19, "p+") == lower_hybrid_frequency(
+            1.3 * u.T, 1e19 * u.m ** -3, "p+"
         )
     assert_can_handle_nparray(lower_hybrid_frequency)