scambio.special

  1r"""Special functions."""
  2
  3import numpy as np
  4from collections.abc import Sequence
  5from functools import lru_cache
  6from mpmath import fp
  7from numbers import Real, Integral
  8from scambio.constant import Constant
  9
 10
 11def rectangle(
 12    x: Sequence[Sequence[Real, ...], ...],
 13    *,
 14    a: Sequence[Real, ...],
 15    b: Sequence[Real, ...],
 16) -> np.ndarray[float]:
 17    r"""Rectangular function.
 18
 19    Args:
 20        x: Independent variable.
 21        a: Rectangle onset.
 22        b: Rectangle offset.
 23
 24    Returns:
 25        Dependent variable.
 26    """
 27    x = np.asarray(x, dtype=float).reshape(1, -1)
 28    a = np.asarray(a, dtype=float).reshape(-1, 1)
 29    b = np.asarray(b, dtype=float).reshape(-1, 1)
 30    y = np.logical_and(a <= x, x <= b).astype(float).squeeze()
 31    return y
 32
 33
 34def bose_einstein_prim(
 35    order: Integral, energy: Real, *, temperature: Real, potential: Real
 36) -> complex | float:
 37    r"""Bose-Einstein primitive.
 38
 39    Note:
 40        See [`ibei` documentation](https://ibei.readthedocs.io/en/stable/)
 41        for details.
 42
 43    Args:
 44        order: Order (n).
 45        energy: Photon energy [$\mathrm{eV}$].
 46        temperature: Temperature [$\mathrm{K}$].
 47        potential: Chemical potential [$\mathrm{eV}$].
 48
 49    Returns:
 50        Bose-Einstein primitive value
 51        [$\mathrm{eV^{n - 2} \, s^{-1} \, cm^{-2}}$].
 52    """
 53    if not isinstance(order, Integral) or not order >= 2:
 54        print("BE order must be an integer greater than or equal to 2!")
 55        return fp.nan
 56
 57    if not isinstance(energy, Real) or not energy >= 0.0:
 58        print("Photon energy must be a nonnegative real number!")
 59        return fp.nan
 60
 61    if not isinstance(temperature, Real) or not temperature > 0.0:
 62        print("Temperature must be a positive real number!")
 63        return fp.nan
 64
 65    if not isinstance(potential, Real):
 66        print("Chemical potential must be a real number!")
 67        return fp.nan
 68
 69    if np.isnan(potential):
 70        return fp.nan
 71
 72    if energy == fp.inf:
 73        return 0.0
 74
 75    factor = -2.0 * fp.pi / Constant.c**2 / Constant.h**3 * fp.factorial(order)
 76    if energy == potential == 0.0:
 77        return (
 78            factor
 79            * (Constant.k * temperature) ** (order + 1)
 80            * fp.polylog(order + 1, 1.0)
 81        )
 82
 83    if energy == potential:
 84        return np.nan
 85
 86    energy_red = (energy - potential) / (Constant.k * temperature)
 87    return factor * sum(
 88        [
 89            fp.factorial(order - i) ** -1
 90            * (Constant.k * temperature) ** (i + 1)
 91            * energy ** (order - i)
 92            * fp.polylog(i + 1, fp.exp(-energy_red))
 93            for i in range(order + 1)
 94        ]
 95    )
 96
 97
 98@lru_cache(maxsize=1000)
 99def bose_einstein_integr(
100    order: Integral, energy: Real, *, temperature: Real, potential: Real
101) -> float:
102    r"""Bose-Einstein integral.
103
104    Note:
105        See [`ibei` documentation](https://ibei.readthedocs.io/en/stable/)
106        for details.
107
108        If the chemical potential is greater than the photon energy, the
109        integrand contains a singularity that makes the integral underfined.
110        However, the Cauchy principal value is finite and equals the real
111        part of the Bose-Einstein primitive, because the imaginary part
112        cancels out.
113
114    Args:
115        order: Order (n).
116        energy: Photon energy [$\mathrm{eV}$].
117        temperature: Temperature [$\mathrm{K}$].
118        potential: Chemical potential [$\mathrm{eV}$].
119
120    Returns:
121        Bose-Einstein integral value
122        [$\mathrm{eV^{n - 2} \, s^{-1} \, cm^{-2}}$].
123    """
124    return -bose_einstein_prim(
125        order, energy, temperature=temperature, potential=potential
126    ).real