code: def debye_integral(y):
          
          # floating point limit is reached around y=155, so values beyond that
          # are set to the limiting value(T-->0, y --> \infty) of
          # 6.4939394 (from wolfram alpha).
          factor = 3. / y ** 3
          if y < 155:
              integral = quadrature(lambda x: x ** 3 / (np.exp(x) - 1.), 0, y)
              return list(integral)[0] * factor
          else:
              return 6.493939 * factor
docs: Debye integral. Eq(5) in  doi.org/10.1016/j.comphy.2003.12.001
  
          Args:
              y (float): debye temperature/T, upper limit
  
          Returns:
              float: unitless