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)))