<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def thermal_state(omega_level, T, return_diagonal=False):
 r"""Return a thermal state for a given set of levels. INPUT: - ``omega_level`` - The angular frequencies of each state. - ``T`` - The temperature of the ensemble (in Kelvin). - ``return_diagonal`` - Whether to return only the populations. [0.0834 0.0834 0.0834 0.0834 0.0834 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833] """

    Ne = len(omega_level)
    E = np.array([hbar*omega_level[i] for i in range(Ne)])
    p = np.exp(-E/k_B/T)
    p = p/sum(p)
    if not return_diagonal:
        return np.diag(p)
    return p