Mercurial > repos > public > sbplib
changeset 1245:0a1c64d3c717 feature/dirac_discr
Avoid indentation of whole function
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 20 Nov 2019 20:30:45 +0100 |
parents | 745dc1f1a069 |
children | 25efceb0c392 |
files | diracDiscr.m |
diffstat | 1 files changed, 34 insertions(+), 35 deletions(-) [+] |
line wrap: on
line diff
diff -r 745dc1f1a069 -r 0a1c64d3c717 diracDiscr.m --- a/diracDiscr.m Wed Nov 20 20:22:53 2019 +0100 +++ b/diracDiscr.m Wed Nov 20 20:30:45 2019 +0100 @@ -36,49 +36,48 @@ if x_s < x(1) || x_s > x(end) ret = zeros(size(x)); return - else - tot_order = m_order+s_order; %This is equiv. to the number of equations solved for - S = []; - M = []; - - % Get interior grid spacing - middle = floor(length(x)/2); - h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids. + end + + tot_order = m_order+s_order; %This is equiv. to the number of equations solved for + S = []; + M = []; - index = sourceIndices(x_s, x, tot_order, h); - - 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))); + % Get interior grid spacing + middle = floor(length(x)/2); + h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids. - quadrature = diag(H); - quadrature_weights = quadrature(index)/h; - - h_polynomial = polynomial(2)-polynomial(1); - b = zeros(tot_order,1); + index = sourceIndices(x_s, x, tot_order, h); - for i = 1:m_order - b(i,1) = x_0^(i-1); - 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))); - for i = 1:tot_order - for j = 1:m_order - M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i); - end - end + quadrature = diag(H); + quadrature_weights = quadrature(index)/h; - for i = 1:tot_order - for j = 1:s_order - S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1); - end - end + h_polynomial = polynomial(2)-polynomial(1); + b = zeros(tot_order,1); - A = [M;S]; - - d = A\b; - ret = x*0; - ret(index) = d/h*h_polynomial; + for i = 1:m_order + b(i,1) = x_0^(i-1); end + for i = 1:tot_order + for j = 1:m_order + M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i); + end + end + + for i = 1:tot_order + for j = 1:s_order + S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1); + end + end + + A = [M;S]; + + d = A\b; + ret = x*0; + ret(index) = d/h*h_polynomial; end