docs: def exclude_states(omega, gamma, r, Lij, states, excluded_states):
      """Exclude states from matrices.
  
      This function takes the matrices and excludes the states listed in
      excluded_states.
      """
      Ne = len(omega)
      excluded_indices = [i for i in range(Ne) if states[i] in excluded_states]
  
      omega_new = []; gamma_new = []; r_new = [[], [], []]; Lij_new = []
      for i in range(Ne):
          row_om = []; row_ga = []; row_L = []
          for j in range(Ne):
              if j not in excluded_indices:
                  row_om += [omega[i][j]]
                  row_ga += [gamma[i][j]]
                  row_L += [Lij[i][j]]
          if i not in excluded_indices:
              omega_new += [row_om]
              gamma_new += [row_ga]
              Lij_new += [row_L]
  
      for p in range(3):
          for i in range(Ne):
              row_r = []
              for j in range(Ne):
                  if j not in excluded_indices:
                      row_r += [r[p][i][j]]
              if i not in excluded_indices:
                  r_new[p] += [row_r]
  
      states_new = [states[i] for i in range(Ne) if i not in excluded_indices]
  
      return omega_new, gamma_new, r_new, Lij_new, states_new