--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscr.m	Tue Nov 14 14:28:35 2017 -0800
@@ -0,0 +1,64 @@
+function ret = diracDiscr(x_0in , x , m_order, s_order, H)
+
+fnorm = diag(H);
+eta = abs(x-x_0in);
+h = x(2)-x(1);
+tot = m_order+s_order;
+S = [];
+M = [];
+poss = find(tot*h/2>eta);
+
+% Use first tot grid points
+if length(poss)<tot && eta(end)>eta(1)
+    index=1:tot;
+    pol=(x(1:tot)-x(1))/(x(tot)-x(1));
+    x_0=(x_0in-x(1))/(x(tot)-x(1));
+    norm=fnorm(1:tot)/h;
+
+% Use last tot grid points
+elseif length(poss)<tot && eta(end)<eta(1)
+    index = length(x)-tot:tot;
+    pol = (x(end-tot:end)-x(end-tot))/(x(end)-x(end-tot));
+    norm = fnorm(end-tot:end)/h;
+    x_0 = (x_0in-x(end-tot))/(x(end)-x(end-tot));
+
+% Interior
+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;
+    index = poss;
+end
+
+h_pol = pol(2)-pol(1);
+b = zeros(m_order+s_order,1);
+
+for i = 1:m_order
+    b(i,1) = x_0^(i-1);
+end
+
+for i = 1:(m_order+s_order)
+    for j = 1:m_order
+        M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
+    end
+end
+
+for i = 1:(m_order+s_order)
+    for j = 1:s_order
+        S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
+    end
+end
+
+
+A = [M;S];
+
+d = A\b;
+ret = x*0;
+ret(index) = d/h*h_pol;
+end
+
+
+
+
+
+