<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def calculate_gamma_matrix(magnetic_states, Omega=1, einsteinA=None, numeric=True):
 ur"""Calculate the matrix of decay between states. This function calculates the matrix :math:`\gamma_{ij}` of decay rates between states :math:`|i\rangle` and :math:`|j\rangle` (in the units specified by the Omega argument). To return the rates in 10^6 rad /s: The :math:`5P_{3/2}, 5S_{1/2}` block of this matrix is [[2.0217 2.0217 2.0217 0. 0. 0. 0. 0. ] [2.5271 2.5271 0. 0.6065 0.3033 0.1011 0. 0. ] [2.5271 0. 2.5271 0. 0.3033 0.4043 0.3033 0. ] [0. 2.5271 2.5271 0. 0. 0.1011 0.3033 0.6065] [3.0325 0. 0. 2.0217 1.0108 0. 0. 0. ] [1.5163 1.5163 0. 1.0108 0.5054 1.5163 0. 0. ] [0.5054 2.0217 0.5054 0. 1.5163 0. 1.5163 0. ] [0. 1.5163 1.5163 0. 0. 1.5163 0.5054 1.0108] [0. 0. 3.0325 0. 0. 0. 1.0108 2.0217] [0. 0. 0. 6.065 0. 0. 0. 0. ] [0. 0. 0. 2.0217 4.0433 0. 0. 0. ] [0. 0. 0. 0.4043 3.2347 2.426 0. 0. ] [0. 0. 0. 0. 1.213 3.639 1.213 0. ] [0. 0. 0. 0. 0. 2.426 3.2347 0.4043] [0. 0. 0. 0. 0. 0. 4.0433 2.0217] [0. 0. 0. 0. 0. 0. 0. 6.065 ]] Let us test if all D2 lines decay at the expected rate (6.065 MHz):
 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 6.065 Let us do this symbolically ⎡ Γ Γ Γ ⎤ ⎢ ─ ─ ─ 0 0 0 0 0 ⎥ ⎢ 3 3 3 ⎥ ⎢ ⎥ ⎢5⋅Γ 5⋅Γ Γ Γ Γ ⎥ ⎢─── ─── 0 ── ── ── 0 0 ⎥ ⎢ 12 12 10 20 60 ⎥ ⎢ ⎥ ⎢5⋅Γ 5⋅Γ Γ Γ Γ ⎥ ⎢─── 0 ─── 0 ── ── ── 0 ⎥ ⎢ 12 12 20 15 20 ⎥ ⎢ ⎥ ⎢ 5⋅Γ 5⋅Γ Γ Γ Γ ⎥ ⎢ 0 ─── ─── 0 0 ── ── ──⎥ ⎢ 12 12 60 20 10⎥ ⎢ ⎥ ⎢ Γ Γ Γ ⎥ ⎢ ─ 0 0 ─ ─ 0 0 0 ⎥ ⎢ 2 3 6 ⎥ ⎢ ⎥ ⎢ Γ Γ Γ Γ Γ ⎥ ⎢ ─ ─ 0 ─ ── ─ 0 0 ⎥ ⎢ 4 4 6 12 4 ⎥ ⎢ ⎥ ⎢Γ Γ Γ Γ Γ ⎥ ⎢── ─ ── 0 ─ 0 ─ 0 ⎥ ⎢12 3 12 4 4 ⎥ ⎢ ⎥ ⎢ Γ Γ Γ Γ Γ ⎥ ⎢ 0 ─ ─ 0 0 ─ ── ─ ⎥ ⎢ 4 4 4 12 6 ⎥ ⎢ ⎥ ⎢ Γ Γ Γ ⎥ ⎢ 0 0 ─ 0 0 0 ─ ─ ⎥ ⎢ 2 6 3 ⎥ ⎢ ⎥ ⎢ 0 0 0 Γ 0 0 0 0 ⎥ ⎢ ⎥ ⎢ Γ 2⋅Γ ⎥ ⎢ 0 0 0 ─ ─── 0 0 0 ⎥ ⎢ 3 3 ⎥ ⎢ ⎥ ⎢ Γ 8⋅Γ 2⋅Γ ⎥ ⎢ 0 0 0 ── ─── ─── 0 0 ⎥ ⎢ 15 15 5 ⎥ ⎢ ⎥ ⎢ Γ 3⋅Γ Γ ⎥ ⎢ 0 0 0 0 ─ ─── ─ 0 ⎥ ⎢ 5 5 5 ⎥ ⎢ ⎥ ⎢ 2⋅Γ 8⋅Γ Γ ⎥ ⎢ 0 0 0 0 0 ─── ─── ──⎥ ⎢ 5 15 15⎥ ⎢ ⎥ ⎢ 2⋅Γ Γ ⎥ ⎢ 0 0 0 0 0 0 ─── ─ ⎥ ⎢ 3 3 ⎥ ⎢ ⎥ ⎣ 0 0 0 0 0 0 0 Γ ⎦ ⎡Γ⎤ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎢Γ⎥ ⎢ ⎥ ⎣Γ⎦ """

    Ne = len(magnetic_states)

    fine_states = []
    fine_map = {}
    ii = 0
    for ei in magnetic_states:
        fine = State(ei.element, ei.isotope, ei.n, ei.l, ei.j)
        if fine not in fine_states:
            fine_states += [fine]
            ii += 1
        fine_map.update({ei: ii-1})

    II = magnetic_states[0].i

    gamma = [[0.0 for j in range(Ne)] for i in range(Ne)]
    for i in range(Ne):
        for j in range(i):
            ei = magnetic_states[i]
            ej = magnetic_states[j]
            if einsteinA is not None:
                iii = fine_map[ei]
                jjj = fine_map[ej]
                einsteinAij = einsteinA[iii, jjj]
            else:
                einsteinAij = Transition(ei, ej).einsteinA

            if einsteinAij != 0:
                ji = ei.j; jj = ej.j
                fi = ei.f; fj = ej.f
                mi = ei.m; mj = ej.m

                gammaij = (2*ji+1)
                gammaij *= (2*fi+1)
                gammaij *= (2*fj+1)

                if numeric:
                    gammaij *= float(wigner_6j(ji, fi, II, fj, jj, 1)**2)

                    gammaij *= sum([float(wigner_3j(fj, 1, fi, -mj, q, mi)**2)
                                    for q in [-1, 0, 1]])

                    gammaij *= einsteinAij/Omega
                    gammaij = float(gammaij)
                else:
                    gammaij *= wigner_6j(ji, fi, II, fj, jj, 1)**2

                    gammaij *= sum([wigner_3j(fj, 1, fi, -mj, q, mi)**2
                                    for q in [-1, 0, 1]])

                    gammaij *= einsteinAij/Omega

                gamma[i][j] = gammaij
                gamma[j][i] = -gammaij

    return gamma