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
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]])]