--- a/diracDiscr.m	Fri Jun 15 13:02:29 2018 -0700
+++ b/diracDiscr.m	Fri Jun 15 14:01:13 2018 -0700
@@ -1,4 +1,4 @@
-function ret = diracDiscr(x_0in , x , m_order, s_order, H)
+function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
 
 m = length(x);
 
@@ -11,12 +11,16 @@
 
     fnorm = diag(H);
     eta = abs(x-x_0in);
-    h = x(2)-x(1);
     tot = m_order+s_order;
     S = [];
     M = [];
+    
+    % Get interior grid spacing
+    middle = floor(m/2);
+    h = x(middle+1) - x(middle);
+
     poss = find(tot*h/2 >= eta);
-    
+
     % Ensure that poss is not too long
     if length(poss) == (tot + 2)
         poss = poss(2:end-1);

--- a/diracDiscrTest.m	Fri Jun 15 13:02:29 2018 -0700
+++ b/diracDiscrTest.m	Fri Jun 15 14:01:13 2018 -0700
@@ -220,11 +220,122 @@
     end
 end
 
+function testHalfGP(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+
+        % Test halfway between all grid points
+        x0s = 1/2*( x(2:end)+x(1:end-1) );
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testAllGPStaggered(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+        % Test all grid points
+        x0s = x;
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testHalfGPStaggered(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+        % Test halfway between all grid points
+        x0s = 1/2*( x(2:end)+x(1:end-1) );
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function testRandomStaggered(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
+
+        % Test random points within grid boundaries
+        x0s = xl + (xr-xl)*rand(1,300);
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+
+% ============== Setup functions =======================
 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
 
     % Grid
     xl = -3;
-    xr = 2;
+    xr = 900;
+    L = xr-xl;
     m = 101;
     h = (xr-xl)/(m-1);
     g = grid.equidistant(m, {xl, xr});
@@ -237,7 +348,30 @@
     % Moment conditions
     fs = cell(mom_cond,1);
     for p = 0:mom_cond-1
-        fs{p+1} = @(x) x.^p;
+        fs{p+1} = @(x) (x/L).^p;
+    end
+
+end
+
+function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond)
+
+    % Grid
+    xl = -3;
+    xr = 900;
+    L = xr-xl;
+    m = 101;
+    [~, g_dual] = grid.primalDual1D(m, {xl, xr});
+    x = g_dual.points();
+    h = g_dual.h;
+
+    % Quadrature
+    ops = sbp.D1Staggered(m, {xl, xr}, order);
+    H = ops.H_dual;
+
+    % Moment conditions
+    fs = cell(mom_cond,1);
+    for p = 0:mom_cond-1
+        fs{p+1} = @(x) (x/L).^p;
     end
 
 end