--- a/diracDiscr.m	Wed Sep 05 14:46:39 2018 -0700
+++ b/diracDiscr.m	Wed Sep 05 18:09:18 2018 -0700
@@ -1,4 +1,37 @@
-function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
+
+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
+
+
+% Helper function for 1D delta functions
+function ret = diracDiscr1D(x_0in , x , m_order, s_order, H)
 
 m = length(x);
 
@@ -14,7 +47,7 @@
     tot = m_order+s_order;
     S = [];
     M = [];
-    
+
     % Get interior grid spacing
     middle = floor(m/2);
     h = x(middle+1) - x(middle);
@@ -53,9 +86,9 @@
         x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
         norm = fnorm(poss)/h;
         index = poss;
-        
+
     % Interior
-    else    
+    else
         pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
         x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
         norm = fnorm(poss)/h;