docs: def angular_momentum_matrix(J, ind="z"):
      ur"""Return the angular momentum operator matrix (divided by hbar) for a\
      given J angular momentum.
  
      INPUT:
  
      -  ``ind`` - A string ("x", "y", "z", "all") indicating which direction \
      to calculate, or to return them all as :math:`(J_x, J_y, J_z)`.
  
      OUTPUT:
  
      - matrix forms of angular momentum operators in the basis \
      :math:`[|J, -J\rangle, \cdot, |J, J\rangle]`.
  
      >>> from sympy import Integer, pprint
      >>> pprint(angular_momentum_matrix(1/Integer(2)))
      ⎡-1/2   0 ⎤
      ⎢         ⎥
      ⎣ 0    1/2⎦
  
      >>> pprint(angular_momentum_matrix(1/Integer(2), "all"))
      ⎛            ⎡     ⅈ⎤             ⎞
      ⎜            ⎢ 0   ─⎥             ⎟
      ⎜⎡ 0   1/2⎤  ⎢     2⎥  ⎡-1/2   0 ⎤⎟
      ⎜⎢        ⎥, ⎢      ⎥, ⎢         ⎥⎟
      ⎜⎣1/2   0 ⎦  ⎢-ⅈ    ⎥  ⎣ 0    1/2⎦⎟
      ⎜            ⎢───  0⎥             ⎟
      ⎝            ⎣ 2    ⎦             ⎠
  
      >>> pprint(angular_momentum_matrix(1, "all"))
      ⎛⎡    √2    ⎤  ⎡         √2⋅ⅈ       ⎤            ⎞
      ⎜⎢0   ──  0 ⎥  ⎢  0      ────    0  ⎥            ⎟
      ⎜⎢    2     ⎥  ⎢          2         ⎥            ⎟
      ⎜⎢          ⎥  ⎢                    ⎥  ⎡-1  0  0⎤⎟
      ⎜⎢√2      √2⎥  ⎢-√2⋅ⅈ           √2⋅ⅈ⎥  ⎢        ⎥⎟
      ⎜⎢──  0   ──⎥, ⎢──────    0     ────⎥, ⎢0   0  0⎥⎟
      ⎜⎢2       2 ⎥  ⎢  2              2  ⎥  ⎢        ⎥⎟
      ⎜⎢          ⎥  ⎢                    ⎥  ⎣0   0  1⎦⎟
      ⎜⎢    √2    ⎥  ⎢        -√2⋅ⅈ       ⎥            ⎟
      ⎜⎢0   ──  0 ⎥  ⎢  0     ──────   0  ⎥            ⎟
      ⎝⎣    2     ⎦  ⎣          2         ⎦            ⎠
  
      """
      MJ = [-J+i for i in range(2*J+1)]
  
      if ind == "x":
          JX = Matrix([[sqrt((J-mj)*(J+mj+1))/2*KroneckerDelta(mi-1, mj)
                       for mj in MJ] for mi in MJ])
          JX += Matrix([[sqrt((J+mj)*(J-mj+1))/2*KroneckerDelta(mi+1, mj)
                        for mj in MJ] for mi in MJ])
          return JX
      elif ind == "y":
          JY = Matrix([[-I*sqrt((J-mj)*(J+mj+1))/2*KroneckerDelta(mi-1, mj)
                       for mj in MJ] for mi in MJ])
          JY += Matrix([[+I*sqrt((J+mj)*(J-mj+1))/2*KroneckerDelta(mi+1, mj)
                        for mj in MJ] for mi in MJ])
          return JY
      elif ind == "z":
          JZ = Matrix([[mi*KroneckerDelta(mi, mj) for mj in MJ] for mi in MJ])
          return JZ
      elif ind == "all":
          JX = angular_momentum_matrix(J, "x")
          JY = angular_momentum_matrix(J, "y")
          JZ = angular_momentum_matrix(J, "z")
          return JX, JY, JZ