Mercurial > repos > public > sbplib

changeset 1129:b29892853daf feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor diracDiscr.m by moving the helper function diracDiscr1D to a separate file.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 May 2019 18:10:06 -0700
parents 3a9262c045d0
children 99fd66ffe714
files diracDiscr.m diracDiscr1D.m
diffstat 2 files changed, 92 insertions(+), 95 deletions(-) [+]
line wrap: on
line diff
--- a/diracDiscr.m	Tue May 21 17:59:30 2019 -0700
+++ b/diracDiscr.m	Tue May 21 18:10:06 2019 -0700
@@ -26,101 +26,6 @@
             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);
-
-% Return zeros if x0 is outside grid
-if(x_0in < x(1) || x_0in > x(end) )
-
-    ret = zeros(size(x));
-
-else
-
-    fnorm = diag(H);
-    eta = abs(x-x_0in);
-    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);
-    elseif length(poss) == (tot + 1)
-        poss = poss(1:end-1);
-    end
-
-    % Use first tot grid points
-    if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
-        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 && x_0in > x(end) - ceil(tot/2)*h;
-        index = length(x)-tot+1:length(x);
-        pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
-        norm = fnorm(end-tot+1:end)/h;
-        x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
-
-    % Interior, compensate for round-off errors.
-    elseif length(poss) < tot
-        if poss(end)<m
-            poss = [poss; poss(end)+1];
-        else
-            poss = [poss(1)-1; poss];
-        end
-        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;
-
-    % 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
-
 end
 
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diracDiscr1D.m	Tue May 21 18:10:06 2019 -0700
@@ -0,0 +1,92 @@
+% Generates discrete 1D delta function
+function ret = diracDiscr1D(x_0in , x , m_order, s_order, H)
+
+    m = length(x);
+
+    % Return zeros if x0 is outside grid
+    if(x_0in < x(1) || x_0in > x(end) )
+
+        ret = zeros(size(x));
+
+    else
+
+        fnorm = diag(H);
+        eta = abs(x-x_0in);
+        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);
+        elseif length(poss) == (tot + 1)
+            poss = poss(1:end-1);
+        end
+
+        % Use first tot grid points
+        if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
+            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 && x_0in > x(end) - ceil(tot/2)*h;
+            index = length(x)-tot+1:length(x);
+            pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
+            norm = fnorm(end-tot+1:end)/h;
+            x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
+
+        % Interior, compensate for round-off errors.
+        elseif length(poss) < tot
+            if poss(end)<m
+                poss = [poss; poss(end)+1];
+            else
+                poss = [poss(1)-1; poss];
+            end
+            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;
+
+        % 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
+
+end
\ No newline at end of file