docs: def fast_inhomo_bloch_equations(Ep, epsilonp, detuning_knob, T, gamma,
                                  omega_level, rm, xi, theta,
                                  unfolding, inhomogeneity, matrix_form=False,
                                  file_name=None, return_code=False):
      r"""Return a fast function that returns the numeric right-hand sides of \
      inhomogeneous Bloch equations.
  
      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.electric_field import PlaneWave
      >>> from fast.symbolic import (define_laser_variables,
      ...                            polarization_vector)
      >>> from fast.atomic_structure import Atom
      >>> from fast.bloch import phase_transformation
  
      >>> 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, 2.4e15]
      >>> gamma21 = 2*np.pi*6e6
      >>> gamma = np.array([[0, -gamma21], [gamma21, 0]])
      >>> theta = phase_transformation(Ne, Nl, rm, xi)
  
      We define symbolic variables to be used as token arguments.
      >>> Ep, omega_laser = define_laser_variables(Nl)
      >>> laser = PlaneWave(0, 0, 0, 0)
      >>> epsilonp = [laser.epsilonp]
      >>> k = [laser.k]
      >>> detuning_knob = [symbols("delta1", real=True)]
  
      A map to unfold the density matrix.
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      We define the Doppler broadening.
  
      >>> shape = [9]
      >>> stds = [[-4, 4]]
      >>> T = 273.15+20
      >>> mass = Atom("Rb", 87).mass
      >>> aux = (shape, stds, T, mass, detuning_knob, k,
      ...        omega_level, xi, theta, unfolding, ["z", "x", "y"],
      ...        True)
  
      >>> doppler_effect = DopplerBroadening(*aux)
      >>> doppler_effect.domain
      [array([-669.86784872, -502.40088654, -334.93392436, -167.46696218,
                0.        ,  167.46696218,  334.93392436,  502.40088654,
              669.86784872])]
  
      We obtain a function to calculate the Bloch equations.
  
      >>> T_symb = symbols("T", positive=True)
      >>> aux = (Ep, epsilonp, detuning_knob, T_symb, gamma,
      ...        omega_level, rm, xi, theta, unfolding, doppler_effect,
      ...        True)
      >>> bloch_equations = fast_inhomo_bloch_equations(*aux)
  
      We calculate an example.
  
      >>> detuning_knobs = [0]
      >>> Eps = electric_field_amplitude_top(0, 1e-3, 1, "SI")
      >>> Eps *= np.exp(1j*np.pi)
      >>> Eps = [Eps]
      >>> A, b = bloch_equations(Eps, detuning_knobs, T)
      >>> print A[:, 2, 1]*1e-6/2/np.pi
      [ 853.49268666  640.12094531  426.74705849  213.37341005    0.
       -213.37317167 -426.74610495 -640.11879984 -853.49125636]
      >>> print b*1e-6
      [[ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0.  0.]]
  
      """
      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
          try:
              T = float(T)
              variable_T = False
          except:
              variable_T = True
      # We obtain code for the homogeneous terms.
      if True:
          if file_name is not None:
              file_name_bloch = file_name+"_bloch"
          else:
              file_name_bloch = file_name
          aux = (Ep, epsilonp, detuning_knob, gamma, omega_level, rm, xi, theta,
                 unfolding, matrix_form, file_name_bloch, True)
          bloch_equations = fast_bloch_equations(*aux)
          code = bloch_equations+"\n\n"
      # We establish the arguments of the output function.
      if True:
          code += "def inhomo_bloch_equations("
          code_args = ""
          if not matrix_form: code_args += "rho, "
          if variable_Ep: code_args += "Ep, "
          if variable_epsilonp: code_args += "epsilonp, "
          if variable_detuning_knob: code_args += "detuning_knob, "
          code += code_args
          if variable_T: code += "T, "
          code += "inhomogeneity=inhomogeneity, "
          code += "bloch_equations=bloch_equations):\n"
          code += '    r"""A fast calculation of inhomogeneous '
          code += 'Bloch equations."""\n'
      # We initialize the output and auxiliaries.
      if True:
          # We introduce the factor that multiplies all terms.
          sha = str(inhomogeneity.shape)[1:-1]+" "
          if matrix_form:
              code += "    A = np.zeros(("+sha+str(Nrho)+", "+str(Nrho)
              if not unfolding.real:
                  code += "), complex)\n\n"
              else:
                  code += "))\n\n"
              if unfolding.normalized:
                  code += "    b = np.zeros(("+sha+str(Nrho)
                  if not unfolding.real:
                      code += "), complex)\n\n"
                  else:
                      code += "))\n\n"
          else:
              code += "    rhs = np.zeros(("+sha+str(Nrho)
              if not unfolding.real:
                  code += "), complex)\n\n"
              else:
                  code += "))\n\n"
      # We calculate the equations for each ensemble.
      if True:
          if variable_T: code += "    inhomogeneity.reset(T)\n"
          if code_args[-2:] == ", ": code_args = code_args[:-2]
          code += "    homogeneous = bloch_equations("+code_args+")\n\n"
          code += "    terms = inhomogeneity.terms\n"
          code += "    shape = inhomogeneity.shape\n"
          code += "    domain = inhomogeneity.domain\n"
  
          shape = inhomogeneity.shape
          dimension = len(shape)
          if dimension == 1:
              code += "    for i in range(shape[0]):\n"
              code += "        result = terms(domain[0][i], detuning_knob)\n"
              if matrix_form:
                  if unfolding.normalized:
                      code += "        A[i] = homogeneous[0]+result[0]\n"
                      code += "        b[i] = homogeneous[1]+result[1]\n"
                  else:
                      code += "        A[i] = homogeneous+result\n"
              else:
                  code += "        rhs[i] = homogeneous+result\n"
      # 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()
  
          inhomo_bloch_equations = code
          if not return_code:
              exec inhomo_bloch_equations
      return inhomo_bloch_equations