Repository: PlasmaPy/PlasmaPy
Pull request #388 (2018-04-21 05:03:13)

[WIP] Fix transport issues
Fixes a bunch of them. I'll update this post.


Patch:
diff --git a/plasmapy/classes/plasma.py b/plasmapy/classes/plasma.py
index 431bbbac..614e9d9d 100644
--- a/plasmapy/classes/plasma.py
+++ b/plasmapy/classes/plasma.py
@@ -1,6 +1,7 @@
 """
 Defines the core Plasma class used by PlasmaPy to represent plasma properties.
 """
+import warnings
 
 import numpy as np
 import astropy.units as u
@@ -23,7 +24,7 @@
                              integer_charge,
                              )
 
-from plasmapy.utils import call_string
+from plasmapy.utils import call_string, CouplingWarning
 
 
 class Plasma3D:
@@ -225,11 +226,15 @@ def coupling(self):
         are important. This compares Coulomb potential energy to thermal
         kinetic energy.
         """
-        return coupling_parameter(self.T_e,
-                                  self.n_e,
-                                  (self.particle, self.particle),
-                                  self.Z)
-
+        couple = coupling_parameter(self.T_e,
+                                    self.n_e,
+                                    (self.particle, self.particle),
+                                    self.Z)
+        if couple < 0.01:
+            warnings.warn(f"Coupling parameter is {couple}, you might have strong coupling effects",
+                          CouplingWarning)
+
+        return couple
     def quantum_theta(self):
         """
         Quantum theta parameter, which compares Fermi kinetic energy to
diff --git a/plasmapy/physics/transport/__init__.py b/plasmapy/physics/transport/__init__.py
index de8b6c5e..5cc4550b 100644
--- a/plasmapy/physics/transport/__init__.py
+++ b/plasmapy/physics/transport/__init__.py
@@ -2,7 +2,7 @@
                          collision_frequency,
                          collision_rate_electron_ion,
                          collision_rate_ion_ion,
-                         b_perp,
+                         impact_parameter_perp,
                          impact_parameter,
                          mean_free_path,
                          mobility,
diff --git a/plasmapy/physics/transport/collisions.py b/plasmapy/physics/transport/collisions.py
index 82d82b60..b986e6ec 100644
--- a/plasmapy/physics/transport/collisions.py
+++ b/plasmapy/physics/transport/collisions.py
@@ -138,6 +138,45 @@ def Coulomb_logarithm(T,
     This means the Coulomb logarithm will not break down for Lambda < 0,
     which occurs for dense, cold plasmas.
 
+    Methods
+    ---
+    Classical
+        classical Landau-Spitzer approach. Fails for large coupling
+        parameter where Lambda can become less than zero.
+    GMS-1
+        1st method listed in Table 1 of reference [3]
+        Landau-Spitzer, but with interpolated bmin instead of bmin
+        selected between deBroglie wavelength and distance of closest
+        approach. Fails for large coupling
+        parameter where Lambda can become less than zero.
+    GMS-2
+        2nd method listed in Table 1 of reference [3]
+        Another Landau-Spitzer like approach, but now bmax is also
+        being interpolated. The interpolation is between the Debye
+        length and the ion sphere radius, allowing for descriptions
+        of dilute plasmas. Fails for large coupling
+        parameter where Lambda can become less than zero.
+        3rd method listed in Table 1 of reference [3]
+        classical Landau-Spitzer fails for argument of Coulomb logarithm
+        Lambda < 0, therefore a clamp is placed at Lambda_min = 2
+    GMS-4
+        4th method listed in Table 1 of reference [3]
+        Spitzer-like extension to Coulomb logarithm by noting that
+        Coulomb collisions take hyperbolic trajectories. Removes
+        divergence for small bmin issue in classical Landau-Spitzer
+        approach, so bmin can be zero. Also doesn't break down as
+        Lambda < 0 is now impossible, even when coupling parameter is large.
+    GMS-5
+        5th method listed in Table 1 of reference [3]
+        Similar to GMS-4, but setting bmin as distance of closest approach
+        and bmax interpolated between Debye length and ion sphere radius.
+        Lambda < 0 impossible.
+    GMS-6
+        6th method listed in Table 1 of reference [3]
+        Similar to GMS-4 and GMS-5, but using interpolation methods
+        for both bmin and bmax.
+    Lambda < 0 impossible.
+
     Examples
     --------
     >>> from astropy import units as u
@@ -172,61 +211,26 @@ def Coulomb_logarithm(T,
                                   z_mean=z_mean,
                                   V=V,
                                   method=method)
-    if method == "classical":
-        # classical Landau-Spitzer approach. Fails for large coupling
-        # parameter where Lambda can become less than zero.
-        ln_Lambda = np.log(bmax / bmin)
-    elif method == "GMS-1":
-        # 1st method listed in Table 1 of reference [3]
-        # Landau-Spitzer, but with interpolated bmin instead of bmin
-        # selected between deBroglie wavelength and distance of closest
-        # approach. Fails for large coupling
-        # parameter where Lambda can become less than zero.
-        ln_Lambda = np.log(bmax / bmin)
-    elif method == "GMS-2":
-        # 2nd method listed in Table 1 of reference [3]
-        # Another Landau-Spitzer like approach, but now bmax is also
-        # being interpolated. The interpolation is between the Debye
-        # length and the ion sphere radius, allowing for descriptions
-        # of dilute plasmas. Fails for large coupling
-        # parameter where Lambda can become less than zero.
+    if method in ["classical", "GMS-1", "GMS-2"]:
         ln_Lambda = np.log(bmax / bmin)
     elif method == "GMS-3":
-        # 3rd method listed in Table 1 of reference [3]
-        # classical Landau-Spitzer fails for argument of Coulomb logarithm
-        # Lambda < 0, therefore a clamp is placed at Lambda_min = 2
         ln_Lambda = np.log(bmax / bmin)
         if ln_Lambda < 2:
             ln_Lambda = 2 * u.dimensionless_unscaled
-    elif method == "GMS-4":
-        # 4th method listed in Table 1 of reference [3]
-        # Spitzer-like extension to Coulomb logarithm by noting that
-        # Coulomb collisions take hyperbolic trajectories. Removes
-        # divergence for small bmin issue in classical Landau-Spitzer
-        # approach, so bmin can be zero. Also doesn't break down as
-        # Lambda < 0 is now impossible, even when coupling parameter is large.
-        ln_Lambda = 0.5 * np.log(1 + bmax ** 2 / bmin ** 2)
-    elif method == "GMS-5":
-        # 5th method listed in Table 1 of reference [3]
-        # Similar to GMS-4, but setting bmin as distance of closest approach
-        # and bmax interpolated between Debye length and ion sphere radius.
-        # Lambda < 0 impossible.
-        ln_Lambda = 0.5 * np.log(1 + bmax ** 2 / bmin ** 2)
-    elif method == "GMS-6":
-        # 6th method listed in Table 1 of reference [3]
-        # Similar to GMS-4 and GMS-5, but using interpolation methods
-        # for both bmin and bmax.
-        # Lambda < 0 impossible.
+    elif method in ["GMS-4", "GMS-5", "GMS-6"]:
         ln_Lambda = 0.5 * np.log(1 + bmax ** 2 / bmin ** 2)
     else:
         raise ValueError("Unknown method! Choose from 'classical' and 'GMS-N', N from 1 to 6.")
     # applying dimensionless units
     ln_Lambda = ln_Lambda.to(u.dimensionless_unscaled).value
-    if ln_Lambda < 4:
+    if ln_Lambda < 0:
+        raise utils.PhysicsError(f"A Coulomb logarithm = {ln_Lambda} < 0 is nonphysical.")
+    elif ln_Lambda < 2 and method in ["classical", "GMS-1", "GMS-2"]:
+        warnings.warn(f"Coulomb logarithm is {ln_Lambda} and {method} relies on weak coupling.",
+                      utils.CouplingWarning)
+    elif ln_Lambda < 4:
         warnings.warn(f"Coulomb logarithm is {ln_Lambda}, you might have strong coupling effects",
-                      utils.PhysicsWarning)
-    # if ln_Lambda < 1:
-    #     raise utils.PhysicsError(f"Coulomb logarithm is {ln_Lambda}, less than 1")
+                      utils.CouplingWarning)
     return ln_Lambda
 
 
@@ -251,6 +255,8 @@ def _boilerPlate(T, particles, V):
     # obtaining reduced mass of 2 particle collision system
     reduced_mass = atomic.reduced_mass(*particles)
 
+    if V == 0:
+        raise utils.exceptions.PhysicsError("You cannot have a collision for zero velocity!")
     # getting thermal velocity of system if no velocity is given
     if np.isnan(V):
         V = parameters.thermal_speed(T, mass=reduced_mass)
@@ -260,9 +266,9 @@ def _boilerPlate(T, particles, V):
 
 @check_quantity({"T": {"units": u.K, "can_be_negative": False}
                  })
-def b_perp(T,
-           particles,
-           V=np.nan * u.m / u.s):
+def impact_parameter_perp(T,
+                          particles,
+                          V=np.nan * u.m / u.s):
     """
     Distance of closest approach for a 90 degree Coulomb collision.
 
@@ -285,7 +291,7 @@ def b_perp(T,
 
     Returns
     -------
-    b_perp : float or numpy.ndarray
+    impact_parameter_perp : float or numpy.ndarray
         The distance of closest approach for a 90 degree Coulomb collision.
 
     Raises
@@ -315,7 +321,7 @@ def b_perp(T,
 
     Notes
     -----
-    The distance of closest approach, b_perp, is given by [1]_
+    The distance of closest approach, impact_parameter_perp, is given by [1]_
 
     .. math::
         b_{\perp} = \frac{Z_1 Z_2}{4 \pi \epsilon_0 m v^2}
@@ -326,7 +332,7 @@ def b_perp(T,
     >>> from astropy import units as u
     >>> T = 1e6*u.K
     >>> particles = ('e', 'p')
-    >>> b_perp(T, particles)
+    >>> impact_parameter_perp(T, particles)
     <Quantity 8.35505011e-12 m>
 
 
@@ -474,9 +480,9 @@ def impact_parameter(T,
     # deBroglie wavelength
     lambdaBroglie = hbar / (2 * reduced_mass * V)
     # distance of closest approach in 90 degree Coulomb collision
-    bPerp = b_perp(T=T,
-                   particles=particles,
-                   V=V)
+    bPerp = impact_parameter_perp(T=T,
+                                  particles=particles,
+                                  V=V)
     # obtaining minimum and maximum impact parameters depending on which
     # method is requested
     if method == "classical":
@@ -627,14 +633,6 @@ def collision_frequency(T,
     taken as the thermal velocity), and :math:`\ln{\Lambda}` is the Coulomb
     logarithm accounting for small angle collisions.
 
-    The collisional cross-section is obtained by
-
-    .. math::
-        \sigma = \pi \rho_{\perp}^2
-
-    where :math:`\rho_{\perp}` is the distance of closest approach for
-    a 90 degree Coulomb collision.
-
     See eq (2.14) in [2]_.
 
     Examples
@@ -666,9 +664,9 @@ def collision_frequency(T,
             V = V_reduced
         # electron-electron collision
         # impact parameter for 90 degree collision
-        bPerp = b_perp(T=T,
-                       particles=particles,
-                       V=V_reduced)
+        bPerp = impact_parameter_perp(T=T,
+                                      particles=particles,
+                                      V=V_reduced)
         # Coulomb logarithm
         cou_log = Coulomb_logarithm(T,
                                     n,
@@ -686,9 +684,9 @@ def collision_frequency(T,
             V = np.sqrt(2 * k_B * T / m_e)
         # need to also correct mass in collision radius from reduced
         # mass to electron mass
-        bPerp = b_perp(T=T,
-                       particles=particles,
-                       V=V) * reduced_mass / m_e
+        bPerp = impact_parameter_perp(T=T,
+                                      particles=particles,
+                                      V=V) * reduced_mass / m_e
         # Coulomb logarithm
         # !!! may also need to correct Coulomb logarithm to be
         # electron-electron version !!!
@@ -704,9 +702,9 @@ def collision_frequency(T,
         if np.isnan(V):
             V = V_reduced
         # ion-ion collision
-        bPerp = b_perp(T=T,
-                       particles=particles,
-                       V=V)
+        bPerp = impact_parameter_perp(T=T,
+                                      particles=particles,
+                                      V=V)
         # Coulomb logarithm
         cou_log = Coulomb_logarithm(T,
                                     n,
@@ -715,7 +713,7 @@ def collision_frequency(T,
                                     V=np.nan * u.m / u.s,
                                     method=method)
     # collisional cross section
-    sigma = np.pi * (2 * bPerp) ** 2
+    sigma = Coulomb_cross_section(bPerp)
     # collision frequency where Coulomb logarithm accounts for
     # small angle collisions, which are more frequent than large
     # angle collisions.
@@ -723,6 +721,51 @@ def collision_frequency(T,
     return freq.to(u.Hz)
 
 
+@check_quantity({
+    'impact_param': {'units': u.m,
+                    'can_be_negative': False}
+    })
+def Coulomb_cross_section(impact_param: u.m):
+    """
+    Cross section for a large angle Coulomb collision
+
+    Parameters
+    ----------
+    impact_param : ~astropy.units.Quantity
+        Impact parameter for the collision.
+
+    Examples
+    --------
+    >>> Coulomb_cross_section(7e-10*u.m)
+    <Quantity 6.1575216e-18 m2>
+    >>> Coulomb_cross_section(0.5*u.m)
+    <Quantity 3.14159265 m2>
+
+    Notes
+    -----
+    The collisional cross-section (see [1]_ for a graphical demonstration)
+    for a 90 degree Coulomb collision is obtained by
+
+    .. math::
+        \sigma = \pi (2 * \rho_{\perp})^2
+
+    where :math:`\rho_{\perp}` is the distance of closest approach for
+    a 90 degree Coulomb collision. This function is a generalization of that
+    calculation. Please note that it is not guaranteed to return the correct
+    results for small angle collisions.
+
+    Returns
+    -------
+    ~astropy.units.Quantity
+        The Coulomb collision cross section area.
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Cross_section_(physics)#Collision_among_gas_particles
+    """
+    sigma = np.pi * (2 * impact_param) ** 2
+    return sigma
+
 @utils.check_quantity({
     'T_e': {'units': u.K, 'can_be_negative': False},
     'n_e': {'units': u.m ** -3, 'can_be_negative': False}
diff --git a/plasmapy/utils/__init__.py b/plasmapy/utils/__init__.py
index d4115d62..bb47e8b2 100644
--- a/plasmapy/utils/__init__.py
+++ b/plasmapy/utils/__init__.py
@@ -15,6 +15,7 @@
                          InvalidParticleError,
                          PlasmaPyWarning,
                          PhysicsWarning,
+                         CouplingWarning,
                          RelativityWarning,
                          AtomicWarning,
                          MissingAtomicDataWarning)
diff --git a/plasmapy/utils/exceptions.py b/plasmapy/utils/exceptions.py
index 4c08ac80..cefb4b26 100644
--- a/plasmapy/utils/exceptions.py
+++ b/plasmapy/utils/exceptions.py
@@ -109,6 +109,12 @@ class RelativityWarning(PhysicsWarning):
     pass
 
 
+class CouplingWarning(PhysicsWarning):
+    """
+    A warning for functions that rely on a particular coupling regime to be valid.
+    """
+
+
 class AtomicWarning(PlasmaPyWarning):
     """The base warning for the `~plasmapy.atomic` subpackage."""
     pass


Test patch:
diff --git a/plasmapy/physics/transport/tests/test_collisions.py b/plasmapy/physics/transport/tests/test_collisions.py
index 5a2ceaa6..53c0721e 100644
--- a/plasmapy/physics/transport/tests/test_collisions.py
+++ b/plasmapy/physics/transport/tests/test_collisions.py
@@ -3,7 +3,7 @@
 from astropy import units as u
 
 from plasmapy.physics.transport import (Coulomb_logarithm,
-                                        b_perp,
+                                        impact_parameter_perp,
                                         impact_parameter,
                                         collision_frequency,
                                         mean_free_path,
@@ -38,6 +38,12 @@ def setup_class(self):
         self.gms6 = 3.635342040477818
         self.gms6_negative = 0.030720859361047514
 
+    def test_symmetry(self):
+        lnLambda = Coulomb_logarithm(self.temperature1, self.density2, self.particles)
+        lnLambdaRev = Coulomb_logarithm(self.temperature1, self.density2, self.particles[::-1])
+        assert lnLambda == lnLambdaRev
+
+
     def test_Chen_Q_machine(self):
         """
         Tests whether Coulomb logarithm gives value consistent with
@@ -60,6 +66,7 @@ def test_Chen_Q_machine(self):
                   f"{lnLambdaChen} and not {lnLambda}.")
         assert testTrue, errStr
 
+
     def test_Chen_lab(self):
         """
         Tests whether Coulomb logarithm gives value consistent with
@@ -176,20 +183,13 @@ def test_GMS1_negative(self):
         Murillo, and Schlanges PRE (2002). This checks for when
         a negative (invalid) Coulomb logarithm is returned.
         """
-        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+        with pytest.raises(exceptions.PhysicsError, match="< 0"):
             methodVal = Coulomb_logarithm(self.temperature2,
                                           self.density2,
                                           self.particles,
                                           z_mean=np.nan * u.dimensionless_unscaled,
                                           V=np.nan * u.m / u.s,
                                           method="GMS-1")
-        testTrue = np.isclose(methodVal,
-                              self.gms1_negative,
-                              rtol=1e-15,
-                              atol=0.0)
-        errStr = (f"Coulomb logarithm for GMS-1 should be "
-                  f"{self.gms1_negative} and not {methodVal}.")
-        assert testTrue, errStr
 
     def test_GMS2(self):
         """
@@ -217,20 +217,13 @@ def test_GMS2_negative(self):
         Murillo, and Schlanges PRE (2002). This checks for when
         a negative (invalid) Coulomb logarithm is returned.
         """
-        with pytest.warns(exceptions.PhysicsWarning, match="strong coupling effects"):
+        with pytest.raises(exceptions.PhysicsError, match="< 0"):
             methodVal = Coulomb_logarithm(self.temperature2,
                                           self.density2,
                                           self.particles,
                                           z_mean=self.z_mean,
                                           V=np.nan * u.m / u.s,
                                           method="GMS-2")
-        testTrue = np.isclose(methodVal,
-                              self.gms2_negative,
-                              rtol=1e-15,
-                              atol=0.0)
-        errStr = (f"Coulomb logarithm for GMS-2 should be "
-                  f"{self.gms2_negative} and not {methodVal}.")
-        assert testTrue, errStr
 
     def test_GMS3(self):
         """
@@ -473,7 +466,7 @@ def test_invalid_particle_error(self):
     particles = ('e', 'p')
 
 
-class Test_b_perp:
+class Test_impact_parameter_perp:
     @classmethod
     def setup_class(self):
         """initializing parameters for tests """
@@ -482,19 +475,24 @@ def setup_class(self):
         self.V = 1e4 * u.km / u.s
         self.True1 = 7.200146594293746e-10
 
+    def test_symmetry(self):
+        result = impact_parameter_perp(self.T, self.particles)
+        resultRev = impact_parameter_perp(self.T, self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
         """
-        methodVal = b_perp(self.T,
-                           self.particles,
-                           V=np.nan * u.m / u.s)
+        methodVal = impact_parameter_perp(self.T,
+                                          self.particles,
+                                          V=np.nan * u.m / u.s)
         testTrue = np.isclose(self.True1,
                               methodVal.si.value,
                               rtol=1e-1,
                               atol=0.0)
         errStr = ("Distance of closest approach for 90 degree Coulomb "
-                  f"collision, b_perp, should be {self.True1} and "
+                  f"collision, impact_parameter_perp, should be {self.True1} and "
                   f"not {methodVal}.")
         assert testTrue, errStr
 
@@ -504,14 +502,14 @@ def test_fail1(self):
         value comparison by some quantity close to numerical error.
         """
         fail1 = self.True1 + 1e-15
-        methodVal = b_perp(self.T,
-                           self.particles,
-                           V=np.nan * u.m / u.s)
+        methodVal = impact_parameter_perp(self.T,
+                                          self.particles,
+                                          V=np.nan * u.m / u.s)
         testTrue = not np.isclose(methodVal.si.value,
                                   fail1,
                                   rtol=1e-16,
                                   atol=0.0)
-        errStr = (f"b_perp value test gives {methodVal} and "
+        errStr = (f"impact_parameter_perp value test gives {methodVal} and "
                   f"should not be equal to {fail1}.")
         assert testTrue, errStr
 
@@ -532,6 +530,11 @@ def setup_class(self):
         self.V = 1e4 * u.km / u.s
         self.True1 = np.array([7.200146594293746e-10, 2.3507660003984624e-08])
 
+    def test_symmetry(self):
+        result = impact_parameter(self.T, self.n_e, self.particles)
+        resultRev = impact_parameter(self.T, self.n_e,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -601,6 +604,11 @@ def setup_class(self):
         self.True_protons = 44450104815.91857
         self.True_zmean = 1346828153985.4646
 
+    def test_symmetry(self):
+        result = collision_frequency(self.T, self.n, self.particles)
+        resultRev = collision_frequency(self.T, self.n, self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -710,6 +718,11 @@ def setup_class(self):
         self.V = 1e4 * u.km / u.s
         self.True1 = 4.4047571877932046e-07
 
+    def test_symmetry(self):
+        result = mean_free_path(self.T, self.n_e, self.particles)
+        resultRev = mean_free_path(self.T, self.n_e,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -763,6 +776,11 @@ def setup_class(self):
         self.True1 = 1.2665402649805445e-3
         self.True_zmean = 0.00020264644239688712
 
+    def test_symmetry(self):
+        result = Spitzer_resistivity(self.T, self.n, self.particles)
+        resultRev = Spitzer_resistivity(self.T, self.n,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -830,6 +848,11 @@ def setup_class(self):
         self.True1 = 0.13066090887074902
         self.True_zmean = 0.32665227217687254
 
+    def test_symmetry(self):
+        result = mobility(self.T, self.n_e, self.particles)
+        resultRev = mobility(self.T, self.n_e,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -900,6 +923,11 @@ def setup_class(self):
         self.V = 1e4 * u.km / u.s
         self.True1 = 440.4757187793204
 
+    def test_symmetry(self):
+        result = Knudsen_number(self.length, self.T, self.n_e, self.particles)
+        resultRev = Knudsen_number(self.length, self.T, self.n_e,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.
@@ -956,6 +984,11 @@ def setup_class(self):
         self.True_zmean = 10.689750083758698
         self.True_quantum = 0.3334662805238162
 
+    def test_symmetry(self):
+        result = coupling_parameter(self.T, self.n_e, self.particles)
+        resultRev = coupling_parameter(self.T, self.n_e,  self.particles[::-1])
+        assert result == resultRev
+
     def test_known1(self):
         """
         Test for known value.