docs: 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).
  
      >>> import numpy as np
      >>> g = State("Rb", 87, 5, 0, 1/Integer(2))
      >>> e = State("Rb", 87, 5, 1, 3/Integer(2))
      >>> magnetic_states = make_list_of_states([g, e], "magnetic")
  
      To return the rates in 10^6 rad /s:
  
      >>> gamma = np.array(calculate_gamma_matrix(magnetic_states, Omega=1e6))
  
      The :math:`5P_{3/2}, 5S_{1/2}` block of this matrix is
  
      >>> print(gamma[8:, :8]/2/np.pi)
      [[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):
  
      >>> Gamma = [sum([gamma[i][j] for j in range(i)])/2/pi
      ...          for i in range(len(magnetic_states))][8:]
      >>> for Gammai in Gamma: print("{:2.3f}".format(Gammai))
      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
  
      >>> from sympy import Matrix, pprint, symbols
      >>> Gamma = symbols("Gamma", positive=True)
      >>> einsteinA = Matrix([[0, -Gamma], [Gamma, 0]])
      >>> gamma = calculate_gamma_matrix(magnetic_states, einsteinA=einsteinA,
      ...                                numeric=False)
  
      >>> pprint(Matrix(gamma)[8:, :8])
      ⎡ Γ    Γ    Γ                        ⎤
      ⎢ ─    ─    ─   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   Γ ⎦
  
  
      >>> Gamma =Matrix([sum([gamma[i][j] for j in range(i)])
      ...         for i in range(len(magnetic_states))][8:])
      >>> pprint(Gamma)
      ⎡Γ⎤
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎢Γ⎥
      ⎢ ⎥
      ⎣Γ⎦
  
  
      """
      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