code: def f0(E, fermi, T):
      
      return 1. / (1. + np.exp((E - fermi) / (_cd("Boltzmann constant in eV/K") * T)))
docs: Returns the equilibrium fermi-dirac.
      Args:
          E (float): energy in eV
          fermi (float): the fermi level in eV
          T (float): the temperature in kelvin