docs: def wigner_d_small(J, beta):
      u"""Return the small Wigner d matrix for angular momentum J.
  
      We use the general formula from [Edmonds74]_, equation 4.1.15.
  
      Some examples form [Edmonds74]_:
  
      >>> from sympy import Integer, symbols, pi
      >>> half = 1/Integer(2)
      >>> beta = symbols("beta", real=True)
      >>> wigner_d_small(half, beta)
      Matrix([
      [ cos(beta/2), sin(beta/2)],
      [-sin(beta/2), cos(beta/2)]])
  
      >>> from sympy import pprint
      >>> pprint(wigner_d_small(2*half, beta), use_unicode=True)
      ⎡        2⎛β⎞              ⎛β⎞    ⎛β⎞           2⎛β⎞     ⎤
      ⎢     cos ⎜─⎟        √2⋅sin⎜─⎟⋅cos⎜─⎟        sin ⎜─⎟     ⎥
      ⎢         ⎝2⎠              ⎝2⎠    ⎝2⎠            ⎝2⎠     ⎥
      ⎢                                                        ⎥
      ⎢       ⎛β⎞    ⎛β⎞       2⎛β⎞      2⎛β⎞        ⎛β⎞    ⎛β⎞⎥
      ⎢-√2⋅sin⎜─⎟⋅cos⎜─⎟  - sin ⎜─⎟ + cos ⎜─⎟  √2⋅sin⎜─⎟⋅cos⎜─⎟⎥
      ⎢       ⎝2⎠    ⎝2⎠        ⎝2⎠       ⎝2⎠        ⎝2⎠    ⎝2⎠⎥
      ⎢                                                        ⎥
      ⎢        2⎛β⎞               ⎛β⎞    ⎛β⎞          2⎛β⎞     ⎥
      ⎢     sin ⎜─⎟        -√2⋅sin⎜─⎟⋅cos⎜─⎟       cos ⎜─⎟     ⎥
      ⎣         ⎝2⎠               ⎝2⎠    ⎝2⎠           ⎝2⎠     ⎦
  
      From table 4 in [Edmonds74]_
  
      >>> wigner_d_small(half, beta).subs({beta:pi/2})
      Matrix([
      [ sqrt(2)/2, sqrt(2)/2],
      [-sqrt(2)/2, sqrt(2)/2]])
  
      >>> wigner_d_small(2*half, beta).subs({beta:pi/2})
      Matrix([
      [       1/2,  sqrt(2)/2,       1/2],
      [-sqrt(2)/2,          0, sqrt(2)/2],
      [       1/2, -sqrt(2)/2,       1/2]])
  
      >>> wigner_d_small(3*half, beta).subs({beta:pi/2})
      Matrix([
      [ sqrt(2)/4,  sqrt(6)/4,  sqrt(6)/4, sqrt(2)/4],
      [-sqrt(6)/4, -sqrt(2)/4,  sqrt(2)/4, sqrt(6)/4],
      [ sqrt(6)/4, -sqrt(2)/4, -sqrt(2)/4, sqrt(6)/4],
      [-sqrt(2)/4,  sqrt(6)/4, -sqrt(6)/4, sqrt(2)/4]])
  
      >>> wigner_d_small(4*half, beta).subs({beta:pi/2})
      Matrix([
      [      1/4,  1/2, sqrt(6)/4,  1/2,       1/4],
      [     -1/2, -1/2,         0,  1/2,       1/2],
      [sqrt(6)/4,    0,      -1/2,    0, sqrt(6)/4],
      [     -1/2,  1/2,         0, -1/2,       1/2],
      [      1/4, -1/2, sqrt(6)/4, -1/2,       1/4]])
  
      """
      def prod(x):
          p = 1
          for i, xi in enumerate(x): p = p*xi
          return p
  
      M = [J-i for i in range(2*J+1)]
      d = []
      for Mi in M:
          row = []
          for Mj in M:
  
              # We get the maximum and minimum value of sigma.
              sigmamax = max([-Mi-Mj, J-Mj])
              sigmamin = min([0, J-Mi])
  
              dij = sqrt(factorial(J+Mi)*factorial(J-Mi) /
                         factorial(J+Mj)/factorial(J-Mj))
              terms = [[(-1)**(J-Mi-s),
                        binomial(J+Mj, J-Mi-s),
                        binomial(J-Mj, s),
                        cos(beta/2)**(2*s+Mi+Mj),
                        sin(beta/2)**(2*J-2*s-Mj-Mi)]
                       for s in range(sigmamin, sigmamax+1)]
  
              terms = [prod(term) if 0 not in term else 0 for term in terms]
  
              dij = dij*sum(terms)
              row += [dij]
          d += [row]
  
      return Matrix(d)