docs: def fast_doppler_terms(v, detuning_knob, k, omega_level, xi, theta,
                         unfolding, axes=["x", "y", "z"], matrix_form=False,
                         file_name=None,
                         return_code=False):
      r"""Return a fast function that returns the Doppler terms.
  
      >>> from sympy import Matrix, symbols
      >>> from scipy.constants import physical_constants
      >>> from fast import PlaneWave
      >>> from fast.bloch import phase_transformation
      >>> Ne = 2
      >>> Nl = 1
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      >>> a0 = physical_constants["Bohr radius"][0]
      >>> rm = [Matrix([[0, 0], [a0, 0]]),
      ...       Matrix([[0, 0], [0, 0]]),
      ...       Matrix([[0, 0], [0, 0]])]
      >>> xi = np.array([[[0, 1], [1, 0]]])
      >>> theta = phase_transformation(Ne, Nl, rm, xi)
  
      >>> omega_level = [1, 2.4e15]
      >>> detuning_knob = [symbols("delta1")]
  
      >>> v = symbols("vx vy vz")
      >>> laser = PlaneWave(0, 0, 0, 0, 1)
      >>> k = [laser.k]
  
      >>> doppler_terms = fast_doppler_terms(v, detuning_knob, k, omega_level,
      ...                                    xi, theta, unfolding,
      ...                                    matrix_form=True)
  
      >>> detuning_knobs = [0]
      >>> A, b = doppler_terms([0, 0, 167], detuning_knobs)
      >>> print A/2/np.pi*1e-9
      [[ 0.          0.          0.        ]
       [ 0.          0.          0.21277821]
       [ 0.         -0.21277821  0.        ]]
  
      """
      # We unpack variables.
      if True:
          Ne = unfolding.Ne
          Nrho = unfolding.Nrho
          Nl = xi.shape[0]
          IJ = unfolding.IJ
          Mu = unfolding.Mu
      # We determine which arguments are constants.
      if 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 establish the arguments of the output function.
      if True:
          code = ""
          code += "def doppler_terms(v, "
          if not matrix_form: code += "rho, "
          if variable_detuning_knob: code += "detuning_knob, "
          if code[-2:] == ", ":
              code = code[:-2]
          code += "):\n"
  
          code += '    r"""A fast calculation of the Doppler terms."""\n'
      # We initialize the output and auxiliaries.
      if True:
          # We introduce the factor that multiplies all terms.
          if unfolding.real:
              code += "    fact = 1.0\n\n"
          else:
              code += "    fact = 1.0j\n\n"
  
          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 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"
      # We put in the code to calculate the Doppler shift
      if True:
          # Delta omega = omega_laser /c k.v
          code += "    # The Doppler shift\n"
          dimension = len(v)
          code += "    c = %s\n" % c_num
          for l in range(Nl):
              lp1 = l+1
              code += "    omega_laser%s = %s" % (lp1, omega_min[l])
              code += "+detuning_knob[%s]\n" % l
              code += "    detuning_knob[%s] = 0\n" % l
  
              axeindices = {"x": 0, "y": 1, "z": 2}
              axeindices = [axeindices[axe] for axe in axes]
              for ii in axeindices[:dimension]:
                  code += "    detuning_knob[%s] += -%s" % (l, k[l][ii])
                  if dimension == 1:
                      code += "*v/c*omega_laser%s\n" % lp1
                  else:
                      code += "*v[%s]/c*omega_laser%s\n" % (ii, lp1)
  
      code += "\n"
  
      code_det = detunings_code(Neu, Nl, pairs, omega_levelu, iu0, ju0)
      code += code_det
      code += "\n"
      #####################################
      # There is a term
      # I * Theta_ij * rho_ij = I * (omega_level_j - omega_level_i
      #                              theta_j - theta_i)
      # for all i != j.
      # This term can be re expressed as
      # re(Theta_ij*rho_ij) = - Theta_ij * im(rho_ij)
      # im(Theta_ij*rho_ij) = + Theta_ij * re(rho_ij)
      _omega_level, omega, gamma = define_frequencies(Ne)
      _omega_levelu, omega, gamma = define_frequencies(Neu)
      E0, omega_laser = define_laser_variables(Nl)
      # We build all combinations.
      combs = detunings_combinations(pairs)
      # We add all terms.
      for mu in range(Nrho):
          s, i, j = IJ(mu)
          if i != j:
              _Thetaij = _omega_levelu[u(j)] - _omega_levelu[u(i)]
              _Thetaij += theta[j] - theta[i]
              aux = (_Thetaij, combs, omega_laser,
                     _omega_levelu, omega_levelu, iu0, ju0)
              assign = detunings_rewrite(*aux)
              if assign != "":
                  if s == 0:
                      nu = mu
                  elif s == 1:
                      assign = "-(%s)" % assign
                      nu = Mu(-s, i, j)
                  elif s == -1:
                      assign = "+(%s)" % assign
                      nu = Mu(-s, i, j)
  
                  if matrix_form:
                      term_code = "    A[%s, %s] = %s\n" % (mu, nu, assign)
                  else:
                      term_code = "    rhs[%s] = (%s)*rho[%s]\n" % (mu,
                                                                    assign, nu)
              else:
                  term_code = ""
              code += term_code
      #####################################
      # 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()
  
          doppler_terms = code
          if not return_code:
              exec doppler_terms
      return doppler_terms