docs: def fast_hamiltonian_terms(Ep, epsilonp, detuning_knob,
                             omega_level, rm, xi, theta,
                             unfolding, matrix_form=False, file_name=None,
                             return_code=False):
      r"""Return a fast function that returns the Hamiltonian terms.
  
      We test a basic two-level system.
  
      >>> import numpy as np
      >>> from scipy.constants import physical_constants
      >>> from sympy import Matrix, symbols
      >>> 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]]])
      >>> omega_level = [0, 1.0e9]
      >>> theta = phase_transformation(Ne, Nl, rm, xi)
  
      We define symbolic variables to be used as token arguments.
      >>> Ep, omega_laser = define_laser_variables(Nl)
      >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]
      >>> detuning_knob = [symbols("delta1", real=True)]
  
      An map to unfold the density matrix.
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      We obtain a function to calculate Hamiltonian terms.
      >>> aux = (Ep, epsilonps, detuning_knob, omega_level, rm, xi, theta,
      ...        unfolding, False, None)
      >>> hamiltonian_terms = fast_hamiltonian_terms(*aux)
  
      Apply this to a density matrix.
      >>> rhos = np.array([[0.6, 3+2j],
      ...                  [3-2j, 0.4]])
      >>> rhosv = unfolding(rhos)
  
      We specify values for the variables
      >>> detuning_knobs = [100e6]
      >>> Eps = electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")
      >>> Eps *= np.exp(1j*np.pi)
      >>> Eps = [Eps]
      >>> print(hamiltonian_terms(rhosv, Eps, detuning_knobs))
      [5.5681e+07 2.0000e+08 2.9722e+08]
  
  
      """
      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)
      Nrho = unfolding.Nrho
      # 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
          try:
              detuning_knob = np.array([float(detuning_knob[l])
                                        for l in range(Nl)])
              variable_detuning_knob = False
          except:
              variable_detuning_knob = True
      # We obtain code for the two parts.
      if True:
          if file_name is not None:
              file_name_rabi = file_name+"_rabi"
              file_name_detuning = file_name+"_detuning"
          else:
              file_name_rabi = file_name
              file_name_detuning = file_name
  
          rabi_terms = fast_rabi_terms(Ep, epsilonp, rm, xi, theta, unfolding,
                                       matrix_form=matrix_form,
                                       file_name=file_name_rabi,
                                       return_code=True)
  
          detuning_terms = fast_detuning_terms(detuning_knob, omega_level, xi,
                                               theta, unfolding,
                                               matrix_form=matrix_form,
                                               file_name=file_name_detuning,
                                               return_code=True)
          code = rabi_terms + "\n\n" + detuning_terms + "\n\n"
  
          # If these functions have 0 arguments, we call them only once!
          if not variable_Ep and not variable_epsilonp and matrix_form:
              code += "rabi_terms = rabi_terms()\n\n"
          if not variable_detuning_knob and matrix_form:
              code += "detuning_terms = detuning_terms()\n\n"
      # We establish the arguments of the output function.
      if True:
          code += "def hamiltonian_terms("
          if not matrix_form: code += "rho, "
          if variable_Ep: code += "Ep, "
          if variable_epsilonp: code += "epsilonp, "
          if variable_detuning_knob: code += "detuning_knob, "
          code += "rabi_terms=rabi_terms, detuning_terms=detuning_terms"
          # if code[-2:] == ", ": code = code[:-2]
          code += "):\n"
  
          code += '    r"""A fast calculation of the hamiltonian terms."""\n'
          # if not variable_Ep and not varia
      # We initialize the output and auxiliaries.
      if True:
          # We introduce the factor that multiplies all terms.
          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 call the Rabi terms.
      if True:
          if not variable_Ep and not variable_epsilonp and matrix_form:
              aux_code = "rabi_terms\n"
          else:
              aux_code = "rabi_terms("
              if not matrix_form: aux_code += "rho, "
              if variable_Ep: aux_code += "Ep, "
              if variable_epsilonp: aux_code += "epsilonp, "
              if aux_code[-2:] == ", ": aux_code = aux_code[:-2]
              aux_code += ")\n"
  
          if matrix_form:
              if unfolding.normalized:
                  code += "    aux = " + aux_code
                  code += "    A += aux[0]\n"
                  code += "    b += aux[1]\n"
              else:
                  code += "    A = " + aux_code
          else:
              code += "    rhs = " + aux_code
      # We call the detuning terms.
      if True:
          if not variable_detuning_knob and matrix_form:
              aux_code = "detuning_terms\n"
          else:
              aux_code = "detuning_terms("
              if not matrix_form: aux_code += "rho, "
              if variable_detuning_knob: aux_code += "detuning_knob, "
              if aux_code[-2:] == ", ": aux_code = aux_code[:-2]
              aux_code += ")\n"
  
          if matrix_form:
              if unfolding.normalized:
                  code += "    aux = " + aux_code
                  code += "    A += aux[0]\n"
                  code += "    b += aux[1]\n"
              else:
                  code += "    A += " + aux_code
          else:
              code += "    rhs += " + aux_code
      # We finish the code.
      if True:
          # code = rabi_code + "\n\n" + code
          if matrix_form:
              if unfolding.normalized:
                  code += "    return A, b\n"
              else:
                  code += "    return A\n"
          else:
              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()
  
          hamiltonian_terms = code
          if not return_code:
              exec hamiltonian_terms
      return hamiltonian_terms