docs: def coupling_matrix_3j(j1, j2, j3):
      ur"""For angular momenta $j_1, j_2, j_3$ the unitary transformation from the \
      uncoupled basis into the $j = (j_1 \oplus j_2)\oplus j_3$ coupled basis.
  
      >>> from sympy import Integer, pprint
      >>> L = 0
      >>> S = 1/Integer(2)
      >>> II = 3/Integer(2)
      >>> pprint(coupling_matrix_3j(L, S, II))
      ⎡                     √3             ⎤
      ⎢0  -1/2   0     0    ──   0    0   0⎥
      ⎢                     2              ⎥
      ⎢                                    ⎥
      ⎢         -√2              √2        ⎥
      ⎢0   0    ────   0     0   ──   0   0⎥
      ⎢          2               2         ⎥
      ⎢                                    ⎥
      ⎢               -√3                  ⎥
      ⎢0   0     0    ────   0   0   1/2  0⎥
      ⎢                2                   ⎥
      ⎢                                    ⎥
      ⎢1   0     0     0     0   0    0   0⎥
      ⎢                                    ⎥
      ⎢    √3                              ⎥
      ⎢0   ──    0     0    1/2  0    0   0⎥
      ⎢    2                               ⎥
      ⎢                                    ⎥
      ⎢          √2              √2        ⎥
      ⎢0   0     ──    0     0   ──   0   0⎥
      ⎢          2               2         ⎥
      ⎢                                    ⎥
      ⎢                              √3    ⎥
      ⎢0   0     0    1/2    0   0   ──   0⎥
      ⎢                              2     ⎥
      ⎢                                    ⎥
      ⎣0   0     0     0     0   0    0   1⎦
  
      """
      idj3 = eye(2*j3+1)
      Jper = perm_j(j1, j2)
      U_Jj3_list = [coupling_matrix_2j(J, j3) for J in Jper]
  
      size = sum([U_Jj3_list[i].shape[0] for i in range(len(Jper))])
      U_Jj3 = zeros(size, size)
      ind0 = 0
      for i, U_Jj3i in enumerate(U_Jj3_list):
          sizeJ = U_Jj3i.shape[0]
          indf = ind0 + sizeJ
          U_Jj3[ind0: indf, ind0: indf] = U_Jj3_list[i]
          ind0 = indf
      return U_Jj3*TensorProduct(coupling_matrix_2j(j1, j2), idj3)