Repository: PlasmaPy/PlasmaPy

Issue:
1D Maxwellian example update
The 1D Maxwellian [example](http://docs.plasmapy.org/en/stable/auto_examples/plot_distribution.html#sphx-glr-auto-examples-plot-distribution-py) in our docs should have a plot for 1D speed distribution with lines that demarcate the various thermal velocities (i.e. "most probable", "mean magnitude", etc.).

Discussion:
pep8speaks: Hello @pheuer! Thanks for submitting your pull request.

* In the file [`plasmapy/formulary/parameters.py`](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/parameters.py):

> [Line 425:9](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/parameters.py#L425): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 425:11](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/parameters.py#L425): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals

* In the file [`plasmapy/formulary/tests/test_parameters.py`](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/tests/test_parameters.py):

> [Line 1097:1](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/tests/test_parameters.py#L1097): [E305](https://duckduckgo.com/?q=pep8%20E305) expected 2 blank lines after class or function definition, found 1
> [Line 1097:21](https://github.com/PlasmaPy/PlasmaPy/blob/7ef6c4e31b091e2eb8cf7477687ed8d35ab42bf7/plasmapy/formulary/tests/test_parameters.py#L1097): [W292](https://duckduckgo.com/?q=pep8%20W292) no newline at end of file

PlasmaPy follows the [PEP8 style guide](https://www.python.org/dev/peps/pep-0008/), which may be checked locally using [pycodestyle](http://pycodestyle.pycqa.org/en/latest/intro.html)
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850?src=pr&el=h1) Report
> Merging [#850](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/8abf3d089c3046383651d0fd9eef4dc96e9f97b7&el=desc) will **increase** coverage by `0.01%`.
> The diff coverage is `87.50%`.

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

```diff
@@            Coverage Diff             @@
##           master     #850      +/-   ##
==========================================
+ Coverage   95.88%   95.89%   +0.01%     
==========================================
  Files          56       56              
  Lines        5127     5144      +17     
==========================================
+ Hits         4916     4933      +17     
  Misses        211      211              
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/parameters.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3BhcmFtZXRlcnMucHk=) | `97.46% <87.50%> (-1.43%)` | :arrow_down: |
| [plasmapy/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvX19pbml0X18ucHk=) | `32.14% <0.00%> (-2.34%)` | :arrow_down: |
| [plasmapy/utils/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvdXRpbHMvX19pbml0X18ucHk=) | `100.00% <0.00%> (ø)` | |
| [plasmapy/classes/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvY2xhc3Nlcy9fX2luaXRfXy5weQ==) | `100.00% <0.00%> (ø)` | |
| [plasmapy/particles/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL2F0b21pYy5weQ==) | `100.00% <0.00%> (ø)` | |
| [plasmapy/utils/exceptions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvdXRpbHMvZXhjZXB0aW9ucy5weQ==) | `100.00% <0.00%> (ø)` | |
| [plasmapy/formulary/quantum.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3F1YW50dW0ucHk=) | `65.85% <0.00%> (ø)` | |
| [plasmapy/particles/nuclear.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL251Y2xlYXIucHk=) | `100.00% <0.00%> (ø)` | |
| [plasmapy/particles/symbols.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL3N5bWJvbHMucHk=) | `100.00% <0.00%> (ø)` | |
| [plasmapy/utils/roman/roman.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree#diff-cGxhc21hcHkvdXRpbHMvcm9tYW4vcm9tYW4ucHk=) | `100.00% <0.00%> (ø)` | |
| ... and [29 more](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/850/diff?src=pr&el=tree-more) | |

------

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

pheuer: @StanczakDominik Thanks for the suggestions! I've implemented the dictionary refactor: that's very elegant. A few other changes: 

- I ran pre-commit and it worked well (once I updated the fork to the current master: before that it wanted to fix lots of files). 

- I updated the example notebook to have a thumbnail. I ran into the same problem here as I did with the Thomson notebook where multi-figure plots are required but make bad thumbnails, so a single plot must also be included. Maybe theres a better way to handle this through a more sophisticated use of nbsphinx? 

StanczakDominik: > I updated the example notebook to have a thumbnail. I ran into the same problem here as I did with the Thomson notebook where multi-figure plots are required but make bad thumbnails, so a single plot must also be included. Maybe theres a better way to handle this through a more sophisticated use of nbsphinx?

I don't think there is, or at least I'm not aware of it yet :(

 I'll give this a quick once-over again and merge it, most likely :)
pheuer: @StanczakDominik  Looks good - most of the tests I've written for both this and the Thomson module have been known value comparisons: I can see the potential problems with this approach but I can't think of a better way to ensure the physics equations haven't been modified?
StanczakDominik: > @StanczakDominik  Looks good - most of the tests I've written for both this and the Thomson module have been known value comparisons: I can see the potential problems with this approach but I can't think of a better way to ensure the physics equations haven't been modified?

Me neither, unless the original derivation provided value samples... Which they usually don't. Oh well.
StanczakDominik: We could potentially try property based testing one of these days... I'd have to pick it up from scratch though.

No matter! Merging :)

Patch:
diff --git a/changelog/850.feature.rst b/changelog/850.feature.rst
new file mode 100644
index 00000000..23bd21f6
--- /dev/null
+++ b/changelog/850.feature.rst
@@ -0,0 +1,1 @@
+Revised parameters.thermal_speed to support 1D and 2D distributions as well as 3D, and added an example notebook for this function.
diff --git a/docs/formulary/distribution.rst b/docs/formulary/distribution.rst
index c805261e..2cc487b7 100644
--- a/docs/formulary/distribution.rst
+++ b/docs/formulary/distribution.rst
@@ -13,3 +13,4 @@ Distribution functions (`plasmapy.formulary.distribution`)
     :caption: Examples
 
     /notebooks/distribution
+    /notebooks/thermal_speed
diff --git a/docs/formulary/parameters.rst b/docs/formulary/parameters.rst
index 76bd8de4..65ba803c 100644
--- a/docs/formulary/parameters.rst
+++ b/docs/formulary/parameters.rst
@@ -11,3 +11,4 @@ Plasma parameters (`plasmapy.formulary.parameters`)
     :caption: Examples
 
     /notebooks/distribution
+    /notebooks/thermal_speed
diff --git a/docs/notebooks/thermal_speed.ipynb b/docs/notebooks/thermal_speed.ipynb
new file mode 100644
index 00000000..4dc999d9
--- /dev/null
+++ b/docs/notebooks/thermal_speed.ipynb
@@ -0,0 +1,156 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Thermal Speed"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import numpy as np\n",
+    "from astropy import units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from plasmapy.formulary import  Maxwellian_speed_1D, Maxwellian_speed_2D, Maxwellian_speed_3D\n",
+    "from plasmapy.formulary.parameters import thermal_speed"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The thermal_speed function can be used to calculate the thermal velocity for a Maxwellian velocity distribution. There are three common definitions of the thermal velocity, which can be selected using the \"method\" keyword, which are defined for a 3D velocity distribution as\n",
+    "\n",
+    "\n",
+    "- 'most_probable' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{2 k_B T}{m}}$\n",
+    "\n",
+    "- 'rms' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{3 k_B T}{m}}$\n",
+    "\n",
+    "- 'mean_magnitude' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{8 k_B T}{m\\pi}}$\n",
+    "\n",
+    "The differences between these velocities can be seen by plotitng them on a 3D Maxwellian speed distribution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+    "tooltip": "Thermal Speeds"
+    }
+    },
+   "outputs": [],
+   "source": [
+    "T = 1e5 * u.K\n",
+    "speeds = np.linspace(0, 8e6, num=600) * u.m/u.s\n",
+    "\n",
+    "pdf_3D = Maxwellian_speed_3D(speeds, T=T, particle='e-')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(4,3))\n",
+    "\n",
+    "v_most_prob = thermal_speed(T=T, particle='e-', method='most_probable', ndim=3)\n",
+    "v_rms = thermal_speed(T=T, particle='e-', method='rms', ndim=3)\n",
+    "v_mean_magnitude = thermal_speed(T=T, particle='e-', method='mean_magnitude', ndim=3)\n",
+    "\n",
+    "ax.plot(speeds/v_rms, pdf_3D, color='black', label='Maxwellian')\n",
+    "    \n",
+    "ax.axvline(x=v_most_prob/v_rms, color='blue', label='Most Probable')\n",
+    "ax.axvline(x=v_rms/v_rms, color='green', label='RMS')\n",
+    "ax.axvline(x=v_mean_magnitude/v_rms, color='red', label='Mean Magnitude')\n",
+    "\n",
+    "ax.set_xlim(-.1, 3)\n",
+    "ax.set_ylim(0, None)\n",
+    "ax.set_title('3D')\n",
+    "ax.set_xlabel(\"|v|/|v$_{rms}|$\")\n",
+    "ax.set_ylabel(\"f(|v|)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similar speeds are defined for 1D and 2D distributions. The differences between these definitions can be illustrated by plotting them on their respective Maxwellian speed distributions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "pdf_1D = Maxwellian_speed_1D(speeds, T=T, particle='e-')\n",
+    "pdf_2D = Maxwellian_speed_2D(speeds, T=T, particle='e-')\n",
+    "\n",
+    "dim = [1,2,3]\n",
+    "pdfs = [pdf_1D, pdf_2D, pdf_3D]\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "fig, ax = plt.subplots(ncols=3, figsize=(10,3))\n",
+    "\n",
+    "for n, pdf in enumerate(pdfs):\n",
+    "    ndim = n+1\n",
+    "    v_most_prob = thermal_speed(T=T, particle='e-', method='most_probable', ndim=ndim)\n",
+    "    v_rms = thermal_speed(T=T, particle='e-', method='rms', ndim=ndim)\n",
+    "    v_mean_magnitude = thermal_speed(T=T, particle='e-', method='mean_magnitude', ndim=ndim)\n",
+    "    \n",
+    "    ax[n].plot(speeds/v_rms, pdf, color='black', label='Maxwellian')\n",
+    "    \n",
+    "    ax[n].axvline(x=v_most_prob/v_rms, color='blue', label='Most Probable')\n",
+    "    ax[n].axvline(x=v_rms/v_rms, color='green', label='RMS')\n",
+    "    ax[n].axvline(x=v_mean_magnitude/v_rms, color='red', label='Mean Magnitude')\n",
+    "    \n",
+    "    ax[n].set_xlim(-.1, 3)\n",
+    "    ax[n].set_ylim(0, None)\n",
+    "    ax[n].set_title('{:d}D'.format(ndim))\n",
+    "    ax[n].set_xlabel(\"|v|/|v$_{rms}|$\")\n",
+    "    ax[n].set_ylabel(\"f(|v|)\")\n",
+    "\n",
+    "\n",
+    "ax[2].legend(bbox_to_anchor=(1.9, .8), loc='upper right')\n",
+    "\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/plasmapy/formulary/parameters.py b/plasmapy/formulary/parameters.py
index 5068e9b9..e1db4115 100644
--- a/plasmapy/formulary/parameters.py
+++ b/plasmapy/formulary/parameters.py
@@ -412,6 +412,16 @@ def ion_sound_speed(
     return V_S
 
 
+# This dictionary defines coefficients for thermal speeds
+# calculated for different methods and values of ndim.
+# Created here to avoid re-instantiating on each call
+_coefficients = {
+    1: {"most_probable": 0, "rms": 1, "mean_magnitude": 2 / np.pi},
+    2: {"most_probable": 1, "rms": 2, "mean_magnitude": np.pi / 2},
+    3: {"most_probable": 2, "rms": 3, "mean_magnitude": 8 / np.pi},
+}
+
+
 @check_relativistic
 @validate_quantities(
     T={"can_be_negative": False, "equivalencies": u.temperature_energy()},
@@ -423,6 +433,7 @@ def thermal_speed(
     particle: particles.Particle = "e-",
     method="most_probable",
     mass: u.kg = np.nan * u.kg,
+    ndim=3,
 ) -> u.m / u.s:
     r"""
     Return the most probable speed for a particle within a Maxwellian
@@ -448,6 +459,10 @@ def thermal_speed(
         anything, but if set, overrides mass acquired from `particle`. Useful
         with relative velocities of particles.
 
+    ndim : int
+        Dimensionality of space in which to calculate thermal velocity. Valid
+        values are 1,2,3.
+
     Returns
     -------
     V : ~astropy.units.Quantity
@@ -479,10 +494,34 @@ def thermal_speed(
     The particle thermal speed is given by:
 
     .. math::
-        V_{th,i} = \sqrt{\frac{2 k_B T_i}{m_i}}
-
-    This function yields the most probable speed within a distribution
-    function.  However, the definition of thermal velocity varies by
+        V_{th,i} = \sqrt{\frac{N k_B T_i}{m_i}}
+
+    where the value of N depends on the dimensionality and the definition of
+    :math:`v_{th}`: most probable, root-mean-square (RMS), or mean magnitude.
+    The value of N in each case is
+
+    .. list-table:: Values of constant N
+       :widths: 50, 25, 25, 25
+       :header-rows: 1
+
+       * - Dim.
+         - Most-Probable
+         - RMS
+         - Mean-Magnitude
+       * - 1D
+         - 0
+         - 1
+         - :math:`2/\pi`
+       * - 2D
+         - 1
+         - 2
+         - :math:`\pi/2`
+       * - 3D
+         - 2
+         - 3
+         - :math:`8/\pi`
+
+    The definition of thermal velocity varies by
     the square root of two depending on whether or not this velocity
     absorbs that factor in the expression for a Maxwellian
     distribution.  In particular, the expression given in the NRL
@@ -509,16 +548,16 @@ def thermal_speed(
     m = mass if np.isfinite(mass) else particles.particle_mass(particle)
 
     # different methods, as per https://en.wikipedia.org/wiki/Thermal_velocity
-    if method == "most_probable":
-        V = np.sqrt(2 * k_B * T / m)
-    elif method == "rms":
-        V = np.sqrt(3 * k_B * T / m)
-    elif method == "mean_magnitude":
-        V = np.sqrt(8 * k_B * T / (m * np.pi))
-    else:
+    try:
+        coef = _coefficients[ndim]
+    except KeyError:
+        raise ValueError("{ndim} is not a supported value for ndim in " "thermal_speed")
+    try:
+        coef = coef[method]
+    except KeyError:
         raise ValueError("Method {method} not supported in thermal_speed")
 
-    return V
+    return np.sqrt(coef * k_B * T / m)
 
 
 @validate_quantities(


Test patch:
diff --git a/plasmapy/formulary/tests/test_parameters.py b/plasmapy/formulary/tests/test_parameters.py
index 7fc8486b..bc5435f8 100644
--- a/plasmapy/formulary/tests/test_parameters.py
+++ b/plasmapy/formulary/tests/test_parameters.py
@@ -365,6 +365,45 @@ def test_thermal_speed():
         thermal_speed(1 * u.MK, particle="p").si.value, 128486.56960876315
     )
 
+    # 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)
+
+    # Regression tests start here!
+    assert np.isclose(
+        thermal_speed(T_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
+    )
+
+    # ndim = 2
+    assert np.isclose(
+        thermal_speed(T_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
+    )
+
+    assert np.isclose(
+        thermal_speed(T_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
+    )
+
+    assert np.isclose(
+        thermal_speed(T_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
+    )
+
     # Case when Z=1 is assumed
     assert thermal_speed(T_i, particle="p") == thermal_speed(T_i, particle="H-1+")
 
@@ -398,9 +437,14 @@ def test_thermal_speed():
         thermal_speed(1e6 * u.K, method="rms").si.value, 6743070.475775486
     )
 
+    # Test invalid method
     with pytest.raises(ValueError):
         thermal_speed(T_i, method="sadks")
 
+    # Test invalid ndim
+    with pytest.raises(ValueError):
+        thermal_speed(T_i, ndim=4)
+
     assert_can_handle_nparray(thermal_speed)