docs: def fast_hamiltonian(Ep, epsilonp, detuning_knob, rm, omega_level, xi, theta,
                       file_name=None):
      r"""Return a fast function that returns a Hamiltonian as an array.
  
      INPUT:
  
      -  ``Ep`` - A list with the electric field amplitudes (real or complex).
      -  ``epsilonp`` - A list of the polarization vectors of the fields \
      (real or complex).
      -  ``detuning_knob`` - A list of the detunings of each field (relative \
      to the transition of lowest energy).
      -  ``rm`` -     The below-diagonal components
          of the position operator in the cartesian basis:
  
          .. math::
              \vec{r}^{(-)}_{i j} = [ x_{ij}, y_{ij}, z_{ij} ]
              \hspace{1cm} \forall \hspace{1cm} 0 < j < i
      -  ``omega_level`` - The angular frequencies of each state.
      -  ``xi`` - An array whose ``xi[l, i, j]`` element is 1 if the \
       transition :math:`|i\rangle \rightarrow |j\rangle`\ is driven by field \
       ``l`` and 0 otherwise.
      -  ``theta`` - A list of symbolic expressions representing a phase \
      transformation.
      -  ``file_name`` - A string indicating a file to save the function's \
      code.
  
      If the arguments Ep, epsilonp, and detuning_knob are symbolic amounts, \
      the returned function will accept numeric values of Ep, epsilonp, and \
      detuning_knob as arguments.
  
      All quantities should be in SI units.
  
      EXAMPLES:
      We build an example using states coupled like this:
  
       --- |4>        --- |5>          --- |6>
        ^              ^                ^
        |              |                |
        |    --- |2>   |     --- |3>    |
      2 |     ^      2 |      ^         | 2
        |   1 |        |    1 |         |
        |     |        |      |         |
      ------------------------------------- |1>
  
      With the numbers on kets labeling states and the plain numbers labeling
      fields.
  
      The number of states and fields:
      >>> Ne = 6
      >>> Nl = 2
  
      We invent some energy levels:
      >>> omega_level = np.array([0.0, 100.0, 100.0, 200.0, 200.0, 300.0])
      >>> omega_level = omega_level*1e6*2*np.pi
  
      We build the symbol xi, that chooses which laser couples which
      transition.
      >>> xi = np.zeros((Nl, Ne, Ne))
      >>> coup = [[(1, 0), (2, 0)], [(3, 0), (4, 0), (5, 0)]]
      >>> for l in range(Nl):
      ...     for pair in coup[l]:
      ...         xi[l, pair[0], pair[1]] = 1.0
      ...         xi[l, pair[1], pair[0]] = 1.0
  
      We invent some electric dipole matrix elements:
      >>> from scipy.constants import physical_constants
      >>> a0 = physical_constants["Bohr radius"][0]
      >>> rm = np.zeros((3, Ne, Ne))
      >>> for l in range(Nl):
      ...     for i in range(Ne):
      ...         for j in range(i):
      ...             if xi[l, i, j] != 0:
      ...                 rm[2, i, j] = float(i)*a0
  
      The phase transformation:
      >>> theta = phase_transformation(Ne, Nl, rm, xi)
  
      We define the possible arguments:
      >>> from sympy import symbols, pi
      >>> from fast.symbolic import polarization_vector
      >>> detuning_knob = symbols("delta1 delta2")
      >>> detuning_knob_vals = np.array([-1.0, 3.0])*1e6*2*np.pi
  
      >>> Ep, omega_laser = define_laser_variables(Nl)
      >>> Ep_vals = [1e2, 1e2]
  
      >>> alpha = symbols("alpha")
      >>> epsilon = polarization_vector(0, pi/2, alpha, 0, 1)
      >>> epsilonp = [epsilon, epsilon]
      >>> epsilonp_vals = [[0.0, 0.0, 1.0], [0.0, 0.0, 1.0]]
  
      There are 8 ways to call fast_hamiltonian:
  
      1 .- Get a function of detunings, field amplitudes, polarizations:
      >>> H1 = fast_hamiltonian(Ep, epsilonp, detuning_knob, rm,
      ...                       omega_level, xi, theta)
  
      2 .- Get a function of field amplitudes, polarizations:
      >>> H2 = fast_hamiltonian(Ep, epsilonp, detuning_knob_vals, rm,
      ...                       omega_level, xi, theta)
  
      3 .- Get a function of detunings, polarizations:
      >>> H3 = fast_hamiltonian(Ep_vals, epsilonp, detuning_knob, rm,
      ...                       omega_level, xi, theta)
  
      4 .- Get a function of detunings, field amplitudes:
      >>> H4 = fast_hamiltonian(Ep, epsilonp_vals, detuning_knob, rm,
      ...                       omega_level, xi, theta)
  
      5 .- Get a function of detunings:
      >>> H5 = fast_hamiltonian(Ep_vals, epsilonp_vals, detuning_knob, rm,
      ...                       omega_level, xi, theta)
  
      6 .- Get a function of field amplitudes:
      >>> H6 = fast_hamiltonian(Ep, epsilonp_vals, detuning_knob_vals, rm,
      ...                       omega_level, xi, theta)
  
      7 .- Get a function of polarizations:
      >>> H7 = fast_hamiltonian(Ep_vals, epsilonp, detuning_knob_vals, rm,
      ...                       omega_level, xi, theta)
  
      8 .- Get a function of nothing:
      >>> H8 = fast_hamiltonian(Ep_vals, epsilonp_vals, detuning_knob_vals, rm,
      ...                       omega_level, xi, theta)
  
  
      We test all of these combinations.
      >>> print(H1(Ep_vals, epsilonp_vals, detuning_knob_vals) \
      ...     /hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H2(Ep_vals, epsilonp_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H3(epsilonp_vals, detuning_knob_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H4(Ep_vals, detuning_knob_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H5(detuning_knob_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H6(Ep_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H7(epsilonp_vals)/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      >>> print(H8()/hbar_num/2/np.pi*1e-6)
      [[ 0.    +0.j  0.6398+0.j  1.2795+0.j  1.9193+0.j  2.5591+0.j  3.1989+0.j]
       [ 0.6398+0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.2795+0.j  0.    +0.j  1.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j]
       [ 1.9193+0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j  0.    +0.j]
       [ 2.5591+0.j  0.    +0.j  0.    +0.j  0.    +0.j -3.    +0.j  0.    +0.j]
       [ 3.1989+0.j  0.    +0.j  0.    +0.j  0.    +0.j  0.    +0.j 97.    +0.j]]
  
      """
      # We determine which arguments are constants.
      if True:
          Nl = len(Ep)
          Ne = np.array(rm[0]).shape[0]
          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 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)])
      # We establish the arguments of the output function.
      if True:
          code = ""
          code += "def hamiltonian("
          if variable_Ep: code += "Ep, "
          if variable_epsilonp: code += "epsilonp, "
          if variable_detuning_knob: code += "detuning_knob, "
          if code[-2:] == ", ":
              code = code[:-2]
          code += "):\n"
  
          code += '    r"""A fast calculation of the hamiltonian."""\n'
          code += "    H = np.zeros(("+str(Ne)+", "+str(Ne)+"), complex)\n\n"
      # We get the code for the below-diagonal elements
      # (Rabi frequencies).
      if True:
          code += "    # We calculate the below-diagonal elements.\n"
          for i in range(Ne):
              for j in range(i):
                  for l in range(Nl):
                      if xi[l, i, j] == 1.0:
                          # We get the below-diagonal terms.
                          code += "    H["+str(i)+", "+str(j)+"] = "
                          # We get the code for Ep.
                          if variable_Ep:
                              code += "0.5*Ep["+str(l)+"]"
  
                          else:
                              code += str(0.5*Ep[l])
                          # We get the code for epsilonp dot rm
                          rmij = rm[:, i, j]
                          if variable_epsilonp:
                              code += "*cartesian_dot_product("
                              code += "epsilonp["+str(l)+"],"
                              code += str(list(rmij*e_num))+" )"
                          else:
                              dp = cartesian_dot_product(epsilonp[l], rmij)
                              dp = dp*e_num
                              code += "*("+str(dp)+")"
  
                          code += "\n"
      # We get the code for the above-diagonal elements
      # (Conjugate Rabi frequencies).
      if True:
          code += "\n"
          code += """    # We calculate the above-diagonal elements.\n"""
          code += """    for i in range("""+str(Ne)+"""):\n"""
          code += """        for j in range(i+1, """+str(Ne)+"""):\n"""
          code += """            H[i, j] = H[j, i].conjugate()\n\n"""
      # We get the code for the diagonal elements (detunings).
      if True:
          code += "    # We calculate the diagonal elements.\n"
          # We build the degeneration simplification and is inverse (to avoid
          # large combinatorics).
          aux = define_simplification(omega_level, xi, Nl)
          u, invu, omega_levelu, Neu, xiu = aux
          # For each field we find the smallest transition frequency, and its
          # simplified indices.
          omega_min, iu0, ju0 = find_omega_min(omega_levelu, Neu, Nl, xiu)
          #####################################
          # We get the code to calculate the non degenerate detunings.
          pairs = detunings_indices(Neu, Nl, xiu)
          if not variable_detuning_knob:
              code += "    detuning_knob = np.zeros("+str(Nl)+")\n"
              for l in range(Nl):
                  code += "    detuning_knob["+str(l)+"] = " +\
                      str(detuning_knob[l])+"\n"
  
          code_det = detunings_code(Neu, Nl, pairs, omega_levelu, iu0, ju0)
          code += code_det
          code += "\n"
          #####################################
          # We find the coefficients a_l that multiply omega_laser_l in
          # H_ii = omega_level_iu + theta_iu = \sum_i a_i varpi_i + remainder
          _omega_level, omega, gamma = define_frequencies(Ne)
          _omega_levelu, omega, gamma = define_frequencies(Neu)
          E0, omega_laser = define_laser_variables(Nl)
          # So we build all combinations.
          combs = detunings_combinations(pairs)
          for i in range(Ne):
              _Hii = theta[i] + _omega_levelu[u(i)]
              aux = (_Hii, combs, omega_laser,
                     _omega_levelu, omega_levelu, iu0, ju0)
              assign = detunings_rewrite(*aux)
  
              if assign != "":
                  code += "    H["+str(i)+", "+str(i)+"] = "+assign+"\n"
  
      code += "\n"
      code += """    for i in range("""+str(Ne)+"""):\n"""
      code += """        H[i, i] = H[i, i]*"""+str(hbar_num)+"\n"
  
      code += "    return H\n"
  
      if file_name is not None:
          f = file(file_name, "w")
          f.write(code)
          f.close()
  
      hamiltonian = code
      exec hamiltonian
      return hamiltonian