docs: def density_matrix_rotation(J_values, alpha, beta, gamma):
      r"""Return a block-wise diagonal Wigner D matrix for that rotates
      a density matrix of an ensemble of particles in definite total
      angular momentum states given by J_values.
  
      >>> from sympy import Integer, pi
      >>> half = 1/Integer(2)
      >>> J_values = [2*half, 0]
      >>> density_matrix_rotation(J_values, 0, pi/2, 0)
      Matrix([
      [       1/2,  sqrt(2)/2,       1/2, 0],
      [-sqrt(2)/2,          0, sqrt(2)/2, 0],
      [       1/2, -sqrt(2)/2,       1/2, 0],
      [         0,          0,         0, 1]])
  
      """
      size = sum([2*J+1 for J in J_values])
      D = zeros(size, size)
      ind0 = 0
      for J in J_values:
          DJ = wigner_d(J, alpha, beta, gamma)
          sizeJ = 2*J+1
          indf = ind0 + sizeJ
          D[ind0: indf, ind0: indf] = DJ
          ind0 += sizeJ
  
      return D