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