def hamiltonian(Ep, epsilonp, detuning_knob, rm, omega_level, omega_laser, xi,
                RWA=True, RF=True):
    r"""Return symbolic Hamiltonian.

    >>> from sympy import zeros, pi, pprint, symbols
    >>> Ne = 3
    >>> Nl = 2
    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> epsilonp = [polarization_vector(0, -pi/2, 0, 0, 1) for l in range(Nl)]
    >>> detuning_knob = symbols("delta1 delta2", real=True)

    >>> xi = [zeros(Ne, Ne) for l in range(Nl)]
    >>> coup = [[(1, 0)], [(2, 0)]]
    >>> for l in range(Nl):
    ...     for pair in coup[l]:
    ...         xi[l][pair[0], pair[1]] = 1
    ...         xi[l][pair[1], pair[0]] = 1

    >>> rm = define_r_components(Ne, xi, explicitly_hermitian=True,
    ...                          helicity=True, p=-1)

    >>> rm = helicity_to_cartesian(rm)
    >>> omega_level, omega, gamma = define_frequencies(Ne, True)

    >>> H = hamiltonian(Ep, epsilonp, detuning_knob, rm, omega_level,
    ...                 omega_laser, xi, RWA=True, RF=False)

    >>> print H[1, 0]
    -E_{01}*e*r_{0;21}*exp(-I*t*varpi_1)/2
    >>> print H[2, 0]
    -E_{02}*e*r_{0;31}*exp(-I*t*varpi_2)/2
    >>> print H[2, 2]
    hbar*omega_3

    """
    # We check what RF is.
    if type(RF) == list:
        theta = RF[:]
        RF = True
    elif type(RF) == Matrix:
        theta = [RF[i, 0] for i in range(RF.shape[0])]
        RF = True
    elif RF:
        # theta should be calculate here!
        s = "We are still missing automatic calculation of phase "
        s += "transformations."
        raise ValueError(s)

    if not RWA and RF:
        s = "The rotating frame does not exist without the rotating wave \
approximation, as far as I know."
        raise ValueError(s)

    # We check that the epsilonp is a list of vectors.
    if not isinstance(epsilonp, list):
        raise ValueError("epsilonp must be a list of polarization vectors.")
    if not isinstance(epsilonp[0], Matrix):
        raise ValueError("epsilonp must be a list of polarization vectors.")

    Ne = len(omega_level)
    Nl = len(omega_laser)

    H = zeros(Ne, Ne)
    hbar, e = symbols("hbar e", positive=True)
    t = symbols("t", real=True)

    for i in range(Ne):
        for j in range(Ne):
            rmij = vector_element(rm, i, j)
            rpij = vector_element(rm, j, i).conjugate()
            for l in range(Nl):
                epsilonpl = epsilonp[l]
                epsilonml = epsilonpl.conjugate()

                if RF:
                    Epl = xi[l][i, j]*Ep[l]
                    Epl *= exp(-I*(theta[j]-theta[i]-t*omega_laser[l]))
                    Eml = xi[l][i, j]*Ep[l].conjugate()
                    Eml *= exp(-I*(theta[j]-theta[i]+t*omega_laser[l]))
                else:
                    Epl = Ep[l]*xi[l][i, j]*exp(-I*omega_laser[l]*t)
                    Eml = Epl.conjugate()

                # The E^(+)r^(-) term
                H[i, j] += -e*Epl/2*cartesian_dot_product(epsilonpl, rmij)
                # The E^(-)r^(+) term
                H[i, j] += -e*Eml/2*cartesian_dot_product(epsilonml, rpij)
                if not RWA:
                    # The E^(+)r^(+) term
                    H[i, j] += -e*Epl/2*cartesian_dot_product(epsilonpl, rpij)
                    # The E^(-)r^(-) term
                    H[i, j] += -e*Eml/2*cartesian_dot_product(epsilonml, rmij)

            if i == j:
                if RF:
                    H[i, j] += hbar*(omega_level[i]+diff(theta[i], t))
                else:
                    H[i, j] += hbar*omega_level[i]

    return H