Repository: PlasmaPy/PlasmaPy
Pull request #831 (2020-06-03 22:32:26)

Thomson Scattering Module Initial Commit
This PR adds a file plasmapy/diagnostics/thomson.py which contain a code for simulating Thomson scattering spectra for Maxwellian plasmas in both the collective and non-collective regimes.

An associated file of tests is also added: plasmapy/diagnostics/tests/test_thomson.py. This file contains two tests that run the program in two regimes (collective and non-collective) and make sure the results match their current output values. 










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

* In the file [`plasmapy/diagnostics/tests/test_thomson.py`](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py):

> [Line 28:26](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L28): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 28:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L28): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 29:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L29): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 29:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L29): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 34:13](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L34): [E222](https://duckduckgo.com/?q=pep8%20E222) multiple spaces after operator
> [Line 35:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L35): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 35:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L35): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 36:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L36): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 36:29](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L36): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 39:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L39): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 39:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L39): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 39:33](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L39): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 39:35](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L39): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 40:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L40): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 40:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L40): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 40:29](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L40): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 41:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L41): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 47:1](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L47): [E303](https://duckduckgo.com/?q=pep8%20E303) too many blank lines (3)
> [Line 55:26](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L55): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 55:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L55): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 56:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L56): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 56:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L56): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 61:13](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L61): [E222](https://duckduckgo.com/?q=pep8%20E222) multiple spaces after operator
> [Line 62:28](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L62): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 62:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L62): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 63:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L63): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 63:29](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L63): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 66:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L66): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 66:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L66): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 66:33](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L66): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 66:35](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L66): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 67:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L67): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 67:27](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L67): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 67:29](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L67): [E251](https://duckduckgo.com/?q=pep8%20E251) unexpected spaces around keyword / parameter equals
> [Line 68:22](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L68): [E128](https://duckduckgo.com/?q=pep8%20E128) continuation line under-indented for visual indent
> [Line 79:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L79): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 80:100](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L80): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (111 > 99 characters)
> [Line 85:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L85): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 86:100](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L86): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (118 > 99 characters)
> [Line 87:100](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L87): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (125 > 99 characters)
> [Line 96:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L96): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 97:100](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L97): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (119 > 99 characters)
> [Line 101:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L101): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 102:100](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/tests/test_thomson.py#L102): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (129 > 99 characters)

* In the file [`plasmapy/diagnostics/thomson.py`](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py):

> [Line 17:1](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L17): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 18:1](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L18): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 19:1](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L19): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 20:1](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L20): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 105:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L105): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 109:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L109): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 118:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L118): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 119:6](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L119): [E225](https://duckduckgo.com/?q=pep8%20E225) missing whitespace around operator
> [Line 119:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L119): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 119:31](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L119): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 120:38](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L120): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 120:39](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L120): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 121:19](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L121): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 121:20](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L121): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 122:19](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L122): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 122:20](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L122): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 123:25](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L123): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 123:26](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L123): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 124:25](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L124): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 124:26](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L124): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 127:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L127): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 130:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L130): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 131:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L131): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 135:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L135): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 138:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L138): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 139:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L139): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 143:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L143): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 144:47](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L144): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 144:48](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L144): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 145:52](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L145): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 145:53](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L145): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 146:30](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L146): [E261](https://duckduckgo.com/?q=pep8%20E261) at least two spaces before inline comment
> [Line 146:31](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L146): [E262](https://duckduckgo.com/?q=pep8%20E262) inline comment should start with '# '
> [Line 148:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L148): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 149:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L149): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 150:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L150): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 154:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L154): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 155:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L155): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 158:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L158): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 162:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L162): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 163:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L163): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 164:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L164): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 171:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L171): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 174:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L174): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '
> [Line 175:36](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L175): [E225](https://duckduckgo.com/?q=pep8%20E225) missing whitespace around operator
> [Line 180:58](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L180): [E225](https://duckduckgo.com/?q=pep8%20E225) missing whitespace around operator
> [Line 183:5](https://github.com/PlasmaPy/PlasmaPy/blob/92c38c7edde8bd8b75ba5dc1644f3797fab5a07d/plasmapy/diagnostics/thomson.py#L183): [E265](https://duckduckgo.com/?q=pep8%20E265) block comment should start with '# '

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/831?src=pr&el=h1) Report
> Merging [#831](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/831?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/eb463877194207b7e1fde025a52fef609a21ed57&el=desc) will **increase** coverage by `0.00%`.
> The diff coverage is `96.00%`.

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

```diff
@@           Coverage Diff           @@
##           master     #831   +/-   ##
=======================================
  Coverage   95.88%   95.88%           
=======================================
  Files          56       57    +1     
  Lines        5123     5173   +50     
=======================================
+ Hits         4912     4960   +48     
- Misses        211      213    +2     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/831?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/diagnostics/thomson.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/831/diff?src=pr&el=tree#diff-cGxhc21hcHkvZGlhZ25vc3RpY3MvdGhvbXNvbi5weQ==) | `96.00% <96.00%> (ø)` | |

------

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

rocco8773: Hi @pheuer, thanks again for doing this.  Like I said in our emails, I'll try to give a proper review by the end of the week.  To get started, look into...

1. Replacing the use of `quantity_input` with [`validate_quantities`](http://docs.plasmapy.org/en/latest/api/plasmapy.utils.decorators.validate_quantities.html#plasmapy.utils.decorators.validate_quantities).
2. The [documentation page for diagnostics](https://github.com/pheuer/PlasmaPy/blob/master/docs/diagnostics/index.rst) will need to be updated.
3. It might be useful to create a jupyter notebook to show how the package works.  This would go in `./docs/notebooks` and we could add it to the `nbgallery` you see on the documentation page.
StanczakDominik: I'd definitely like to also take a look at this functionality :)
pheuer: Thanks for the comments! In addition to the changes noted above, I have also created an example Jupyter notebook and *attempted* to edit the diagnostics page in the documentation (although I'm not sure if there's more I need to do in the actual function docstrings to make that work properly?).
StanczakDominik: I'm also getting some errors with tests on older Astropy versions at https://dev.azure.com/plasmapy/PlasmaPy/_build/results?buildId=993&view=logs&j=8955e93f-e95e-5ac2-4406-a311ed91aa3e&t=adda2fec-11c4-56e4-c47d-199954e67231&l=5683 (seems to point to some `rad/cm` vs `rad/s` inconsistency in the code?), and the documentation build points to warnings on section titles (as mentioned before): https://app.circleci.com/pipelines/github/PlasmaPy/PlasmaPy/287/workflows/044fac81-3940-4f64-8970-dda8018fbd7d/jobs/6586/parallel-runs/0/steps/0-107. I'm not sure Numpydoc has a separate header for keyword-only arguments, but just cleaning up Arguments may help.
StanczakDominik: Oh yeah, if the fix to the units thing turns out to be hard to figure out, I think I'm fine with bumping up the minimal astropy version. They've been pumping those out like crazy.
pheuer: Ok I'm still looking into the astropy.units issue: sifting through the logs the issue seems to be at (currently) line 146: w_e = w - k*np.dot(fluid_vel, k_vec). The units are w->rad/s, k->rad/cm, fluid_vel->cm/s, and k_vec is dimensionless. It looks like, on the older astropy version, the second term is ending up with units of rad/cm and then the error throws when it tries to do the subtraction?
StanczakDominik: > The numpydoc issue should be resolved: There is now only one heading for "parameters" which is their standard. Although, sidebar, I at least find it hard to distinguish between positional and keyword args in this case in the documentation?

Yeah, it is pretty hard to distinguish... @rocco8773, given that you usually have Opinions about documentation layout (and good ones at that!) - any suggestions here?
StanczakDominik: Wrong hotkey, posted too quickly...

> Ok I'm still looking into the astropy.units issue: sifting through the logs the issue seems to be at (currently) line 146: w_e = w - k*np.dot(fluid_vel, k_vec). The units are w->rad/s, k->rad/cm, fluid_vel->cm/s, and k_vec is dimensionless. It looks like, on the older astropy version, the second term is ending up with units of rad/cm and then the error throws when it tries to do the subtraction?

It rang no bells at first, but looking at the formula and seeing `np.dot` reminded me about #59, coming back to bite us in the posterior ~700 issues/PRs later! Given that this was fixed in some later Astropy version, I vote to just bump the minimal req for that. I'm still going to test it to make sure this is it.
StanczakDominik: Hahahaha, well done, good timing :rofl: 4.0 it is, then.
rocco8773: > > The numpydoc issue should be resolved: There is now only one heading for "parameters" which is their standard. Although, sidebar, I at least find it hard to distinguish between positional and keyword args in this case in the documentation?
> 
> Yeah, it is pretty hard to distinguish... @rocco8773, given that you usually have Opinions about documentation layout (and good ones at that!) - any suggestions here?

Yeah, this list is of parameters is fairly long.  If there's some rarely used keywords, then you can add them into an [**Other Parameters**](https://numpydoc.readthedocs.io/en/latest/format.html) section.  Otherwise, I'd stick to just the **Parameters** section for consistency.  To distiguish keywords from argument, I usually explicitly state the default value in their definition.

@pheuer To help yourself distinguish arguments vs keywords in the future, look at the function signature.  Arguments will always be first in line and keywords will always have default values.
StanczakDominik: > > > The numpydoc issue should be resolved: There is now only one heading for "parameters" which is their standard. Although, sidebar, I at least find it hard to distinguish between positional and keyword args in this case in the documentation?
> > 
> > 
> > Yeah, it is pretty hard to distinguish... @rocco8773, given that you usually have Opinions about documentation layout (and good ones at that!) - any suggestions here?
> 
> Yeah, this list is of parameters is fairly long. If there's some rarely used keywords, then you can add them into an [**Other Parameters**](https://numpydoc.readthedocs.io/en/latest/format.html) section. Otherwise, I'd stick to just the **Parameters** section for consistency. To distiguish keywords from argument, I usually explicitly state the default value in their definition.
> 
> @pheuer To help yourself distinguish arguments vs keywords in the future, look at the function signature. Arguments will always be first in line and keywords will always have default values.

To clarify - we're talking about having a single positional argument here (wavelength) and a few being required (non-optional) keyword-only arguments. What would you recommend for that?

Re: astropy bump, I'm making a separate PR for that. Looking at the test results, it'll definitely fix the bug here.
pheuer: Great, thanks for resolving that issue. 

Unfortunately all of the arguments for this function are required (except those with default arguments I suppose). Now that I've corrected the numpydoc formatting the documentation is a little more clear online so I think wavelength is now distinguishable as a positional argument there.
rocco8773: > Hahahaha, well done, good timing 🤣 4.0 it is, then.

My original take was to setup a `try-except` statement to handle to the two version, but the more I think about it I'm good with changing our dependencies to `astropy >= 4.0`.  It's been out for nearly half a year already.
rocco8773: > > > > The numpydoc issue should be resolved: There is now only one heading for "parameters" which is their standard. Although, sidebar, I at least find it hard to distinguish between positional and keyword args in this case in the documentation?
> > > 
> > > 
> > > Yeah, it is pretty hard to distinguish... @rocco8773, given that you usually have Opinions about documentation layout (and good ones at that!) - any suggestions here?
> > 
> > 
> > Yeah, this list is of parameters is fairly long. If there's some rarely used keywords, then you can add them into an [**Other Parameters**](https://numpydoc.readthedocs.io/en/latest/format.html) section. Otherwise, I'd stick to just the **Parameters** section for consistency. To distiguish keywords from argument, I usually explicitly state the default value in their definition.
> > @pheuer To help yourself distinguish arguments vs keywords in the future, look at the function signature. Arguments will always be first in line and keywords will always have default values.
> 
> To clarify - we're talking about having a single positional argument here (wavelength) and a few being required (non-optional) keyword-only arguments. What would you recommend for that?

My take, keyword parameters should NEVER be required.  Their implementation makes them by default optional.  So, something like this is just bad practice...

```python
def foo(x=None):
    if x is None:
        raise ValueError("x can not be a None value")
```
Note though, in the current `spectral_density` function all arguments are positional and keyword, so you can still call the function like `spectral_density(wavelength=wavelength)`.

In summary, if a parameter is required then make it a positional argument (which can still be called like a keyword argument).  If a parameter has default behavior if not given, then make it a keyword argument.
StanczakDominik: https://github.com/PlasmaPy/PlasmaPy/pull/831#discussion_r435447584 I should probably have linked this in the first place, resolved discussions can be hard to find and easy to miss, sorry :( https://www.python.org/dev/peps/pep-3102/#id4 has this sort of thing:

```python
def compare(a, b, *, key):
    ...
```
No default value on `key`, but it can still only be input explicitly with the variable name. This is what I was going for in review, and completely forgot about it later...
rocco8773: > [#831 (comment)](https://github.com/PlasmaPy/PlasmaPy/pull/831#discussion_r435447584) I should probably have linked this in the first place, resolved discussions can be hard to find and easy to miss, sorry :( https://www.python.org/dev/peps/pep-3102/#id4 has this sort of thing:
> 
> ```python
> def compare(a, b, *, key):
>     ...
> ```
> 
> No default value on `key`, but it can still only be input explicitly with the variable name. This is what I was going for in review, and completely forgot about it later...

?? Do you mean?

```python
def compare(a, b, *, key=None):
    ...
```

This syntax just enforces that everything after `*` is a keyword-only argument. ([types of paramters](https://docs.python.org/3/glossary.html#term-parameter))
StanczakDominik: Nope, actually, it gets cooler than that!
```python
[ins] In [1]: def compare(a, b, *, key): 
         ...:     print(key) 
         ...:     return a ==b 
         ...:                                                                                          

[ins] In [2]: compare(1, 2)                                                                            
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-5c7da606c74f> in <module>
----> 1 compare(1, 2)

TypeError: compare() missing 1 required keyword-only argument: 'key'

[ins] In [3]: compare(1, 2, 3)                                                                         
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-b0c1df51f863> in <module>
----> 1 compare(1, 2, 3)

TypeError: compare() takes 2 positional arguments but 3 were given

[ins] In [4]: compare(1, 2, key = "I AM A POTATO")                                                     
I AM A POTATO
Out[4]: False

```
rocco8773: Oh weird! I've never seen that before.  Not sure how I feel about it, but I'm warming up to it.  I'm surprised Python does not showcase that better.  This is the only line from [PEP 3102](https://www.python.org/dev/peps/pep-3102/#specification) I've ever seen on this...

> Keyword-only arguments are not required to have a default value. Since Python requires that all arguments be bound to a value, and since the only way to bind a value to a keyword-only argument is via keyword, such arguments are therefore 'required keyword' arguments. Such arguments must be supplied by the caller, and they must be supplied via keyword.
rocco8773: @pheuer How do you typically use this in your workflow?  Do you manually generate it to overlay onto collected data or do you try to fit it?

I'm trying to decide if we should be a classed based approach, kinda like this [comment](https://github.com/PlasmaPy/PlasmaPy/pull/816#discussion_r432699884).  If there is ever a need for fitting or integrating, then a unit-less access could be useful while maintaining the common API as a unit based one.
pheuer: @rocco8773 

*Disclaimer: I've never actually processed Thomson scattering data in the wild, I've been using this to help plan out a new LAPD diagnostic* 

I think the use-case is 50/50 predictive and fitting. I have done some experimenting fitting synthetic data (using curve_fit and a lambda function to isolate just a few fitting parameters), so that would be good functionality to support so I guess  a unit-less function would be worth having. 

I'm not sure I understand the advantage of the class structure you linked to - why not just include a wrapper  function that takes float inputs and assumes some units when it calls the main function (or vice-versa)? 


StanczakDominik: Maybe let's keep it simple for now? You could apply that logic to most of `formulary`. It probably doesn't need to be that hyper optimized.

Besides, if you're integrating this density, seems like you can quickly calculate it first, then dump units and integrate numerically. It's only a problem for automatic integration routines, and I'm sure we'll figure that out sooner or later.

Off to bed for me! Sorry if I seem grumpy.
lemmatum: Hey, just chiming in. You can greatly simplify the code in `thomson.py` by using `permittivity_1D_Maxwellian` from `dielectric.py`. This calculates your chi_e and chi_i for you. Thomson scattering simulations is one of the reasons I wrote that code up... 

Also, I've got a little Thomson scattering code I've been working on. It might help add functionality to this, whenever this merge is done. Only issue with it is I've got a bug somewhere which is breaking detailed balance (the upshifted and downshifted plasmon peaks don't look right).

I agree with @pheuer that this module can be used for both prediction and fitting. I mostly see people doing Thomson scattering fits "by hand", meaning they vary the parameters and eyeball the best fit. Sometimes they do a chi-squared. Ideally, they should be doing something that can handle multi-dimensional fits and give error bars, like Markov Chain Monte Carlo. It might be too much to ask for fitting right now... having the code simulate spectra is already quite powerful and users can build their own fitting routines around it as necessary.

Last point: I see you have an IPython notebook. I think you should translate that into an example for the gallery.

Let me know if you guys need any more help with this. I've done Thomson scattering experiments and analysis before (mainly x-ray).
rocco8773: Sleep tight Mr Grumpy. 🤣

If you're going to have numerous functions and data that are highly dependent on each other, then classes do a better job of encapsulating everything (both for functionality and readability).  I think that's the direction this would go in the future, but, as @StanczakDominik said, it does not need to happen for this PR.

Additionally, I'm currently prototyping the Diagnostic Framework, which will hopefully come online in a few months.  Once that does, this will need to be refactored to merge into the framework.
rocco8773: > I agree with @pheuer that this module can be used for both prediction and fitting. I mostly see people doing Thomson scattering fits "by hand", meaning they vary the parameters and eyeball the best fit. Sometimes they do a chi-squared. Ideally, they should be doing something that can handle multi-dimensional fits and give error bars, like Markov Chain Monte Carlo. It might be too much to ask for fitting right now... having the code simulate spectra is already quite powerful and users can build their own fitting routines around it as necessary.

Yeah, this is what I was trying to get at.  I think functionality like this would be very useful to the community, but does not need to be implemented now.
rocco8773: > Last point: I see you have an IPython notebook. I think you should translate that into an example for the gallery.

Yes, it is currently setup like that.

https://plasmapy--831.org.readthedocs.build/en/831/examples.html
https://plasmapy--831.org.readthedocs.build/en/831/diagnostics/thomson.html
lemmatum: > > Last point: I see you have an IPython notebook. I think you should translate that into an example for the gallery.
> 
> Yes, it is currently setup like that.
> 
> https://plasmapy--831.org.readthedocs.build/en/831/examples.html
> https://plasmapy--831.org.readthedocs.build/en/831/diagnostics/thomson.html

Ah sweet! I missed that!
lemmatum: Looking through that example in the gallery, it looks like detailed balance isn't implemented in this. The upshifted and downshifted peaks are symmetric (same shape and signal level) when in reality the downshifted peak is stronger than the upshifted peak. We don't need to add this now, but it definitely needs to be added in the future.
pheuer: @lemmatum Thanks for pointing out the dielectric functions! They didn't shorten the code much (because I still need to compute the phase velocities), but it's better to not re-invent the wheel anyway.  I also think I found a sqrt(2) error in my previous code while converting over to your functions.

While I was changing that, I realized it was more convenient to use the plasmapy.Particle class to allow users to input ion species, so I've changed over to doing that.

Regarding the height of the peaks: my understanding is that this arises from the geometric factor G(k,w) ( [https://aip.scitation.org/doi/10.1063/1.3478975](as described here)). I didn't include this in S(k,w) because it seems to be written separately in most of theory and, in the regimes I'm used to working in, it's effectively unity.  However I know for x-ray TS this factor becomes very important (as v/c gets larger), so we should find a way to deal with it. I could include the first-order correction easily, but I think for x-ray spectra you really need a fully relativistic treatment of everything?

I hadn't heard the connection to detailed balance before: if you have a reference on that I'd be interested! 


pheuer: @rocco8773 Thanks for the suggestions: I've accepted the ones i could and manually edited the __init__.py file. I did add the metadata to the notebook, but I don't see a thumbnail yet.
rocco8773: > @rocco8773 Thanks for the suggestions: I've accepted the ones i could and manually edited the **init**.py file. I did add the metadata to the notebook, but I don't see a thumbnail yet.

Yeah, I was playing around with that myself.  It looks like the output is being treated as one image, versus separate images for each sub-figure.  You can either remove the `"output-index"` metadata or re-work the plots so each sub-figure is its own image.  The former sounds way easier at the moment.
lemmatum: @pheuer Oh yea, the phase velocities and all that stuff does take up a bunch of space... Now that I think about it, perhaps those should be separated out into their own functions instead of hiding in them in one monolithic function? They may be useful on their own. Same goes for the scatter vectors.

Yea separating it into a geometrical factor or something else makes sense. I just meant that we should implement detailed balance somewhere, and not that it necessarily has to go into the DSF. The Ross paper shows optical Thomson scattering spectra with asymmetric peaks and explicitly mentions that the relativistic corrections are important in this regime, even for low temperatures, so this is not just an XRTS term.

Equation 11.2.12 in the Froula, Glenzer, Luhmann, Sheffield book shows the detailed balance term along with a couple references particular to XRTS. The gist of it is that in LTE you have detailed balance and that balance means that the probe beam is more likely to deposit energy into the plasma than the plasma organizing itself in such a way that it donates energy to the scattered probe beam. This can break when you don't have LTE (I wrote a paper on that if you're interested).

The geometrical factor presented in Ross looks a little different from the detailed balance factor for XRTS, so whenever we implement this we should give the user options :smile: 
StanczakDominik: Nice to see you back here, @lemmatum :tada: 

>  I mostly see people doing Thomson scattering fits "by hand", meaning they vary the parameters and eyeball the best fit. Sometimes they do a chi-squared. Ideally, they should be doing something that can handle multi-dimensional fits and give error bars, like Markov Chain Monte Carlo. It might be too much to ask for fitting right now..

That's kind of sort of what my masters' thesis is turning out to be. [PyMC3](https://docs.pymc.io/) is a neat library! I should be able to contribute this back here soon-ish, and for now.
rocco8773: I think the notebook is a good start, but it needs to me more than just informational.  It also needs to be educational. Can you updated it with some images and more physics?
lemmatum: 
> That's kind of sort of what my masters' thesis is turning out to be. [PyMC3](https://docs.pymc.io/) is a neat library! I should be able to contribute this back here soon-ish.

Sweet! Teach me your ways! :mage_man: 

pheuer: @rocco8773 

Thanks for the suggestions, here are a few questions/comments:

1) validate_quantities throws an ImplicitUnitWarningConversion warning if you use units if eV: is there a way to silence this? Entering temperatures in eV is almost definitely going to be the norm, so it should not print warnings. 

2) I've added some more details to the docstring, including an equation that I hope helps make clear what S(k,w) is. The full derivation of this expression is the subject of several chapters in the Sheffield book, so ultimately users will have to look there for a fundamental understanding. 

3) I have added some more explanation to the notebook: I cannot currently think of any other plots that I think would be informative. Eventually it would be nice to have other functions in this module that treat other cases (magnetized plasmas, relativistic effects, calculating power spectra) which would then invite a more detailed explanation of how all these parts fit together. 
pheuer: @lemmatum Thank you for your explanation: that makes perfect sense. I will check out your paper too. 

I think the G(k,w) in Ross is a 1st order limit in v/c, but I gather that a more complete treatment is necessary for XRTS. Assuming the full version is not too computationally taxing and reduces correctly in the non-relativistic limit, I suppose that would be the thing to implement, but I have a lot more to learn before I could start such a project...

Ideally I think this module would have different functions for each term in the scattered power spectrum (as well as a wrapper function that uses them all to produce a predicted spectrum).
rocco8773: > 1. validate_quantities throws an ImplicitUnitWarningConversion warning if you use units if eV: is there a way to silence this? Entering temperatures in eV is almost definitely going to be the norm, so it should not print warnings.

The is put there on purpose to indicate that a non-standard unit conversion is occurring.  With the way warnings work, that will only appear for the first occurrence in a python session.  It's up to the user to filter that out if he or she doesn't want it. If you don't want it to occur in your own scripts then you can add something like `warnings.simplefilter("ignore", category=ImplicitUnitConversionWarning)` to the top of your script. (see [Python's warning suppression](https://docs.python.org/3/library/warnings.html#temporarily-suppressing-warnings)).
rocco8773: > 2. I've added some more details to the docstring, including an equation that I hope helps make clear what S(k,w) is. The full derivation of this expression is the subject of several chapters in the Sheffield book, so ultimately users will have to look there for a fundamental understanding.

Ok, I'll look it over.  I'm not looking for anything detailed, just enough so users understand what's being calculated and a few references so users have somewhere to start searching for more.
StanczakDominik: I've merged #835 in, so that should solve the minimal-conda test failure. Please ignore the flake8 and pre-commit test failures; they're a side effect of work in #845/#846.

I wouldn't want to make this review cycle too long (ekhm transport - @pheuer if you don't understand this little in-joke, be very glad that you don't :laughing:), so I think we should be heading towards merging this, and possibly adding additional fixes/features in follow-ups. I'll give this another look-through now.
pheuer: @StanczakDominik thanks, I updated the notebook as you suggested. I also merged the master branch in, but it doesn't look like it solved the py36-minimal-conda test problem? 

Regarding the module overall, I agree that it would be best to merge this function into master before contemplating further functionality. There's a whole mess of Thomson-related functionality that could be added, but that will take time and should mostly be independent of this function. Maybe we can create a new PR for supporting the relativistic correction/detailed balance that @lemmatum was mentioning after this one is resolved. 
StanczakDominik: > I also merged the master branch in, but it doesn't look like it solved the py36-minimal-conda test problem?

Okay, that's pretty weird. I'll give it another look.
StanczakDominik: I guess the problem was in Numpy, and my previous attempt at bumping its requirements aimed too low.

I also merged `master` again (something seems to have gone wrong and you didn't have the actual current version) and applied `pre-commit run --all-files` to shut the style linter up. It applies the `black` code formatter.
rocco8773: @StanczakDominik Yeah, I still want to do another once over before merging.  I've been busy the last few days, but I'll trying to review by end of Monday.

Patch:
diff --git a/changelog/831.feature.rst b/changelog/831.feature.rst
new file mode 100644
index 00000000..6ae1b3b9
--- /dev/null
+++ b/changelog/831.feature.rst
@@ -0,0 +1,5 @@
+Added a new diagnostics module `thomson` containing a function
+`spectral_density` that calculates Thomson scattering spectra for
+Maxwellian plasmas in both the collective and non-collective regimes. As
+a followup to PR #835, set the minimal required Numpy version to 1.18.1 to
+finally fix unit dropping bugs.
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/about/credits.rst b/docs/about/credits.rst
index 886b091c..604b47f6 100644
--- a/docs/about/credits.rst
+++ b/docs/about/credits.rst
@@ -42,6 +42,7 @@ in parentheses are `ORCID author identifiers <https://orcid.org>`__.
 * `Graham Goudeau <https://github.com/GrahamGoudeau>`__
 * `Silvina Guidoni <https://www.american.edu/cas/faculty/guidoni.cfm>`__
 * `Colby Haggerty <https://github.com/colbych>`__
+* `Peter Heuer <https://github.com/pheuer>`__ (`0000-0001-5050-6606 <https://orcid.org/0000-0001-5050-6606>`__)
 * `Julien Hillairet <https://github.com/jhillairet>`__ (`0000-0002-1073-6383 <https://orcid.org/0000-0002-1073-6383>`__)
 * `Poh Zi How <https://github.com/pohzipohzi>`__
 * `Yi-Min Huang <https://github.com/yopology>`__ (`0000-0002-4237-2211 <https://orcid.org/0000-0002-4237-2211>`__)
diff --git a/docs/diagnostics/index.rst b/docs/diagnostics/index.rst
index 09bbd749..d5d32128 100644
--- a/docs/diagnostics/index.rst
+++ b/docs/diagnostics/index.rst
@@ -4,14 +4,12 @@
 Plasma diagnostics (`plasmapy.diagnostics`)
 *******************************************
 
-.. currentmodule:: plasmapy.diagnostics.langmuir
+.. automodapi:: plasmapy.diagnostics
+    :no-heading:
 
-.. automodapi:: plasmapy.diagnostics.langmuir
-   :no-heading:
+.. toctree::
+   :maxdepth: 1
+   :caption: Sub-Packages and Modules
 
-|
-
-.. nbgallery::
-    :caption: Examples
-
-    /notebooks/langmuir_analysis
+   langmuir
+   thomson
diff --git a/docs/diagnostics/langmuir.rst b/docs/diagnostics/langmuir.rst
new file mode 100644
index 00000000..3d71ec16
--- /dev/null
+++ b/docs/diagnostics/langmuir.rst
@@ -0,0 +1,17 @@
+.. _langmuir:
+
+*****************
+Langmuir analysis
+*****************
+
+.. currentmodule:: plasmapy.diagnostics.langmuir
+
+.. automodapi:: plasmapy.diagnostics.langmuir
+   :no-heading:
+
+|
+
+.. nbgallery::
+    :caption: Examples
+
+    /notebooks/langmuir_analysis
diff --git a/docs/diagnostics/thomson.rst b/docs/diagnostics/thomson.rst
new file mode 100644
index 00000000..96493b3e
--- /dev/null
+++ b/docs/diagnostics/thomson.rst
@@ -0,0 +1,15 @@
+.. _thomson:
+
+******************
+Thomson scattering
+******************
+
+.. currentmodule:: plasmapy.diagnostics.thomson
+
+.. automodapi:: plasmapy.diagnostics.thomson
+   :no-heading:
+
+.. nbgallery::
+    :caption: Examples
+
+    /notebooks/thomson
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/install.rst b/docs/install.rst
index 72b776a9..4aad2653 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -18,7 +18,7 @@ Requirements
 PlasmaPy requires Python version 3.6 or newer.
 PlasmaPy requires the following packages for installation:
 
-- `NumPy <http://www.numpy.org/>`_ 1.16.6 or newer
+- `NumPy <http://www.numpy.org/>`_ 1.18.1 or newer
 - `SciPy <https://www.scipy.org/>`_ 1.2 or newer
 - `Astropy <http://www.astropy.org/>`_ 4.0 or newer
 - `colorama <https://pypi.org/project/colorama/>`_ 0.3 or newer
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/docs/notebooks/thomson.ipynb b/docs/notebooks/thomson.ipynb
new file mode 100644
index 00000000..a26720cb
--- /dev/null
+++ b/docs/notebooks/thomson.ipynb
@@ -0,0 +1,348 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Thomson Scattering: Spectral Density"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[thomson]: ../diagnostics/thomson.rst\n",
+    "[spectral-density]: ../api/plasmapy.diagnostics.thomson.spectral_density.rst#spectral-density\n",
+    "[sheffield]: https://www.sciencedirect.com/book/9780123748775/plasma-scattering-of-electromagnetic-radiation\n",
+    "\n",
+    "The [thomson.spectral_density][spectral-density] function calculates the [spectral density function S(k,w)][sheffield], which is one of several terms that determine the scattered power spectrum for the Thomson scattering of a probe laser beam by a plasma. In particular, this function calculates $S(k,w)$ for the case of a plasma consisting of one or more ion species and a neutralizing electron fluid under the assumption that all of the ion species and the electron fluid have Maxwellian velocity distribution functions. In this regime, the spectral density is given by the equation:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "S(k,\\omega) = \\frac{2\\pi}{k} \\bigg |1 - \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{e0}\\bigg ( \\frac{\\omega}{k} \\bigg ) + \\sum_i \\frac{2\\pi Z_i}{k} \\bigg | \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{i0, i} \\bigg ( \\frac{\\omega}{k} \\bigg )\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\chi_e$ is the electron component susceptibility of the plasma and $\\epsilon = 1 + \\chi_e + \\sum_i \\chi_i$ is the total plasma dielectric function (with $\\chi_i$ being the ion component of the susceptibility), $Z_i$ is the charge of each ion, $k$ is the scattering wavenumber, $\\omega$ is the scattering frequency, and the functions $f_{e0}$ and $f_{i0,i}$ are the Maxwellian velocity distributions for the electrons and ion species respectively.\n",
+    "\n",
+    "Thomson scattering can be either non-collective (the scattered spectrum is a linear sum of the light scattered by individual particles) or collective (the scattered spectrum is dominated by scattering off of collective plasma waves). The [thomson.spectral_density][spectral-density] function can be used in both cases. These regimes are delineated by the dimensionless constant $\\alpha$:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\alpha = \\frac{1}{k \\lambda_{De}}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\lambda_{De}$ is the Debye length. $\\alpha > 1$ corresponds to collective scattering, while $\\alpha < 1$ corresponds to non-collective scattering. Depending on which of these regimes applies, fitting the scattered spectrum can provide the electron (and sometimes ion) density and temperature. Doppler shifting of the spectrum can also provide a measurement of the drift velocity of each plasma species.\n",
+    "\n",
+    "For a detailed explanation of the underlying physics (and derivations of these expressions), see [\"Plasma Scattering of Electromagnetic Radiation\" by Sheffield et al.][sheffield]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import astropy.units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import warnings\n",
+    "\n",
+    "from plasmapy.diagnostics import thomson\n",
+    "from plasmapy.utils.exceptions import ImplicitUnitConversionWarning"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Construct parameters that define the Thomson diagnostic setup, the probing beam and scattering collection.  These parameters will be used for all examples."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The probe wavelength can in theory be anything, but in practice integer frequency multiples of the Nd:YAG wavelength\n",
+    "# 1064 nm are used (532 corresponds to a frequency-doubled probe beam from such a laser).\n",
+    "probe_wavelength = 532*u.nm\n",
+    "\n",
+    "# Array of wavelengths over which to calcualte the spectral distribution\n",
+    "wavelengths = np.arange(probe_wavelength.value-60, probe_wavelength.value+60, 0.01)*u.nm\n",
+    "\n",
+    "# The scattering geometry is defined by unit vectors for the orientation of the probe laser beam (probe_n) and\n",
+    "# the path from the scattering volume (where the measurement is made) to the detector (scatter_n).\n",
+    "# These can be setup for any experimental geometry.\n",
+    "probe_vec = np.array([1, 0, 0])\n",
+    "scattering_angle = np.deg2rad(63)\n",
+    "scatter_vec = np.array([np.cos(scattering_angle), np.sin(scattering_angle), 0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In order to calcluate the scattered spectrum, we must also include some information about the plasma. For this plot we'll allow the ``fract``, ``ion_species``, ``fluid_vel``, and ``ion_vel`` keywords to keep their default values, describing a single-species H+ plasma at rest in the laboratory frame. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+     "tooltip": "Thomson Scattering"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "ne = 2e17*u.cm**-3\n",
+    "Te = 12*u.eV\n",
+    "Ti = 10*u.eV\n",
+    "\n",
+    "# This warning filter catches an ImplicitUnitConversionWarning that results from specifying\n",
+    "# temperatures in eV instead of Kelvin. \n",
+    "with warnings.catch_warnings():\n",
+    "    warnings.simplefilter(\"ignore\", ImplicitUnitConversionWarning)\n",
+    "    \n",
+    "    alpha, Skw = thomson.spectral_density(wavelengths, probe_wavelength,\n",
+    "                                          ne, Te, Ti, probe_vec=probe_vec, \n",
+    "                                          scatter_vec=scatter_vec)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(wavelengths, Skw, lw=2)\n",
+    "ax.set_xlim(probe_wavelength.value-10, probe_wavelength.value+10)\n",
+    "ax.set_ylim(0, 1e-13)\n",
+    "ax.set_xlabel('$\\lambda$ (nm)')\n",
+    "ax.set_ylabel('S(k,w)')\n",
+    "ax.set_title('Thomson Scattering Spectral Density')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Example Cases in Different Scattering Regimes\n",
+    "\n",
+    "We will now consider several example cases in different scattering regimes. In order to facilitate this, we'll set up each example as a dictionary of plasma parameters:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that results in a scattering spectrum dominated by scattering off of single electrons."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "non_collective = {\n",
+    "    'name': 'Non-Collective Regime', \n",
+    "    'ne': 5e15*u.cm**-3, \n",
+    "    'Te': 40*u.eV, \n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]), \n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s, \n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that result in weakly collective scattering (scattering paramter $\\alpha$ approaching 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "weakly_collective = {\n",
+    "    'name': 'Weakly Collective Regime',\n",
+    "    'ne': 2e17*u.cm**-3,\n",
+    "    'Te': 20*u.eV,\n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that result in a spectrum dominated by multi-particle scattering, including scattering off of ions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "collective = {\n",
+    "    'name': 'Collective Regime',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([4])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A case identical to the collective example above, except that now the electron fluid has a substantial drift velocity parallel to the probe laser and the ions have a drift (relative to the electrons) at an angle."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "drifts = {\n",
+    "    'name': 'Drift Velocities',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([700, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[-600, -100, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A case identical to the collective example, except that now the plasma consists 25% He+1 and 75% C+5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "two_species = {\n",
+    "    'name': 'Two Ion Species',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([10, 50])*u.eV,\n",
+    "    'fract': np.array([.25, .75]),\n",
+    "    'ion_species': ['He-4 1+', 'C-12 5+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0],[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "examples = [non_collective, weakly_collective, collective, drifts, two_species]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For each example, plot the the spectral distribution function over a large range to show the broad electron scattering feature (top row) and a narrow range around the probe wavelength to show the ion scattering feature (bottom row)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(ncols=len(examples), nrows=2, figsize=[25,10])\n",
+    "fig.subplots_adjust( wspace=0.4, hspace=0.4)\n",
+    "\n",
+    "lbls = 'abcdefg'\n",
+    "\n",
+    "for i, x in enumerate(examples):\n",
+    "    \n",
+    "    with warnings.catch_warnings():\n",
+    "        warnings.simplefilter(\"ignore\", ImplicitUnitConversionWarning)\n",
+    "        \n",
+    "        alpha, Skw = thomson.spectral_density(wavelengths, probe_wavelength,\n",
+    "                         x['ne'], x['Te'], x['Ti'], fract=x['fract'],\n",
+    "                         ion_species=x['ion_species'], fluid_vel=x['fluid_vel'],\n",
+    "                         probe_vec=probe_vec, scatter_vec=scatter_vec)\n",
+    "    \n",
+    "    ax[0][i].axvline(x=probe_wavelength.value, color='red')  # Mark the probe wavelength\n",
+    "    ax[0][i].plot(wavelengths, Skw)\n",
+    "    ax[0][i].set_xlim(probe_wavelength.value-15, probe_wavelength.value+15)\n",
+    "    ax[0][i].set_ylim(0, 1e-13)\n",
+    "    ax[0][i].set_xlabel('$\\lambda$ (nm)')\n",
+    "    \n",
+    "    ax[0][i].set_title( lbls[i] + ') ' + x['name'] + '\\n$\\\\alpha$={:.4f}'.format(alpha))\n",
+    "\n",
+    "    \n",
+    "    ax[1][i].axvline(x=probe_wavelength.value, color='red')  # Mark the probe wavelength\n",
+    "    ax[1][i].plot(wavelengths, Skw)\n",
+    "    ax[1][i].set_xlim(probe_wavelength.value-1, probe_wavelength.value+1)\n",
+    "    ax[1][i].set_ylim(0, 1.1*np.max(Skw.value))\n",
+    "    ax[1][i].set_xlabel('$\\lambda$ (nm)')\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plots of the spectral density function (Skw) which determines the amount of light scattered into different wavelengths.\n",
+    "\n",
+    "a. In the non-collective regime only the electron feature is visible.\n",
+    "\n",
+    "b. In the weakly collective regime (alpha approaches 1) an ion feature starts to appear and the electron feature is distorted\n",
+    "\n",
+    "c. In the collective regime both features split into two peaks, corresponding to scattering off of forward and backwards propagating plasma oscillations. \n",
+    "\n",
+    "d. The introduction of drift velocities introduces several Doppler shifts in the calculations, resulting in a shifted spectrum.\n",
+    "\n",
+    "e. Changing the ion species present changes the resulting ion feature.\n"
+   ]
+  }
+ ],
+ "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/diagnostics/__init__.py b/plasmapy/diagnostics/__init__.py
index 9d8b717c..79305478 100644
--- a/plasmapy/diagnostics/__init__.py
+++ b/plasmapy/diagnostics/__init__.py
@@ -1,3 +1,7 @@
-"""The diagnostics subpackage contains tools for experimental research.
-Currently, we have functionality for analyzing data from Langmuir probes.
 """
+The diagnostics subpackage contains functionality for defining diagnostic parameters
+and processing data collected by diagnostics, synthetic or experimental.
+"""
+__all__ = ["langmuir", "thomson"]
+
+from plasmapy.diagnostics import langmuir, thomson
diff --git a/plasmapy/diagnostics/thomson.py b/plasmapy/diagnostics/thomson.py
new file mode 100644
index 00000000..41d29e6c
--- /dev/null
+++ b/plasmapy/diagnostics/thomson.py
@@ -0,0 +1,277 @@
+"""
+Defines the Thomson scattering analysis module as
+part of the diagnostics package.
+"""
+
+__all__ = [
+    "spectral_density",
+]
+
+import astropy.constants as const
+import astropy.units as u
+import numpy as np
+
+from plasmapy.formulary.parameters import plasma_frequency, thermal_speed
+from plasmapy.formulary.dielectric import permittivity_1D_Maxwellian
+from plasmapy.particles import Particle
+from plasmapy.utils.decorators import validate_quantities
+from typing import List, Tuple, Union
+
+
+# TODO: interface for inputting a multi-species configuration could be
+# simplified using the plasmapy.classes.plasma_base class if that class
+# included ion and electron drift velocities and information about the ion
+# atomic species.
+
+
+@validate_quantities(
+    wavelengths={"can_be_negative": False},
+    probe_wavelength={"can_be_negative": False},
+    ne={"can_be_negative": False},
+    Te={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+    Ti={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+)
+def spectral_density(
+    wavelengths: u.nm,
+    probe_wavelength: u.nm,
+    ne: u.m ** -3,
+    Te: u.K,
+    Ti: u.K,
+    fract: np.ndarray = np.ones(1),
+    ion_species: Union[str, List[str], Particle, List[Particle]] = "H+",
+    fluid_vel: u.m / u.s = np.zeros(3) * u.m / u.s,
+    ion_vel: u.m / u.s = None,
+    probe_vec=np.array([1, 0, 0]),
+    scatter_vec=np.array([0, 1, 0]),
+) -> Tuple[Union[np.floating, np.ndarray], np.ndarray]:
+    r"""
+    Calculate the spectral density function for Thomson scattering of a
+    probe laser beam by a multi-species Maxwellian plasma.
+
+    This function calculates the spectral density function for Thomson
+    scattering of a probe laser beam by a plasma consisting of one or more ion
+    species and a neutralizing electron fluid:
+
+    .. math::
+        S(k,\omega) = 1 - \frac{2\pi}{k}
+        \bigg |1 - \frac{\chi_e}{\epsilon} \bigg |^2
+        f_{e0} \bigg (\frac{\omega}{k} \bigg ) +
+        \sum_i \frac{2\pi Z_i}{k}
+        \bigg |\frac{\chi_e}{\epsilon} \bigg |^2 f_{i0,i}
+        \bigg ( \frac{\omega}{k} \bigg )
+
+    where :math:`\chi_e` is the electron component susceptibility of the
+    plasma and :math:`\epsilon = 1 + \chi_e + \sum_i \chi_i` is the total
+    plasma dielectric  function (with :math:`\chi_i` being the ion component
+    of the susceptibility), :math:`Z_i` is the charge of each ion, :math:`k`
+    is the scattering wavenumber, :math:`\omega` is the scattering frequency,
+    and :math:`f_{e0}` and :math:`f_{i0,i}` are the electron and ion velocity
+    distribution functions respectively. In this function the electron and ion
+    velocity distribution functions are assumed to be Maxwellian, making this
+    function equivalent to Eq. 3.4.6 in `Sheffield`_.
+
+    Parameters
+    ----------
+
+    wavelengths : `~astropy.units.Quantity`
+        Array of wavelengths over which the spectral density function
+        will be calculated. (convertible to nm)
+
+    probe_wavelength : `~astropy.units.Quantity`
+        Wavelength of the probe laser. (convertible to nm)
+
+    ne : `~astropy.units.Quantity`
+        Mean (0th order) electron density of all plasma components combined.
+        (convertible to cm^-3.)
+
+    Te : `~astropy.units.Quantity`
+        Temperature of the electron component. (in K or convertible to eV)
+
+    Ti : `~astropy.units.Quantity`, shape (N, )
+        Temperature of each ion component. Shape (N, ) must be equal to the
+        number of ion components N. (in K or convertible to eV)
+
+    fract : array_like, shape (N, ), optional
+        An array-like object where each element represents the fraction (or ratio)
+        of the ion component number density to the total ion number density.
+        Must sum to 1.0. Default is a single ion species.
+
+    ion_species : str or `~plasmapy.particles.Particle`, shape (N, ), optional
+        A list or single instance of `~plasmapy.particles.Particle`, or strings
+        convertible to `~plasmapy.particles.Particle`. Default is `'H+'`
+        corresponding to a single species of hydrogen ions.
+
+    fluid_vel : `~astropy.units.Quantity`, shape (3, ), optional
+        Electron fluid velocity in the rest frame. (convertible to m/s)
+        Defaults to a stationary plasma [0, 0, 0] m/s.
+
+    ion_vel : `~astropy.units.Quantity`, shape (N, 3), optional
+        Velocity vectors for each ion population relative to the
+        fluid velocity. (convertible to m/s) Defaults zero drift
+        for all specified ion species.
+
+    probe_vec : float `~numpy.ndarray`, shape (3, )
+        Unit vector in the direction of the probe laser. Defaults to
+        [1, 0, 0].
+
+    scatter_vec : float `~numpy.ndarray`, shape (3, )
+        Unit vector pointing from the scattering volume to the detector.
+        Defaults to [0, 1, 0] which, along with the default `probe_vec`,
+        corresponds to a 90 degree scattering angle geometry.
+
+    Returns
+    -------
+    alpha : float
+        Mean scattering parameter, where `alpha` > 1 corresponds to collective
+        scattering and `alpha` < 1 indicates non-collective scattering.
+
+    Skw : `~astropy.units.Quantity`
+        Computed spectral density function over the input `wavelengths` array
+        with units of s/rad.
+
+    Notes
+    -----
+
+    For details, see "Plasma Scattering of Electromagnetic Radiation" by
+    Sheffield et al. `ISBN 978\\-0123748775`_. This code is a modified version
+    of the program described therein.
+
+    For a concise summary of the relevant physics, see Chapter 5 of Derek
+    Schaeffer's thesis, DOI: `10.5281/zenodo.3766933`_.
+
+    .. _`ISBN 978\\-0123748775`: https://www.sciencedirect.com/book/9780123748775/plasma-scattering-of-electromagnetic-radiation
+    .. _`10.5281/zenodo.3766933`: https://doi.org/10.5281/zenodo.3766933
+    .. _`Sheffield`: https://doi.org/10.1016/B978-0-12-374877-5.00003-8
+    """
+    fract = np.asarray(fract, dtype=np.float64)
+
+    # If ion drift velocity is not specified, create an array corresponding
+    # to zero drift
+    if ion_vel is None:
+        ion_vel = np.zeros([fract.size, 3]) * u.m / u.s
+
+    # Condition ion_species
+    if isinstance(ion_species, (str, Particle)):
+        ion_species = [ion_species]
+    if len(ion_species) == 0:
+        raise ValueError("At least one ion species needs to be defined.")
+    for ii, ion in enumerate(ion_species):
+        if isinstance(ion, Particle):
+            continue
+        ion_species[ii] = Particle(ion)
+
+    # Condition Ti
+    if Ti.size == 1:
+        # If a single quantity is given, put it in an array so it's iterable
+        # If Ti.size != len(ion_species), assume same temp. for all species
+        Ti = [Ti.value] * len(ion_species) * Ti.unit
+    elif Ti.size != len(ion_species):
+        raise ValueError(
+            f"Got {Ti.size} ion temperatures and expected " f"{len(ion_species)}."
+        )
+
+    if (
+        (len(ion_species) != fract.size)
+        or (ion_vel.shape[0] != fract.size)
+        or (Ti.size != fract.size)
+    ):
+        raise ValueError(
+            "Inconsistent number of species in fract, "
+            "ion_species, Ti, and/or ion_vel."
+        )
+
+    # Ensure unit vectors are normalized
+    probe_vec = probe_vec / np.linalg.norm(probe_vec)
+    scatter_vec = scatter_vec / np.linalg.norm(scatter_vec)
+
+    # Define some constants
+    C = const.c.si  # speed of light
+
+    # Calculate plasma parameters
+    vTe = thermal_speed(Te, particle="e-")
+    vTi, ion_z = [], []
+    for T, ion in zip(Ti, ion_species):
+        vTi.append(thermal_speed(T, particle=ion).value)
+        ion_z.append(ion.integer_charge * u.dimensionless_unscaled)
+    vTi = vTi * vTe.unit
+    zbar = np.sum(fract * ion_z)
+    ni = fract * ne / zbar  # ne/zbar = sum(ni)
+    wpe = plasma_frequency(n=ne, particle="e-")
+
+    # Compute the ion velocity in the rest frame
+    ion_vel = fluid_vel + ion_vel
+
+    # Convert wavelengths to angular frequencies (electromagnetic waves, so
+    # phase speed is c)
+    ws = (2 * np.pi * u.rad * C / wavelengths).to(u.rad / u.s)
+    wl = (2 * np.pi * u.rad * C / probe_wavelength).to(u.rad / u.s)
+
+    # Compute the frequency shift (required by energy conservation)
+    w = ws - wl
+
+    # Compute the wavenumbers in the plasma
+    # See Sheffield Sec. 1.8.1 and Eqs. 5.4.1 and 5.4.2
+    ks = np.sqrt(ws ** 2 - wpe ** 2) / C
+    kl = np.sqrt(wl ** 2 - wpe ** 2) / C
+
+    # Compute the wavenumber shift (required by momentum conservation)
+    scattering_angle = np.arccos(np.dot(probe_vec, scatter_vec))
+    # Eq. 1.7.10 in Sheffield
+    k = np.sqrt(ks ** 2 + kl ** 2 - 2 * ks * kl * np.cos(scattering_angle))
+    # Normal vector along k
+    k_vec = (scatter_vec - probe_vec) * u.dimensionless_unscaled
+
+    # Compute Doppler-shifted frequencies for both the ions and electrons
+    # Matmul is simultaneously conducting dot product over all wavelengths
+    # and ion components
+    w_e = w - k * np.dot(fluid_vel, k_vec)
+    w_i = w - np.matmul(ion_vel, np.outer(k, k_vec).T)
+
+    # Compute the scattering parameter alpha
+    # expressed here using the fact that v_th/w_p = root(2) * Debye length
+    alpha = np.sqrt(2) * wpe / (k * vTe)
+
+    # Calculate the normalized phase velocities (Sec. 3.4.2 in Sheffield)
+    xe = (w_e / (k * vTe)).to(u.dimensionless_unscaled)
+    xi = (np.outer(1 / vTi, 1 / k) * w_i).to(u.dimensionless_unscaled)
+
+    # Calculate the susceptibilities
+    chiE = permittivity_1D_Maxwellian(w_e, k, Te, ne, "e-")
+
+    # Treatment of multiple species is an extension of the discussion in
+    # Sheffield Sec. 5.1
+    chiI = np.zeros([fract.size, w.size], dtype=np.complex128)
+    for i, ion in enumerate(ion_species):
+        chiI[i, :] = permittivity_1D_Maxwellian(
+            w_i[i, :], k, Ti[i], ni[i], ion, z_mean=ion_z[i]
+        )
+
+    # Calculate the longitudinal dielectric function
+    epsilon = 1 + chiE + np.sum(chiI, axis=0)
+
+    # Calculate the contributions to the spectral density function
+    econtr = (
+        2
+        * np.sqrt(np.pi)
+        / k
+        / vTe
+        * np.power(np.abs(1 - chiE / epsilon), 2)
+        * np.exp(-(xe ** 2))
+    )
+
+    icontr = np.zeros([fract.size, w.size], dtype=np.complex128) * u.s / u.rad
+    for m in range(fract.size):
+        icontr[m, :] = (
+            2
+            * np.sqrt(np.pi)
+            * ion_z[m]
+            / k
+            / vTi[m]
+            * np.power(np.abs(chiE / epsilon), 2)
+            * np.exp(-xi[m, :] ** 2)
+        )
+
+    # Recast as real: imaginary part is already zero
+    Skw = np.real(econtr + np.sum(icontr, axis=0))
+
+    return np.mean(alpha), Skw
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(
diff --git a/requirements/requirements.txt b/requirements/requirements.txt
index 0a968d30..62af843d 100644
--- a/requirements/requirements.txt
+++ b/requirements/requirements.txt
@@ -5,6 +5,6 @@ h5py >= 2.8
 lmfit >= 1.0.1
 matplotlib >= 2.0
 mpmath >= 1.0.0
-numpy >= 1.16.6
+numpy >= 1.18.1
 scipy >= 1.2
 setuptools >= 41.2
diff --git a/setup.cfg b/setup.cfg
index 7f25475c..d8b3f845 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -32,7 +32,7 @@ setup_requires = setuptools_scm
 install_requires =
   astropy >= 4.0
   colorama >= 0.3
-  numpy >= 1.16.6
+  numpy >= 1.18.1
   scipy >= 1.2
   setuptools >= 41.2
 
diff --git a/tox.ini b/tox.ini
index e090249c..c46a57a4 100644
--- a/tox.ini
+++ b/tox.ini
@@ -45,7 +45,7 @@ deps =
   pytest-cov
   pytest-xdist
 conda_deps =
-  numpy>=1.16.6
+  numpy>=1.18.1
   scipy>=1.2
   astropy>=4.0
   colorama>=0.3
@@ -68,7 +68,7 @@ deps =
   pytest-cov
   pytest-xdist
 conda_deps =
-  numpy=1.16.6
+  numpy=1.18.1
   scipy=1.2
   astropy=4.0
   colorama=0.3


Test patch:
diff --git a/plasmapy/diagnostics/tests/test_thomson.py b/plasmapy/diagnostics/tests/test_thomson.py
new file mode 100644
index 00000000..e16b095e
--- /dev/null
+++ b/plasmapy/diagnostics/tests/test_thomson.py
@@ -0,0 +1,251 @@
+"""
+Tests for Thomson scattering analysis functions
+"""
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+from plasmapy.diagnostics import thomson
+from plasmapy.particles import Particle
+
+
+def width_at_value(x, y, val):
+    """
+    Calculates the width of a curve at a given value.
+    """
+    above = np.where(y > val, x, np.nan)
+    return np.abs(np.nanmax(above) - np.nanmin(above))
+
+
+def gen_collective_spectrum():
+    """
+    Generates an example Thomson scattering spectrum in the collective regime
+    """
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 10 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = ["C-12 5+"]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def gen_multiple_ion_species_spectrum():
+    """
+    Generates an example Thomson scattering spectrum for multiple ion species
+    that also have drift velocities. Parameters are set to be in the
+    collective regime where ion species are important.
+    """
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([0.7, 0.3])
+    Te = 10 * u.eV
+    Ti = np.array([5, 5]) * u.eV
+    fluid_vel = np.array([300, 0, 0]) * u.km / u.s
+    ion_vel = np.array([[-500, 0, 0], [0, 500, 0]]) * u.km / u.s
+
+    # Use this to also test passing in ion species as Particle objects
+    ion_species = [Particle("p+"), Particle("C-12 5+")]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+        fluid_vel=fluid_vel,
+        ion_vel=ion_vel,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def gen_non_collective_spectrum():
+    """
+    Generates an example Thomson scattering spectrum in the non-collective
+    regime
+    """
+    wavelengths = np.arange(500, 570, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e15 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 100 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = ["H+"]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def test_different_input_types():
+
+    # Define some constants
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 10 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = "C-12 5+"
+
+    # Raise a ValueError with inconsistent ion array lengths
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            Ti,
+            fract=np.array([0.5, 0.5]),
+            ion_species=ion_species,
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+    # Raise a ValueError with inconsistent ion temperature array
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            np.array([5, 5]) * u.eV,
+            fract=fract,
+            ion_species=ion_species,
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+    # Raise a ValueError with empty ion_species
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            Ti,
+            fract=fract,
+            ion_species=[],
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+
+def test_collective_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+    alpha, wavelength, Skw = gen_collective_spectrum()
+
+    # Check that alpha is correct
+    assert np.isclose(alpha.value, 1.801, atol=0.01), (
+        "Collective case alpha " f"returns {alpha} instead of " "expected 1.801"
+    )
+
+    i_width = width_at_value(wavelength.value, Skw.value, 2e-13)
+    e_width = width_at_value(wavelength.value, Skw.value, 0.2e-13)
+
+    # Check that the widths of the ion and electron features match expectations
+    assert np.isclose(i_width, 0.1599, 1e-3), (
+        "Collective case ion feature "
+        f"width is {i_width}"
+        "instead of expected 0.1599"
+    )
+
+    assert np.isclose(e_width, 17.7899, 1e-3), (
+        "Collective case electron "
+        f"feature width is {e_width} "
+        "instead of expected 17.7899"
+    )
+
+
+def test_non_collective_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+    alpha, wavelength, Skw = gen_non_collective_spectrum()
+
+    # Check that alpha is correct
+    assert np.isclose(alpha.value, 0.05707, atol=0.01), (
+        "Non-collective case alpha " f"returns {alpha} instead of " "expected 0.05707"
+    )
+
+    e_width = width_at_value(wavelength.value, Skw.value, 0.2e-13)
+
+    # Check that the widts of the electron feature matchs expectations
+    assert np.isclose(e_width, 22.6699, 1e-3), (
+        "Non-collective case electron "
+        f"feature width is {e_width} "
+        "instead of expected 22.6699"
+    )
+
+
+def test_multiple_ion_species_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+
+    alpha, wavelength, Skw = gen_multiple_ion_species_spectrum()
+
+    # Compute the width and max of the spectrum, and the wavelength
+    # of the max (sensitive to ion vel)
+    width = width_at_value(wavelength.value, Skw.value, 0.2e-11)
+    max_skw = np.max(Skw.value)
+    max_wavelength = wavelength.value[np.argmax(Skw.value)]
+
+    # Check width
+    assert np.isclose(width, 0.1599, 1e-3), (
+        f"Multiple ion species case spectrum width is {width} instead of "
+        "expected 0.1599"
+    )
+
+    # Check max value
+    assert np.isclose(max_skw, 2.4e-11, 1e-11), (
+        f"Multiple ion species case spectrum max is {max_skw} instead of "
+        "expected 2.4e-11"
+    )
+
+    # Check max peak location
+    assert np.isclose(max_wavelength, 531.549, 1e-3), (
+        "Multiple ion species case spectrum peak wavelength is "
+        f"{max_wavelength} instead of expected 531.549"
+    )
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)