function d = diracDiscr(x_s, x, m_order, s_order, H)
    % n-dimensional delta function
    % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z].
    % x: cell array of grid point column vectors for each dimension.
    % m_order: Number of moment conditions
    % s_order: Number of smoothness conditions
    % H: cell array of 1D norm matrices

    dim = length(x_s);
    d_1D = cell(dim,1);

    % If 1D, non-cell input is accepted
    if dim == 1 && ~iscell(x)
        d = diracDiscr1D(x_s, x, m_order, s_order, H);

    else
        for i = 1:dim
            d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i});
        end

        d = d_1D{dim};
        for i = dim-1: -1: 1
            % Perform outer product, transpose, and then turn into column vector
            d = (d_1D{i}*d')';
            d = d(:);
        end
    end
end