docs: 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