Mercurial > repos > public > sbplib

changeset 1236:3722c2579818 feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt to factor out a function for finding indecies of the source
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 20 Nov 2019 00:10:30 +0100
parents 48c9a83260c8
children 6e4cc4b66de0
files diracDiscr.m
diffstat 1 files changed, 30 insertions(+), 50 deletions(-) [+]
line wrap: on
line diff
--- a/diracDiscr.m	Tue Nov 19 23:51:32 2019 +0100
+++ b/diracDiscr.m	Wed Nov 20 00:10:30 2019 +0100
@@ -41,7 +41,6 @@
         ret = zeros(size(x));
 
     else
-
         fnorm = diag(H);
         tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
         S = [];
@@ -51,65 +50,26 @@
         middle = floor(m/2);
         h = x(middle+1) - x(middle);
 
-        % Find the indices that are within range of of the point source location
-        ind_delta = find(tot_order*h/2 >= abs(x-x_s));
-
-        % Ensure that ind_delta is not too long
-        if length(ind_delta) == (tot_order + 2)
-            ind_delta = ind_delta(2:end-1);
-        elseif length(ind_delta) == (tot_order + 1)
-            ind_delta = ind_delta(1:end-1);
-        end
-
-        % Use first tot_order grid points
-        if length(ind_delta)<tot_order && x_s < x(1) + ceil(tot_order/2)*h;
-            index=1:tot_order;
-            polynomial=(x(1:tot_order)-x(1))/(x(tot_order)-x(1));
-            x_0=(x_s-x(1))/(x(tot_order)-x(1));
-            norm=fnorm(1:tot_order)/h;
+        index = sourceIndecies(x_s, x, tot_order, h)
 
-        % Use last tot_order grid points
-        elseif length(ind_delta)<tot_order && x_s > x(end) - ceil(tot_order/2)*h;
-            index = length(x)-tot_order+1:length(x);
-            polynomial = (x(end-tot_order+1:end)-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
-            norm = fnorm(end-tot_order+1:end)/h;
-            x_0 = (x_s-x(end-tot_order+1))/(x(end)-x(end-tot_order+1));
-
-        % Interior, compensate for round-off errors.
-        elseif length(ind_delta) < tot_order
-            if ind_delta(end)<m
-                ind_delta = [ind_delta; ind_delta(end)+1];
-            else
-                ind_delta = [ind_delta(1)-1; ind_delta];
-            end
-            
-            index = ind_delta;
-            polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
-            x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
-            norm = fnorm(ind_delta)/h;
-
-        % Interior
-        else
-            index = ind_delta;
-            polynomial = (x(ind_delta)-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
-            x_0 = (x_s-x(ind_delta(1)))/(x(ind_delta(end))-x(ind_delta(1)));
-            norm = fnorm(ind_delta)/h;
-        end
+        polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
+        x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
+        norm = fnorm(index)/h;
 
         h_polynomial = polynomial(2)-polynomial(1);
-        b = zeros(m_order+s_order,1);
+        b = zeros(tot_order,1);
 
         for i = 1:m_order
             b(i,1) = x_0^(i-1);
         end
 
-        for i = 1:(m_order+s_order)
+        for i = 1:tot_order
             for j = 1:m_order
                 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i);
             end
         end
 
-        for i = 1:(m_order+s_order)
+        for i = 1:tot_order
             for j = 1:s_order
                 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
             end
@@ -125,7 +85,27 @@
 end
 
 
-
+function I = sourceIndecies(x_s, x, tot_order, h)
+    % Find the indices that are within range of of the point source location
+    I = find(tot_order*h/2 >= abs(x-x_s));
 
-
-
+    if length(I) > tot_order
+        if length(I) == tot_order + 2
+            I = I(2:end-1);
+        elseif length(I) == tot_order + 1
+            I = I(1:end-1);
+        end
+    elseif length(I) < tot_order
+        if x_s < x(1) + ceil(tot_order/2)*h;
+            I = 1:tot_order;
+        elseif x_s > x(end) - ceil(tot_order/2)*h;
+            I = length(x)-tot_order+1:length(x);
+        else
+            if I(end) < length(x)
+                I = [I; I(end)+1];
+            else
+                I = [I(1)-1; I];
+            end
+        end
+    end
+end