--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscrTest.m	Tue Nov 14 15:35:28 2017 -0800
@@ -0,0 +1,218 @@
+function tests = diracDiscrTest()
+	    tests = functiontests(localfunctions);
+end
+
+function testLeftGP(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 left boundary grid points
+        x0s = xl + [0, h, 2*h];
+
+        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 testLeftRandom(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 random points near left boundary
+        x0s = xl + 2*h*rand(1,10);
+
+        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 testRightGP(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 right boundary grid points
+        x0s = xr-[0, h, 2*h];
+
+        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 testRightRandom(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 random points near right boundary
+        x0s = xr - 2*h*rand(1,10);
+
+        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 testInteriorGP(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 interior grid points
+        m_half = round(m/2);
+        x0s = xl + (m_half-1:m_half+1)*h;
+
+        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 testInteriorRandom(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 random points in interior
+        x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
+
+        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
+
+% x0 outside grid should yield 0 integral!
+function testX0OutsideGrid(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 points outisde grid
+        x0s = [xl-1.1*h, xr+1.1*h];
+
+        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 - 0);
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
+function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
+
+    % Grid
+    xl = -3;
+    xr = 2;
+    m = 21;
+    h = (xr-xl)/(m-1);
+    g = grid.equidistant(m, {xl, xr});
+    x = g.points();
+
+    % Quadrature
+    ops = sbp.D2Standard(m, {xl, xr}, order);
+    H = ops.H;
+
+    % Moment conditions
+    fs = cell(mom_cond,1);
+    for p = 0:mom_cond-1
+        fs{p+1} = @(x) x.^p;
+    end
+
+end
+
+