docs: def define_r_components(Ne, xi=None, explicitly_hermitian=False,
                          helicity=False, real=True, p=None):
      r"""Define the components of the position operators.
  
      In general, these are representations of the position operators x, y, z
  
      >>> define_r_components(2)
      [Matrix([
      [     0, x_{12}],
      [x_{21},      0]]), Matrix([
      [     0, y_{12}],
      [y_{21},      0]]), Matrix([
      [     0, z_{12}],
      [z_{21},      0]])]
  
      We can make these operators explicitly hermitian
  
      >>> define_r_components(2, explicitly_hermitian=True)
      [Matrix([
      [     0, x_{21}],
      [x_{21},      0]]), Matrix([
      [     0, y_{21}],
      [y_{21},      0]]), Matrix([
      [     0, z_{21}],
      [z_{21},      0]])]
  
      Make them real
  
      >>> r = define_r_components(2, real=True, explicitly_hermitian=True)
      >>> print [r[p]-r[p].transpose() for p in range(3)]
      [Matrix([
      [0, 0],
      [0, 0]]), Matrix([
      [0, 0],
      [0, 0]]), Matrix([
      [0, 0],
      [0, 0]])]
  
      We can get the components of the operator in the helicity basis
  
      >>> define_r_components(2, helicity=True)
      [Matrix([
      [        0, r_{-1;12}],
      [r_{-1;21},         0]]), Matrix([
      [       0, r_{0;12}],
      [r_{0;21},        0]]), Matrix([
      [        0, r_{+1;12}],
      [r_{+1;21},         0]])]
  
      And combinations thereof. For instance, let us check that the components
      in the helicity basis produce hermitian operators in the cartesian basis.
  
      >>> r_helicity = define_r_components(2, helicity=True,
      ...                                  explicitly_hermitian=True)
  
      [Matrix([
      [        0, -r_{+1;21}],
      [r_{-1;21},          0]]), Matrix([
      [       0, r_{0;21}],
      [r_{0;21},        0]]), Matrix([
      [        0, -r_{-1;21}],
      [r_{+1;21},          0]])]
  
      >>> r_cartesian = helicity_to_cartesian(r_helicity)
      >>> r_cartesian[0]
      Matrix([
      [                                 0, sqrt(2)*(-r_{+1;21} + r_{-1;21})/2],
      [sqrt(2)*(-r_{+1;21} + r_{-1;21})/2,                                  0]])
  
  
      >>> [(r_cartesian[p]-r_cartesian[p].adjoint()).expand() for p in range(3)]
      [Matrix([
      [0, 0],
      [0, 0]]), Matrix([
      [0, 0],
      [0, 0]]), Matrix([
      [0, 0],
      [0, 0]])]
  
      """
      frequency_sign = p
      if Ne > 9: comma = ","
      else: comma = ""
  
      if helicity:
          names = ["r_{-1;", "r_{0;", "r_{+1;"]
      else:
          names = ["x", "y", "z"]
  
      r = []
      if helicity:
          for p in range(3):
              r_comp = []
              for i in range(Ne):
                  r_row = []
                  for j in range(Ne):
                      if i == j:
                          r_row += [0]
                      elif i > j:
                          r_row += [Symbol(names[p]+str(i+1)+comma+str(j+1)+"}",
                                    real=real)]
                      elif explicitly_hermitian:
                          sign = int((-1)**(p-1))
                          r_row += [sign*conjugate(Symbol(names[2-p]+str(j+1) +
                                                          comma+str(i+1)+"}",
                                                          real=real))]
                      else:
                          r_row += [Symbol(names[p]+str(i+1)+comma+str(j+1)+"}",
                                           real=real)]
                  r_comp += [r_row]
              r_comp = Matrix(r_comp)
              r += [r_comp]
  
      else:
          for p in range(3):
              r_comp = []
              for i in range(Ne):
                  r_row = []
                  for j in range(Ne):
                      if i == j:
                          r_row += [0]
                      elif i > j:
                          r_row += [Symbol(names[p]+r"_{"+str(i+1) +
                                           comma+str(j+1)+"}", real=real)]
                      elif explicitly_hermitian:
                          r_row += [conjugate(Symbol(names[p]+r"_{"+str(j+1) +
                                                     comma+str(i+1)+"}",
                                                     real=real))]
                      else:
                          r_row += [Symbol(names[p]+r"_{"+str(i+1) +
                                           comma+str(j+1)+"}", real=real)]
                  r_comp += [r_row]
              r_comp = Matrix(r_comp)
              r += [r_comp]
  
      # We select only the upper diagonal or lower diagonal components according
      # to the sign r^(+) or r^(-) provided.
      if frequency_sign == 1:
          r = [Matrix([[r[p][i, j]*delta_lesser(i, j)
                       for j in range(Ne)] for i in range(Ne)])
               for p in range(3)]
      elif frequency_sign == -1:
          r = [Matrix([[r[p][i, j]*delta_greater(i, j)
                       for j in range(Ne)] for i in range(Ne)])
               for p in range(3)]
  
      if xi is not None:
          Nl = len(xi)
          for p in range(3):
              for i in range(Ne):
                  for j in range(Ne):
                      zero = True
                      for l in range(Nl):
                          if xi[l][i, j] != 0:
                              zero = False
                      if zero:
                          r[p][i, j] = 0
      return r