docs: def fast_steady_state(Ep, epsilonp, detuning_knob, gamma,
                        omega_level, rm, xi, theta,
                        file_name=None, return_code=False):
      r"""Return a fast function that returns a steady state.
  
      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.symbolic import (define_laser_variables,
      ...                            polarization_vector)
  
      >>> 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, 1.0e9]
      >>> 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)
      >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]
      >>> detuning_knob = [symbols("delta1", real=True)]
  
      An map to unfold the density matrix.
      >>> unfolding = Unfolding(Ne, True, True, True)
  
      We obtain a function to calculate Hamiltonian terms.
      >>> aux = (Ep, epsilonps, detuning_knob, gamma,
      ...        omega_level, rm, xi, theta)
      >>> steady_state = fast_steady_state(*aux)
  
      We specify values for the variables
      >>> detuning_knobs = [100e6]
      >>> Eps = electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")
      >>> Eps *= np.exp(1j*np.pi)
      >>> Eps = [Eps]
      >>> print(steady_state(Eps, detuning_knobs))
      [ 0.018   0.1296 -0.0244]
  
      """
      # We unpack variables.
      if True:
          Ne = len(omega_level)
          Nl = xi.shape[0]
          unfolding = Unfolding(Ne, True, True, True)
      # 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
      # We obtain code for the three parts.
      if True:
          args = (Ep, epsilonp, detuning_knob, gamma,
                  omega_level, rm, xi, theta,
                  unfolding, True, None, True)
  
          bloch_equations = fast_bloch_equations(*args)
  
          code = bloch_equations+"\n\n"
  
          if ((not variable_Ep) and
             (not variable_epsilonp) and
             (not variable_detuning_knob)):
              # We can call bloch_equations here!
              code += "bloch_equations = bloch_equations()\n"
      # We establish the arguments of the output function.
      if True:
          code += "def steady_state("
          if variable_Ep: code += "Ep, "
          if variable_epsilonp: code += "epsilonp, "
          if variable_detuning_knob: code += "detuning_knob, "
          code += "bloch_equations=bloch_equations):\n"
          code += '    r"""A fast calculation of the steady state."""\n'
      # We call the Bloch equations.
      if True:
          code += r"""    A, b = bloch_equations"""
          if ((not variable_Ep) and
             (not variable_epsilonp) and
             (not variable_detuning_knob)):
  
              code += "\n"
          else:
              code += "("
              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 += """    rhox = np.linalg.solve(A, b)\n"""
          code += """    return rhox\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()
  
          steady_state = code
          if not return_code:
              exec steady_state
      return steady_state