docs: def detunings_rewrite(expr, combs, omega_laser, symb_omega_levelu,
                        omega_levelu, iu0, ju0):
      r"""Rewrite a symbolic expression in terms of allowed transition detunings.
  
      >>> Ne = 6
      >>> Nl = 2
      >>> omega_level = [0.0, 100.0, 100.0, 200.0, 200.0, 300.0]
      >>> 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
  
      >>> aux = define_simplification(omega_level, xi, Nl)
      >>> u, invu, omega_levelu, Neu, xiu = aux
      >>> omega_min, iu0, ju0 = find_omega_min(omega_levelu, Neu, Nl, xiu)
      >>> pairs = detunings_indices(Neu, Nl, xiu)
      >>> combs = detunings_combinations(pairs)
      >>> symb_omega_levelu, omega, gamma = define_frequencies(Neu)
      >>> E0, omega_laser = define_laser_variables(Nl)
  
      Most times it is possible to express these combinations of optical
      frequencies in terms of allowed transition detunings.
  
      >>> expr = +(omega_laser[0]-(symb_omega_levelu[1]-symb_omega_levelu[0]))
      >>> expr += -(omega_laser[1]-(symb_omega_levelu[3]-symb_omega_levelu[0]))
      >>> expr
      -omega_2 + omega_4 + varpi_1 - varpi_2
      >>> detunings_rewrite(expr, combs, omega_laser, symb_omega_levelu,
      ...                   omega_levelu, iu0, ju0)
      '+delta1_2_1-delta2_4_1'
  
      But some times it is not possible:
  
      >>> expr = +(omega_laser[1]-(symb_omega_levelu[1]-symb_omega_levelu[0]))
      >>> expr += -(omega_laser[0]-(symb_omega_levelu[3]-symb_omega_levelu[0]))
      >>> expr
      -omega_2 + omega_4 - varpi_1 + varpi_2
      >>> detunings_rewrite(expr, combs, omega_laser, symb_omega_levelu,
      ...                   omega_levelu, iu0, ju0)
      '300.000000000000-detuning_knob[0]+detuning_knob[1]'
  
      """
      Nl = len(omega_laser)
      Neu = len(symb_omega_levelu)
      # We find the coefficients a_i of the field frequencies.
      a = [diff(expr, omega_laser[l]) for l in range(Nl)]
  
      # We look for a combination of the detunings obtained with the
      # function detunings_code. For each combination we sum the
      # detunings weighed by a_i.
      success = False
      for comb in combs:
          expr_try = 0
          for l in range(Nl):
              expr_try += a[l]*(omega_laser[l] -
                                symb_omega_levelu[comb[l][0]] +
                                symb_omega_levelu[comb[l][1]])
          if expr-expr_try == 0:
              success = True
              break
  
      assign = ""
      if success:
          for l in range(Nl):
              if a[l] != 0:
                  if a[l] == 1:
                      assign += "+"
                  elif a[l] == -1:
                      assign += "-"
                  assign += "delta"+str(l+1)
                  assign += "_"+str(comb[l][0]+1)
                  assign += "_"+str(comb[l][1]+1)
      else:
          # We get the code for Hii using detuning knobs.
          # We find out the remainder terms.
          _remainder = expr - sum([a[l]*omega_laser[l]
                                   for l in range(Nl)])
          # We find the coefficients of the remainder.
          b = [diff(_remainder, symb_omega_levelu[j]) for j in range(Neu)]
          # We calculate the remainder numerically.
          remainder = sum([b[j]*omega_levelu[j] for j in range(Neu)])
          # We add the contributions from the detuning knobs.
          remainder += sum([a[l]*(omega_levelu[iu0[l]] -
                                  omega_levelu[ju0[l]])
                            for l in range(Nl)])
          assign = str(remainder)
          # We get the code for Hii using detuning knobs.
          for l in range(Nl):
              if a[l] != 0:
                  if a[l] == 1:
                      assign += "+"
                  elif a[l] == -1:
                      assign += "-"
                  assign += "detuning_knob["+str(l)+"]"
  
      return assign