code: def _frank_help(alpha, tau):
          
  
          def debye(t):
              return t / (np.exp(t) - 1)
  
          debye_value = integrate.quad(debye, EPSILON, alpha)[0] / alpha
          return 4 * (debye_value - 1) / alpha + 1 - tau
docs: Compute first order debye function to estimate theta.