docs: def bethe_find_crystalfield(populations, hopping):
      """Return the orbital energies to have the system populates as
         desired by the given individual populations"""
  
      zero = lambda orb: [bethe_filling_zeroT(-em, tz) - pop \
                          for em, tz, pop in zip(orb, hopping, populations)]
  
      return fsolve(zero, np.zeros(len(populations)))