def calculate_r_matrices(fine_states, reduced_matrix_elements, q=None,
                         numeric=True, convention=1):
    ur"""Calculate the matrix elements of the electric dipole (in the helicity
    basis).

    We calculate all matrix elements for the D2 line in Rb 87.

    >>> from sympy import symbols, pprint
    >>> red = symbols("r", positive=True)
    >>> reduced_matrix_elements = [[0, -red], [red, 0]]
    >>> g = State("Rb", 87, 5, 0, 1/Integer(2))
    >>> e = State("Rb", 87, 5, 1, 3/Integer(2))
    >>> fine_levels = [g, e]
    >>> r = calculate_r_matrices(fine_levels, reduced_matrix_elements,
    ...     numeric=False)
    >>> pprint(r[0][8:,:8])
    ⎡                √3⋅r                                  ⎤
    ⎢ 0       0      ────    0    0       0     0      0   ⎥
    ⎢                 6                                    ⎥
    ⎢                                                      ⎥
    ⎢      -√15⋅r                       √15⋅r              ⎥
    ⎢ 0    ───────     0     0    0     ─────   0      0   ⎥
    ⎢         12                          60               ⎥
    ⎢                                                      ⎥
    ⎢               -√15⋅r                     √5⋅r        ⎥
    ⎢ 0       0     ───────  0    0       0    ────    0   ⎥
    ⎢                  12                       20         ⎥
    ⎢                                                      ⎥
    ⎢                                                √10⋅r ⎥
    ⎢ 0       0        0     0    0       0     0    ───── ⎥
    ⎢                                                  20  ⎥
    ⎢                                                      ⎥
    ⎢√2⋅r                       -√6⋅r                      ⎥
    ⎢────     0        0     0  ──────    0     0      0   ⎥
    ⎢ 4                           12                       ⎥
    ⎢                                                      ⎥
    ⎢         r                          -r                ⎥
    ⎢ 0       ─        0     0    0      ───    0      0   ⎥
    ⎢         4                           4                ⎥
    ⎢                                                      ⎥
    ⎢                √3⋅r                      -r          ⎥
    ⎢ 0       0      ────    0    0       0    ───     0   ⎥
    ⎢                 12                        4          ⎥
    ⎢                                                      ⎥
    ⎢                                                -√6⋅r ⎥
    ⎢ 0       0        0     0    0       0     0    ──────⎥
    ⎢                                                  12  ⎥
    ⎢                                                      ⎥
    ⎢ 0       0        0     0    0       0     0      0   ⎥
    ⎢                                                      ⎥
    ⎢                        r                             ⎥
    ⎢ 0       0        0     ─    0       0     0      0   ⎥
    ⎢                        2                             ⎥
    ⎢                                                      ⎥
    ⎢                            √6⋅r                      ⎥
    ⎢ 0       0        0     0   ────     0     0      0   ⎥
    ⎢                             6                        ⎥
    ⎢                                                      ⎥
    ⎢                                   √10⋅r              ⎥
    ⎢ 0       0        0     0    0     ─────   0      0   ⎥
    ⎢                                     10               ⎥
    ⎢                                                      ⎥
    ⎢                                          √5⋅r        ⎥
    ⎢ 0       0        0     0    0       0    ────    0   ⎥
    ⎢                                           10         ⎥
    ⎢                                                      ⎥
    ⎢                                                √15⋅r ⎥
    ⎢ 0       0        0     0    0       0     0    ───── ⎥
    ⎢                                                  30  ⎥
    ⎢                                                      ⎥
    ⎢ 0       0        0     0    0       0     0      0   ⎥
    ⎢                                                      ⎥
    ⎣ 0       0        0     0    0       0     0      0   ⎦
    >>> pprint(r[1][8:,:8])
    ⎡       -√3⋅r                                                 ⎤
    ⎢  0    ──────     0      0      0        0       0       0   ⎥
    ⎢         6                                                   ⎥
    ⎢                                                             ⎥
    ⎢√15⋅r                         -√5⋅r                          ⎥
    ⎢─────    0        0      0    ──────     0       0       0   ⎥
    ⎢  12                            20                           ⎥
    ⎢                                                             ⎥
    ⎢                                      -√15⋅r                 ⎥
    ⎢  0      0        0      0      0     ───────    0       0   ⎥
    ⎢                                         30                  ⎥
    ⎢                                                             ⎥
    ⎢               -√15⋅r                          -√5⋅r         ⎥
    ⎢  0      0     ───────   0      0        0     ──────    0   ⎥
    ⎢                  12                             20          ⎥
    ⎢                                                             ⎥
    ⎢                        √3⋅r                                 ⎥
    ⎢  0      0        0     ────    0        0       0       0   ⎥
    ⎢                         6                                   ⎥
    ⎢                                                             ⎥
    ⎢  r                            √3⋅r                          ⎥
    ⎢  ─      0        0      0     ────      0       0       0   ⎥
    ⎢  4                             12                           ⎥
    ⎢                                                             ⎥
    ⎢        √3⋅r                                                 ⎥
    ⎢  0     ────      0      0      0        0       0       0   ⎥
    ⎢         6                                                   ⎥
    ⎢                                                             ⎥
    ⎢                  r                            -√3⋅r         ⎥
    ⎢  0      0        ─      0      0        0     ──────    0   ⎥
    ⎢                  4                              12          ⎥
    ⎢                                                             ⎥
    ⎢                                                       -√3⋅r ⎥
    ⎢  0      0        0      0      0        0       0     ──────⎥
    ⎢                                                         6   ⎥
    ⎢                                                             ⎥
    ⎢  0      0        0      0      0        0       0       0   ⎥
    ⎢                                                             ⎥
    ⎢                        √3⋅r                                 ⎥
    ⎢  0      0        0     ────    0        0       0       0   ⎥
    ⎢                         6                                   ⎥
    ⎢                                                             ⎥
    ⎢                              √30⋅r                          ⎥
    ⎢  0      0        0      0    ─────      0       0       0   ⎥
    ⎢                                15                           ⎥
    ⎢                                                             ⎥
    ⎢                                       √15⋅r                 ⎥
    ⎢  0      0        0      0      0      ─────     0       0   ⎥
    ⎢                                         10                  ⎥
    ⎢                                                             ⎥
    ⎢                                               √30⋅r         ⎥
    ⎢  0      0        0      0      0        0     ─────     0   ⎥
    ⎢                                                 15          ⎥
    ⎢                                                             ⎥
    ⎢                                                        √3⋅r ⎥
    ⎢  0      0        0      0      0        0       0      ──── ⎥
    ⎢                                                         6   ⎥
    ⎢                                                             ⎥
    ⎣  0      0        0      0      0        0       0       0   ⎦
    >>> pprint(r[2][8:,:8])
    ⎡√3⋅r                                           ⎤
    ⎢────     0     0      0     0      0     0    0⎥
    ⎢ 6                                             ⎥
    ⎢                                               ⎥
    ⎢                    √10⋅r                      ⎥
    ⎢  0      0     0    ─────   0      0     0    0⎥
    ⎢                      20                       ⎥
    ⎢                                               ⎥
    ⎢√15⋅r                      √5⋅r                ⎥
    ⎢─────    0     0      0    ────    0     0    0⎥
    ⎢  12                        20                 ⎥
    ⎢                                               ⎥
    ⎢       √15⋅r                     √15⋅r         ⎥
    ⎢  0    ─────   0      0     0    ─────   0    0⎥
    ⎢         12                        60          ⎥
    ⎢                                               ⎥
    ⎢  0      0     0      0     0      0     0    0⎥
    ⎢                                               ⎥
    ⎢                    √6⋅r                       ⎥
    ⎢  0      0     0    ────    0      0     0    0⎥
    ⎢                     12                        ⎥
    ⎢                                               ⎥
    ⎢√3⋅r                        r                  ⎥
    ⎢────     0     0      0     ─      0     0    0⎥
    ⎢ 12                         4                  ⎥
    ⎢                                               ⎥
    ⎢         r                         r           ⎥
    ⎢  0      ─     0      0     0      ─     0    0⎥
    ⎢         4                         4           ⎥
    ⎢                                               ⎥
    ⎢              √2⋅r                      √6⋅r   ⎥
    ⎢  0      0    ────    0     0      0    ────  0⎥
    ⎢               4                         12    ⎥
    ⎢                                               ⎥
    ⎢  0      0     0      0     0      0     0    0⎥
    ⎢                                               ⎥
    ⎢  0      0     0      0     0      0     0    0⎥
    ⎢                                               ⎥
    ⎢                    √15⋅r                      ⎥
    ⎢  0      0     0    ─────   0      0     0    0⎥
    ⎢                      30                       ⎥
    ⎢                                               ⎥
    ⎢                           √5⋅r                ⎥
    ⎢  0      0     0      0    ────    0     0    0⎥
    ⎢                            10                 ⎥
    ⎢                                               ⎥
    ⎢                                 √10⋅r         ⎥
    ⎢  0      0     0      0     0    ─────   0    0⎥
    ⎢                                   10          ⎥
    ⎢                                               ⎥
    ⎢                                        √6⋅r   ⎥
    ⎢  0      0     0      0     0      0    ────  0⎥
    ⎢                                         6     ⎥
    ⎢                                               ⎥
    ⎢                                              r⎥
    ⎢  0      0     0      0     0      0     0    ─⎥
    ⎣                                              2⎦

    """
    magnetic_states = make_list_of_states(fine_states, 'magnetic', verbose=0)
    aux = calculate_boundaries(fine_states, magnetic_states)
    index_list_fine, index_list_hyperfine = aux

    Ne = len(magnetic_states)

    r = [[[0 for j in range(Ne)] for i in range(Ne)] for p in range(3)]

    II = fine_states[0].i

    for p in [-1, 0, 1]:
        for i in range(Ne):
            ei = magnetic_states[i]
            ii = fine_index(i, index_list_fine)

            for j in range(Ne):
                ej = magnetic_states[j]
                jj = fine_index(j, index_list_fine)

                reduced_matrix_elementij = reduced_matrix_elements[ii][jj]
                if reduced_matrix_elementij != 0:

                    ji = ei.j; jj = ej.j
                    fi = ei.f; fj = ej.f
                    mi = ei.m; mj = ej.m

                    rpij = matrix_element(ji, fi, mi, jj, fj, mj,
                                          II, reduced_matrix_elementij, p,
                                          numeric=numeric,
                                          convention=convention)

                    if q == 1:
                        r[p+1][i][j] = rpij*delta_lesser(i, j)
                    elif q == -1:
                        r[p+1][i][j] = rpij*delta_greater(i, j)
                    else:
                        r[p+1][i][j] = rpij

    if not numeric:
        r = [Matrix(ri) for ri in r]
    return r