--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracPrimDiscr.m	Tue May 21 18:44:01 2019 -0700
@@ -0,0 +1,34 @@
+
+function d = diracPrimDiscr(x_s, x, m_order, s_order, H, derivDir)
+    % n-dimensional delta function, with derivative in direction number derivDir
+    % 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
+    default_arg('derivDir', 1);
+
+    dim = length(x_s);
+    d_1D = cell(dim,1);
+
+    % If 1D, non-cell input is accepted
+    if dim == 1 && ~iscell(x)
+        d = diracPrimDiscr1D(x_s, x, m_order, s_order, H);
+    else
+        for i = 1:dim
+            if i == derivDir
+                d_1D{i} = diracPrimDiscr1D(x_s(i), x{i}, m_order, s_order, H{i});
+            else
+                d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i});
+            end
+        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
\ No newline at end of file