Repository: PlasmaPy/PlasmaPy
Pull request #860 (2020-06-23 15:59:40)

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


codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/860?src=pr&el=h1) Report
> Merging [#860](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/860?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/983c2f841d7690643e9568b71e649bb4859ff70f&el=desc) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

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

```diff
@@           Coverage Diff           @@
##           master     #860   +/-   ##
=======================================
  Coverage   96.00%   96.01%           
=======================================
  Files          59       59           
  Lines        5209     5216    +7     
=======================================
+ Hits         5001     5008    +7     
  Misses        208      208           
```


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

------

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

StanczakDominik: I kind of sort of thought the two - ionization balance and the Saha equation were the same thing? :open_mouth: 
StanczakDominik: I went ahead and added `ionization` to the `formulary` doc index, so that that module is actually visible in there.

Patch:
diff --git a/changelog/860.feature.rst b/changelog/860.feature.rst
new file mode 100644
index 00000000..e26088d2
--- /dev/null
+++ b/changelog/860.feature.rst
@@ -0,0 +1,1 @@
+ Create :func:`~plasmapy.formulary.ionization.Saha` function.
diff --git a/changelog/862.trivial.rst b/changelog/862.trivial.rst
new file mode 100644
index 00000000..95d8678e
--- /dev/null
+++ b/changelog/862.trivial.rst
@@ -0,0 +1,2 @@
+Update :class:`~plasmapy.particles.particle_class.CustomParticle` so input parameters
+`mass` and `charge` can accept string representations of astropy `Quantities`.
diff --git a/docs/formulary/index.rst b/docs/formulary/index.rst
index 86ac44dd..bfdb7580 100644
--- a/docs/formulary/index.rst
+++ b/docs/formulary/index.rst
@@ -24,6 +24,7 @@ physical quantities helpful for plasma physics.
    relativity
    braginskii
    collisions
+   ionization
 
 The subpackage makes heavy use of `astropy.units.Quantity` for handling
 conversions between different unit systems. This is especially important
diff --git a/docs/formulary/ionization.rst b/docs/formulary/ionization.rst
new file mode 100644
index 00000000..4de60666
--- /dev/null
+++ b/docs/formulary/ionization.rst
@@ -0,0 +1,10 @@
+.. _ionization:
+
+******************************************************************
+Ionization related functionality (`plasmapy.formulary.ionization`)
+******************************************************************
+
+.. currentmodule:: plasmapy.formulary.ionization
+
+.. automodapi:: plasmapy.formulary.ionization
+   :no-heading:
diff --git a/plasmapy/formulary/__init__.py b/plasmapy/formulary/__init__.py
index 82a77298..ef4f4efc 100644
--- a/plasmapy/formulary/__init__.py
+++ b/plasmapy/formulary/__init__.py
@@ -15,3 +15,4 @@
 from .parameters import *
 from .quantum import *
 from .relativity import *
+from .ionization import *
diff --git a/plasmapy/formulary/ionization.py b/plasmapy/formulary/ionization.py
index 5909b861..d999fc29 100644
--- a/plasmapy/formulary/ionization.py
+++ b/plasmapy/formulary/ionization.py
@@ -1,15 +1,16 @@
 """ This module gathers functions relating to ionization states and the properties thereof.
 """
-__all__ = ["ionization_balance"]
+__all__ = ["ionization_balance", "Saha", "Z_bal_"]
+
+import astropy.units as u
+
+from astropy.constants import a0, k_B
+from numpy import pi, exp, sqrt, log
 
-from astropy import units as u
 from plasmapy.utils.decorators import (
-    angular_freq_to_hz,
     check_relativistic,
     validate_quantities,
 )
-from astropy.constants import c, a0, k_B
-from numpy import pi, exp, sqrt, log
 
 
 @validate_quantities(
@@ -36,6 +37,8 @@ def ionization_balance(n: u.m ** -3, T_e: u.K) -> u.dimensionless_unscaled:
     :math:`a_0` is the Bohr radius, and
     :math:`E_H` is the ionization energy of Hydrogen
 
+    **Aliases:** `Z_bal_`
+
     Parameters
     ----------
     T_e : `~astropy.units.Quantity`
@@ -46,7 +49,7 @@ def ionization_balance(n: u.m ** -3, T_e: u.K) -> u.dimensionless_unscaled:
 
     Warns
     -----
-    ~astropy.units.UnitsWarning
+    : `~astropy.units.UnitsWarning`
         If units are not provided, SI units are assumed.
 
     Raises
@@ -64,11 +67,11 @@ def ionization_balance(n: u.m ** -3, T_e: u.K) -> u.dimensionless_unscaled:
     >>> T_e = 5000 * u.K
     >>> n = 1e19 * u.m ** -3
     >>> ionization_balance(n, T_e)
-    <Quantity 0.27478417>
+    <Quantity 0.274...>
     >>> T_e = 50 * u.eV
     >>> n = 1e10 * u.m ** -3
     >>> ionization_balance(n, T_e)
-    <Quantity 12.61576334>
+    <Quantity 12.615...>
 
     Returns
     -------
@@ -86,5 +89,96 @@ def ionization_balance(n: u.m ** -3, T_e: u.K) -> u.dimensionless_unscaled:
     return Z
 
 
-Z_bal = ionization_balance
-"""alias for :func:`ionization_balance`"""
+Z_bal_ = ionization_balance
+"""Alias for :func:`ionization_balance`."""
+
+
+@validate_quantities(
+    T_e={"can_be_negative": False, "equivalencies": u.temperature_energy()}
+)
+def Saha(g_j, g_k, n_e: u.m ** -3, E_jk: u.J, T_e: u.K) -> u.dimensionless_unscaled:
+    r"""
+    The Saha equation, derived in statistical mechanics, gives an approximation of
+    the ratio of population of ions in two different ionization states in a plasma.
+    This approximation applies to plasmas in thermodynamic equilibrium where ionization
+    and recombination of ions with electrons are balanced.
+
+    .. math::
+        \frac{N_j}{N_k} = \frac{1}{n_e} \frac{g_j}{4 g_k a_0^{3}} \left( \frac{k_B T_e}{\pi E_H} \right)^{\frac{3}{2}}
+         \exp\left( \frac{-E_{jk}}{k_B T_e} \right)
+
+    Where :math:`k_B` is the Boltzmann constant,
+    :math:`a_0` is the Bohr radius,
+    :math:`E_H` is the ionization energy of Hydrogen,
+    :math:`N_j` and :math:`N_k` are the population of ions in the j and k states respectively.
+    This function is equivalent to Eq. 3.47 in `Drake`_.
+
+    Parameters
+    ----------
+    T_e : `~astropy.units.Quantity`
+        The electron temperature.
+
+    g_j : `int`
+        The degeneracy of the 'j'th ionization state.
+
+    g_k : `int`
+        The degeneracy of the 'k'th ionization state.
+
+    E_jk : `~astropy.units.Quantity`
+        The energy difference between ionization states j and k.
+
+    n_e : `~astropy.units.Quantity`
+        The electron number density of the plasma.
+
+    Warns
+    -----
+    : `~astropy.units.UnitsWarning`
+        If units are not provided, SI units are assumed.
+
+    Raises
+    ------
+    TypeError
+        The `T_e`, `E_jk`, or `n_e` is not a `~astropy.units.Quantity` and cannot be
+        converted into a ~astropy.units.Quantity.
+
+    ~astropy.units.UnitConversionError
+        If the `T_e`, `E_jk`, or `n` not in appropriate units.
+
+    Returns
+    -------
+    ratio: `~astropy.units.Quantity`
+        The ratio of population of ions in ionization state j to state k.
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> T_e = 5000 * u.K
+    >>> n = 1e19 * u.m ** -3
+    >>> g_j = 2
+    >>> g_k = 2
+    >>> E_jk = 1 * u.Ry
+    >>> Saha(g_j, g_k, n, E_jk, T_e)
+    <Quantity 3.299...e-06>
+    >>> T_e = 1 * u.Ry
+    >>> n = 1e23 * u.m ** -3
+    >>> Saha(g_j, g_k, n, E_jk, T_e)
+    <Quantity 1114595.586...>
+
+    Notes
+    -----
+    For reference to this function and for more information regarding the Saha equation,
+    see chapter 3 of R Paul Drake's book, "High-Energy-Density Physics: Foundation of
+    Inertial Fusion and Experimental Astrophysics" (`DOI: 10.1007/978-3-319-67711-8_3`_).
+
+    .. _`Drake`: https://doi.org/10.1007/978-3-319-67711-8
+    .. _`DOI: 10.1007/978-3-319-67711-8_3`: https://doi.org/10.1007/978-3-319-67711-8_3
+    """
+    E_h = 1 * u.Ry
+
+    degeneracy_factor = (1 / n_e) * g_j / (4 * g_k * a0 ** 3)
+    physical_constants = (k_B * T_e / (pi * E_h)) ** (3 / 2)
+    boltzmann_factor = exp(-E_jk / (k_B * T_e))
+
+    ratio = degeneracy_factor * physical_constants * boltzmann_factor
+
+    return ratio
diff --git a/plasmapy/particles/particle_class.py b/plasmapy/particles/particle_class.py
index 179dfe37..60fc1452 100644
--- a/plasmapy/particles/particle_class.py
+++ b/plasmapy/particles/particle_class.py
@@ -1832,15 +1832,18 @@ def charge(self) -> u.C:
     @mass.setter
     def mass(self, m: u.kg):
         if m is None:
-            self._mass = np.nan * u.kg
+            m = np.nan * u.kg
             warnings.warn("CustomParticle mass set to NaN kg", MissingAtomicDataWarning)
+        elif isinstance(m, str):
+            m = u.Quantity(m)
         elif not isinstance(m, u.Quantity):
             raise TypeError(
                 "The mass of a custom particle must be a nonnegative Quantity "
                 "with units of mass."
             )
+        if np.isnan(m):
+            self._mass = m
         else:
-
             if not isinstance(m.value, Real):
                 raise TypeError(
                     "The mass of a custom particle must be a real number "
@@ -1858,10 +1861,15 @@ def mass(self, m: u.kg):
     @charge.setter
     def charge(self, q: Optional[Union[u.Quantity, Real]]):
         if q is None:
-            self._charge = np.nan * u.C
+            q = np.nan * u.C
             warnings.warn(
                 "CustomParticle charge set to NaN C", MissingAtomicDataWarning
             )
+        elif isinstance(q, str):
+            q = u.Quantity(q)
+
+        if np.isnan(q):
+            self._charge = q
         elif isinstance(q, Real):
             self._charge = q * const.e.si
             warnings.warn(


Test patch:
diff --git a/plasmapy/formulary/tests/test_ionization.py b/plasmapy/formulary/tests/test_ionization.py
index 3df0aa45..1d87c4e3 100644
--- a/plasmapy/formulary/tests/test_ionization.py
+++ b/plasmapy/formulary/tests/test_ionization.py
@@ -1,19 +1,41 @@
 import astropy.units as u
 import numpy as np
 import pytest
-from astropy.constants import c
 
-from plasmapy.formulary.ionization import Z_bal
+from plasmapy.formulary.ionization import Saha, ionization_balance, Z_bal_
 
 
-def test_Z_bal():
+def test_ionization_balance():
     n = 1e19 * u.m ** -3
     T_e = 5000 * u.K
 
-    assert (Z_bal(n, T_e) * u.dimensionless_unscaled).unit == u.dimensionless_unscaled
+    val = ionization_balance(n, T_e)
+    assert val.unit == u.dimensionless_unscaled
+    assert np.isclose(val.value, 0.27478417234035346)
+
+    # aliases must be the same function
+    assert Z_bal_ is ionization_balance
+
+    with pytest.warns(u.UnitsWarning):
+        ionization_balance(1e19, T_e)
+
+    with pytest.raises(u.UnitTypeError):
+        ionization_balance(2e19 * u.kg, T_e)
+
+
+def test_Saha():
+    g_j = 2
+    g_k = 1
+    n_e = 1e19 * u.m ** -3
+    E_jk = 10 * u.Ry
+    T_e = 5000 * u.K
+
+    val = Saha(g_j, g_k, n_e, E_jk, T_e)
+    assert val.unit == u.dimensionless_unscaled
+    assert np.isclose(val.value, 2.4775608756800435e-129)
 
     with pytest.warns(u.UnitsWarning):
-        Z_bal(1e19, T_e)
+        Saha(g_j, g_k, 1e19, E_jk, T_e)
 
     with pytest.raises(u.UnitTypeError):
-        Z_bal(2e19 * u.kg, T_e)
+        Saha(g_j, g_k, 1e19 * u.kg, E_jk, T_e)
diff --git a/plasmapy/particles/tests/test_particle_class.py b/plasmapy/particles/tests/test_particle_class.py
index 25ae4b92..89dd7c16 100644
--- a/plasmapy/particles/tests/test_particle_class.py
+++ b/plasmapy/particles/tests/test_particle_class.py
@@ -905,7 +905,9 @@ def test_that_object_can_be_dict_key(key):
     (CustomParticle, {"charge": -0.1 * u.C}, "charge", -0.1 * u.C),
     (CustomParticle, {"charge": -2}, "charge", -2 * const.e.si),
     (CustomParticle, {"mass": np.inf * u.g}, "mass", np.inf * u.kg),
+    (CustomParticle, {"mass": "100.0 g"}, "mass", 100.0 * u.g),
     (CustomParticle, {"charge": -np.inf * u.kC}, "charge", -np.inf * u.C),
+    (CustomParticle, {"charge": "5.0 C"}, "charge", 5.0 * u.C),
 ]
 
 
@@ -947,12 +949,16 @@ def test_customized_particles(cls, kwargs, attr, expected):
     (CustomParticle, {"charge": np.complex128(5 + 2j)}, InvalidParticleError),
     (CustomParticle, {"mass": np.complex128(5 + 2j)}, InvalidParticleError),
     (CustomParticle, {"mass": -1e-36 * u.kg}, InvalidParticleError),
+    (CustomParticle, {"mass": "not a mass"}, InvalidParticleError),
+    (CustomParticle, {"mass": "5.0 km"}, InvalidParticleError),
     (CustomParticle, {"mass": np.array([1, 1]) * u.kg}, InvalidParticleError),
     (CustomParticle, {"charge": np.array([1, 1]) * u.C}, InvalidParticleError),
     (CustomParticle, {"charge": (5 + 2j) * u.C}, InvalidParticleError),
     (CustomParticle, {"mass": (5 + 2j) * u.kg}, InvalidParticleError),
     (CustomParticle, {"charge": np.complex128(5 + 2j) * u.C}, InvalidParticleError),
     (CustomParticle, {"mass": np.complex128(5 + 2j) * u.kg}, InvalidParticleError),
+    (CustomParticle, {"charge": "not a charge"}, InvalidParticleError),
+    (CustomParticle, {"charge": "5.0 km"}, InvalidParticleError),
 ]