docs: def fast_rabi_terms(Ep, epsilonp, rm, xi, theta, unfolding,
                      matrix_form=False, file_name=None, return_code=False):
      r"""Return a fast function that returns the Rabi frequency terms.
  
      We test a basic two-level system.
  
      >>> import numpy as np
      >>> from scipy.constants import physical_constants
      >>> from sympy import Matrix
      >>> from fast.electric_field import electric_field_amplitude_top
      >>> from fast.symbolic import define_laser_variables, polarization_vector
  
      >>> Ne = 2
      >>> Nl = 1
      >>> a0 = physical_constants["Bohr radius"][0]
      >>> rm = [np.array([[0, 0], [a0, 0]]),
      ...       np.array([[0, 0], [0, 0]]),
      ...       np.array([[0, 0], [0, 0]])]
      >>> xi = np.array([[[0, 1], [1, 0]]])
      >>> theta = phase_transformation(Ne, Nl, rm, xi)
  
      We define symbolic variables to be used as token arguments.
      >>> Eps = [electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")]
      >>> Ep, omega_laser = define_laser_variables(Nl)
      >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]
  
      An map to unfold the density matrix.
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      We obtain a function to calculate Rabi frequency terms.
      >>> rabi_terms = fast_rabi_terms(Ep, epsilonps, rm, xi, theta, unfolding)
  
      Apply this to a density matrix.
      >>> rhos = np.array([[0.6, 3+2j],
      ...                  [3-2j, 0.4]])
      >>> rhosv = unfolding(rhos)
      >>> rhs_rabi = rabi_terms(rhosv, Eps)
      >>> print(rhs_rabi)
      [-55680831.474          0.       2784041.5737]
  
      """
      if not unfolding.lower_triangular:
          mes = "It is very inefficient to solve using all components of the "
          mes += "density matrix. Better set lower_triangular=True in Unfolding."
          raise NotImplementedError(mes)
      if matrix_form and (not unfolding.real) and (unfolding.lower_triangular):
          mes = "It is not possible to express the equations in matrix form "
          mes += "for complex lower triangular components only."
          raise ValueError(mes)
      Nl = len(Ep)
      Ne = unfolding.Ne
      # We determine which arguments are constants.
      if True:
          try:
              Ep = np.array([complex(Ep[l]) for l in range(Nl)])
              variable_Ep = False
          except:
              variable_Ep = True
  
          try:
              epsilonp = [np.array([complex(epsilonp[l][i]) for i in range(3)])
                          for l in range(Nl)]
              variable_epsilonp = False
          except:
              variable_epsilonp = True
      # We unpack variables.
      if True:
          Nrho = unfolding.Nrho
          # Mu = unfolding.Mu
          IJ = unfolding.IJ
          normalized = unfolding.normalized
          lower_triangular = unfolding.lower_triangular
  
          # The conjugate stuff.
          Em = [ii.conjugate() for ii in Ep]
          if variable_epsilonp:
              epsilonm = Vector3D(epsilonp.args[0].conjugate())
          else:
              epsilonm = [ii.conjugate() for ii in epsilonp]
          # We convert rm to a numpy array
          rm = np.array([[[complex(rm[k][i, j])
                         for j in range(Ne)] for i in range(Ne)]
                         for k in range(3)])
          rp = np.array([rm[ii].T.conjugate() for ii in range(3)])
          rm_aux = Vector3D(IndexedBase("rm", shape=(Ne, Ne)))
          rp_aux = Vector3D(IndexedBase("rp", shape=(Ne, Ne)))
  
          # We define needed matrices.
          rho = define_density_matrix(Ne, explicitly_hermitian=lower_triangular,
                                      normalized=normalized)
  
          if variable_epsilonp:
              Omega = [[sum([xi[l, i, j] *
                            (Ep[l]*dot(epsilonp[l], rm_aux[i, j]) *
                            delta_greater(i, j) +
                             Em[l]*dot(epsilonm[l], rp_aux[i, j]) *
                            delta_lesser(i, j))
                             for l in range(Nl)])
                       for j in range(Ne)] for i in range(Ne)]
          else:
              Omega = [[sum([xi[l, i, j] *
                            (Ep[l]*dot(epsilonp[l], rm[:, i, j]) +
                             Em[l]*dot(epsilonm[l], rp[:, i, j]))
                             for l in range(Nl)])
                       for j in range(Ne)] for i in range(Ne)]
      # We establish the arguments of the output function.
      if True:
          code = ""
          code += "def rabi_terms("
          if not matrix_form: code += "rho, "
          if variable_Ep: code += "Ep, "
          if variable_epsilonp: code += "epsilonp, "
          if code[-2:] == ", ":
              code = code[:-2]
          code += "):\n"
          code += '    r"""A fast calculation of the Rabi terms."""\n'
      # We initialize the output and auxiliaries.
      if True:
          # We introduce the factor that multiplies all terms.
          if unfolding.real:
              code += "    fact = "+str(e_num/hbar_num)+"\n\n"
          else:
              code += "    fact = "+str(1j*e_num/hbar_num)+"\n\n"
          if variable_epsilonp:
              # We put rm and rp into the code
              np.set_printoptions(threshold=np.nan)
              code += "    rm = np."+rm.__repr__()+"\n\n"
              code += "    rp = np."+rp.__repr__()+"\n\n"
              code += "    def dot(epsilon, rij):\n"
              code += "        return epsilon[0]*rij[0]"
              code += " + epsilon[1]*rij[1]"
              code += " + epsilon[2]*rij[2]\n\n"
              np.set_printoptions(threshold=1000)
  
          if matrix_form:
              code += "    A = np.zeros(("+str(Nrho)+", "+str(Nrho)
              if not unfolding.real:
                  code += "), complex)\n\n"
              else:
                  code += "))\n\n"
              if unfolding.normalized:
                  code += "    b = np.zeros(("+str(Nrho)
                  if not unfolding.real:
                      code += "), complex)\n\n"
                  else:
                      code += "))\n\n"
          else:
              code += "    rhs = np.zeros(("+str(Nrho)
              if not unfolding.real:
                  code += "), complex)\n\n"
              else:
                  code += "))\n\n"
  
      # We write code for the linear terms.
      rho00_terms = []
      for mu in range(Nrho):
          s, i, j = IJ(mu)
          for l in range(Nl):
              for k in range(Ne):
                  if xi[l, i, k] == 1:
                      # There is a I* Omega_l,i,k * rho_k,j term.
                      u = k; v = j
                      args = (0.5*Omega[i][k], rho[k, j], s, i, j, k, u, v,
                              unfolding, matrix_form)
                      term_list = linear_get_coefficients(*args)
                      for term in term_list:
                          code += term_code(*term)
  
                      # We keep note that there was a term with rho00.
                      if k == 0 and j == 0:
                          rho00_terms += [args]
  
                  if xi[l, k, j] == 1:
                      # There is a -I * Omega_l,k,j * rho_i,k term.
                      u = i; v = k
                      args = (-0.5*Omega[k][j], rho[i, k], s, i, j, k, u, v,
                              unfolding, matrix_form)
                      term_list = linear_get_coefficients(*args)
                      for term in term_list:
                          code += term_code(*term)
                      # We keep note that there was a term with rho00.
                      if i == 0 and k == 0:
                          rho00_terms += [args]
  
      # We write code for the independent terms.
      if unfolding.normalized:
          code += "\n    # Independent terms:\n"
          for term in rho00_terms:
              coef_list = independent_get_coefficients(*term)
              for coef in coef_list:
                  code += term_code(*coef, linear=False)
      # We finish the code.
      if True:
          if matrix_form:
              if unfolding.normalized:
                  code += "    A *= fact\n"
                  code += "    b *= fact\n"
                  code += "    return A, b\n"
              else:
                  code += "    A *= fact\n"
                  code += "    return A\n"
          else:
              code += "    rhs *= fact\n"
              code += "    return rhs\n"
      # We write the code to file if provided, and execute it.
      if True:
          if file_name is not None:
              f = file(file_name+".py", "w")
              f.write(code)
              f.close()
  
          rabi_terms = code
          if not return_code:
              exec rabi_terms
      return rabi_terms