docs: def fast_lindblad_terms(gamma, unfolding, matrix_form=False, file_name=None,
                          return_code=False):
      r"""Return a fast function that returns the Lindblad terms.
  
      We test a basic two-level system.
  
      >>> import numpy as np
      >>> Ne = 2
      >>> gamma21 = 2*np.pi*6e6
      >>> gamma = np.array([[0.0, -gamma21],
      ...                   [gamma21, 0.0]])
      >>> rhos = np.array([[0.6, 3+2j],
      ...                  [3-2j, 0.4]])
  
      An map to unfold the density matrix.
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      We obtain a function to calculate Lindblad terms.
      >>> lindblad_terms = fast_lindblad_terms(gamma, unfolding)
  
      Apply this to a density matrix.
      >>> rhos = np.array([[0.6, 3+2j],
      ...                  [3-2j, 0.4]])
      >>> rhosv = unfolding(rhos)
      >>> rhs_lindblad = lindblad_terms(rhosv)
      >>> print(rhs_lindblad)
      [-15079644.7372 -56548667.7646  37699111.8431]
  
      """
      Ne = unfolding.Ne
      Nrho = unfolding.Nrho
      Mu = unfolding.Mu
  
      # We establish the arguments of the output function.
      if True:
          code = ""
          code += "def lindblad_terms("
          if not matrix_form: code += "rho, "
          if code[-2:] == ", ": code = code[:-2]
          code += "):\n"
      # 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"
  
      for a in range(Ne):
          for b in range(a):
              # The first term is of the from
              # gamma_ab * rho_aa |b><b|
              if not (unfolding.normalized and b == 0):
                  coef = gamma[a, b]
                  if unfolding.real:
                      mu = Mu(1, b, b)
                      nu = Mu(1, a, a)
                  else:
                      mu = Mu(0, b, b)
                      nu = Mu(0, a, a)
                  code += term_code(mu, nu, coef, matrix_form, False)
  
              # The second term is of the form
              #  sum_j -gamma_ab/2 rho_aj |a><j|
              # for a lower triangular unfolding, this j runs from 1 to a.
              for j in range(a):
                  coef = -gamma[a, b]*0.5
                  if unfolding.real:
                      mur = Mu(1, a, j)
                      code += term_code(mur, mur, coef, matrix_form, False)
                      mui = Mu(-1, a, j)
                      code += term_code(mui, mui, coef, matrix_form, False)
                  else:
                      mu = Mu(0, a, j)
                      code += term_code(mu, mu, coef, matrix_form, False)
  
              # The third term is of the form
              #  - sum_i 1/2 rho_ia |i><a|
              # for a lower triangular unfolding, this i runs from a to Ne.
              for i in range(a+1, Ne):
                  coef = -gamma[a, b]*0.5
                  if unfolding.real:
                      mur = Mu(1, i, a)
                      code += term_code(mur, mur, coef, matrix_form, False)
                      mui = Mu(-1, i, a)
                      code += term_code(mui, mui, coef, matrix_form, False)
                  else:
                      mu = Mu(0, i, a)
                      code += term_code(mu, mu, coef, matrix_form, False)
  
              # We missed one term in each of the previous fors, that together
              # correspond to
              # -gamma_ab * rho_aa |a><a|
              coef = -gamma[a, b]
              if unfolding.real:
                  mu = Mu(1, a, a)
              else:
                  mu = Mu(0, a, a)
              code += term_code(mu, mu, coef, matrix_form, False)
  
      # We finish the code.
      if True:
          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()
  
          lindblad_terms = code
          if not return_code:
              exec lindblad_terms
      return lindblad_terms