code: def thermal_distribution(
      x, amplitude=1.0, center=0.0, kt=1.0, form="bose"
  ):
      """Return a thermal distribution function.
  
      The variable `form` defines the kind of distribution:
  
      - ``form='bose'`` (default) is the Bose-Einstein distribution:
  
          thermal_distribution(x, amplitude=1.0, center=0.0, kt=1.0) =
              1/(amplitude*exp((x - center)/kt) - 1)
  
      - ``form='maxwell'`` is the Maxwell-Boltzmann distribution:
          thermal_distribution(x, amplitude=1.0, center=0.0, kt=1.0) =
              1/(amplitude*exp((x - center)/kt))
  
      - ``form='fermi'`` is the Fermi-Dirac distribution:
          thermal_distribution(x, amplitude=1.0, center=0.0, kt=1.0) =
              1/(amplitude*exp((x - center)/kt) + 1)
  
      Notes
      -----
      - `kt` should be defined in the same units as `x`. The Boltzmann
        constant is ``kB = 8.617e-5 eV/K``.
      - set ``kt<0`` to implement the energy loss convention common in
        scattering research.
  
      For more information, see:
      http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/disfcn.html
  
      """
      form = form.lower()
      if form.startswith("bose"):
          offset = -1
      elif form.startswith("maxwell"):
          offset = 0
      elif form.startswith("fermi"):
          offset = 1
      else:
          msg = (
              f"Invalid value ('{form}') for argument 'form'; should be one "
              "of 'maxwell', 'fermi', or 'bose'."
          )
          raise ValueError(msg)
  
      return real(
          1
          / (
              amplitude * exp((x - center) / not_zero(kt))
              + offset
              + tiny * 1j
          )
      )
  
imports: from numpy import arctan, copysign, cos, exp, isclose, isnan, log, log1p, maximum, minimum, pi, real, sin, sqrt, where
  from scipy.special import betaln
  from scipy.special import erf, erfc
  from scipy.special import gamma
  from scipy.special import loggamma
  from scipy.special import wofz
label_desc: Use image node_id 23 for calling a global function with example usage: thermal_distribution(x, amplitude, center, kt, form) and returns: real