Repository: jrsmith3/ibei

Issue:
Add functionality to compute full Bose-Einstein integral
# Overview
Functionality to complete the full Bose-Einstein integral would be a good compliment to the upper-incomplete functionality.


# Related issues
* Depends on #68.

Discussion:

Patch:
diff --git a/CHANGELOG b/CHANGELOG
index e568d1c..158497f 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -7,6 +7,7 @@
 * 50: Leverage `tox` to run tests.
 * 49: Refactor tests to use `pytest`.
 * 48: Replace `conda` with `hatch`.
+* 43: Add method to compute full Bose-Einstein integral.
 * 39: Update for basic compatibility with Python 3.
 
 
diff --git a/src/ibei/models.py b/src/ibei/models.py
index 3c0d0b3..a8e2cb6 100644
--- a/src/ibei/models.py
+++ b/src/ibei/models.py
@@ -151,34 +151,112 @@ class BEI():
         ----------
         .. [1] :cite:`10.1016/j.sse.2006.06.017`
         """
-        kT = self.temperature * astropy.constants.k_B
-
-        reduced_energy_bound = self.energy_bound / kT
-        reduced_chemical_potential = self.chemical_potential / kT
-
-        prefactor = (2 * np.pi * np.math.factorial(self.order) * kT**(self.order + 1)) / \
-            (astropy.constants.h**3 * astropy.constants.c**2)
-
-        expt = (reduced_chemical_potential - reduced_energy_bound).decompose()
+        expt = self.reduced_chemical_potential - self.reduced_energy_bound
         real_arg = np.exp(expt.value)
 
-        if reduced_chemical_potential == 0 and reduced_energy_bound == 0:
+        if self.reduced_chemical_potential == 0 and self.reduced_energy_bound == 0:
             # Specify this condition just to skip the next condition.
             term = float(mpmath.polylog(self.order + 1, real_arg))
-            bei = term * prefactor
+            bei = term * self.prefactor * np.math.factorial(self.order)
 
-        elif reduced_chemical_potential >= reduced_energy_bound:
-            bei = 0 * prefactor
+        elif self.reduced_chemical_potential >= self.reduced_energy_bound:
+            bei = 0 * self.prefactor * np.math.factorial(self.order)
 
         else:
             summand = 0
             for indx in range(1, self.order + 2):
                 index = self.order - indx + 1
 
-                term = reduced_energy_bound**index * float(mpmath.polylog(indx, real_arg)) / np.math.factorial(index)
+                term = self.reduced_energy_bound**index * float(mpmath.polylog(indx, real_arg)) / np.math.factorial(index)
 
                 summand += term
 
-            bei = prefactor * summand
+            bei = self.prefactor * summand * np.math.factorial(self.order)
 
         return bei
+
+
+    def full(self) -> astropy.units.Quantity:
+        """
+        Full Bose-Einstein integral.
+
+        Returns
+        -------
+        astropy.units.Quantity
+            Value of the Bose-Einstein integral.
+        """
+        expt = self.reduced_chemical_potential
+        real_arg = np.exp(expt.value)
+
+        if self.reduced_chemical_potential > 0:
+            bei = 0 * self.prefactor
+
+        else:
+            bei = self.prefactor * float(mpmath.gamma(self.order + 1)) * float(mpmath.polylog(self.order + 1, real_arg))
+
+        return bei
+
+
+    def photon_flux(self) -> astropy.units.Quantity[astropy.units.Unit("1/(m2 s)")]:
+        """
+        Number of photons radiated per unit time per unit area.
+
+        Notes
+        -----
+        This convenience method is a special case of `BEI.full`. This method
+        assumes the value of `order` is 2.
+        """
+        flux = (4 * np.pi * float(mpmath.zeta(3)) * self.kT**3) / \
+            (astropy.constants.h**3 * astropy.constants.c**2)
+
+        return flux.to("1/(m2 s)")
+
+
+    def radiant_power_flux(self) -> astropy.units.Quantity[astropy.units.Unit("J/(m2 s)")]:
+        """
+        Energy radiated per unit time per unit area.
+
+        Notes
+        -----
+        This convenience method implements the Stefan-Boltzmann law and is a
+        special case of the `BEI.full`. This method assumes the value of
+        `order` is 3.
+        """
+        power_flux = astropy.constants.sigma_sb * self.temperature**4
+
+        return power_flux.to("J/(m2 s)")
+
+
+    @property
+    def kT(self) -> astropy.units.Quantity[astropy.units.dimensionless_unscaled]:
+        """
+        Product of Boltzmann's constant and temperature.
+        """
+        return (self.temperature * astropy.constants.k_B).decompose()
+
+    @property
+    def reduced_energy_bound(self) -> astropy.units.Quantity[astropy.units.dimensionless_unscaled]:
+        """
+        Dimensionless energy bound in units of kT.
+        """
+        return (self.energy_bound / self.kT).decompose()
+
+    @property
+    def reduced_chemical_potential(self) -> astropy.units.Quantity[astropy.units.dimensionless_unscaled]:
+        """
+        Dimensionless chemical potential in units of kT.
+        """
+        return (self.chemical_potential / self.kT).decompose()
+
+    @property
+    def prefactor(self) -> astropy.units.Quantity:
+        """
+        Factor preceding Bose-Einstein integral integration result.
+
+        Note
+        ----
+        For `order=3`, this factor is NOT equal to the Stefan-Boltzmann
+        constant. This factor is missing the Riemann-Zeta term.
+        """
+        return (2 * np.pi * self.kT**(self.order + 1)) / \
+            (astropy.constants.h**3 * astropy.constants.c**2)


Test patch:
diff --git a/test/test_BEI.py b/test/test_BEI.py
index 189842b..12c013b 100644
--- a/test/test_BEI.py
+++ b/test/test_BEI.py
@@ -216,6 +216,76 @@ class TestIssues():
                 "upper",
                 astropy.units.Quantity(10549122.240303338, "1/(m2 s)"),
             ),
+            (
+                {
+                    "order": 3,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "full",
+                astropy.units.Quantity(459.30032795, "W/m2"),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "photon_flux",
+                astropy.units.Quantity(4.1052443203614687e+22, "1/(m2 s)"),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "radiant_power_flux",
+                astropy.units.Quantity(459.30032795, "W/m2"),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "kT",
+                astropy.units.Quantity(4.141947e-21, "J"),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "reduced_energy_bound",
+                astropy.units.Quantity(42.54989978, astropy.units.dimensionless_unscaled),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 2.1,
+                    "temperature": 300,
+                    "chemical_potential": 1.1,
+                },
+                "reduced_chemical_potential",
+                astropy.units.Quantity(42.54989978, astropy.units.dimensionless_unscaled),
+            ),
+            (
+                {
+                    "order": 2,
+                    "energy_bound": 1.1,
+                    "temperature": 300,
+                    "chemical_potential": 0,
+                },
+                "prefactor",
+                astropy.units.Quantity(1.70759151e+22, "1/(m2 s)"),
+            ),
         ]
     )
 def test_methods_regression(args, method_under_test, expected_output):
@@ -227,7 +297,12 @@ def test_methods_regression(args, method_under_test, expected_output):
     This test tests each of the methods of a `BEI` instance at least once.
     """
     bei = ibei.models.BEI(**args)
-    output = getattr(bei, method_under_test)()
+    method = getattr(bei, method_under_test)
+
+    if callable(method):
+        output = method()
+    else:
+        output = method
 
     assert astropy.units.allclose(expected_output, output)
 
@@ -237,7 +312,12 @@ def test_methods_regression(args, method_under_test, expected_output):
             (3, "J/(m2 s)"),
         ]
     )
-@pytest.mark.parametrize("method_under_test", ("upper",))
+@pytest.mark.parametrize("method_under_test", 
+        (
+            "upper",
+            "full",
+        )
+    )
 def test_methods_units(order, expected_unit, method_under_test, valid_constructor_quantity_args):
     """
     Methods' units should match known units for low orders
@@ -258,6 +338,38 @@ def test_methods_units(order, expected_unit, method_under_test, valid_constructo
     assert output.unit.is_equivalent(expected_unit)
 
 
+def test_consistency_upper_and_full_methods(valid_constructor_quantity_args):
+    """
+    When `energy_bound` is 0, `BEI.upper` should equal `BEI.full`.
+    """
+    valid_constructor_quantity_args["energy_bound"] = 0
+
+    bei = ibei.models.BEI(**valid_constructor_quantity_args)
+
+    assert astropy.units.allclose(bei.upper(), bei.full())
+
+
+@pytest.mark.parametrize("order,helper_method_name",[
+            (2, "photon_flux"),
+            (3, "radiant_power_flux"),
+        ]
+    )
+def test_consistency_full_and_helper_methods(order, helper_method_name, valid_constructor_quantity_args):
+    """
+    `BEI.full` should equal `BEI.photon_flux` and `BEI.radiant_power_flux` for order 2 and 3, respectively.
+
+    Notes
+    -----
+    This condition holds as long as `chemical_potential` equals zero.
+    """
+    valid_constructor_quantity_args["order"] = order
+    valid_constructor_quantity_args["chemical_potential"] = 0.
+    bei = ibei.models.BEI(**valid_constructor_quantity_args)
+    output = getattr(bei, helper_method_name)()
+
+    assert astropy.units.allclose(bei.full(), output)
+
+
 # Pytest fixture definitions
 # ==========================
 @pytest.fixture