Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 646:0990765e3e4d feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug in diracDiscr
author Martin Almquist <malmquist@stanford.edu>
date Tue, 14 Nov 2017 15:01:50 -0800
parents a5307adc6177
children 9e5dd0d3cf60
comparison
equal deleted inserted replaced
645:a5307adc6177 646:0990765e3e4d
4 eta = abs(x-x_0in); 4 eta = abs(x-x_0in);
5 h = x(2)-x(1); 5 h = x(2)-x(1);
6 tot = m_order+s_order; 6 tot = m_order+s_order;
7 S = []; 7 S = [];
8 M = []; 8 M = [];
9 poss = find(tot*h/2>eta); 9 poss = find(tot*h/2 >= eta);
10
11 % Ensure that poss is not too long
12 if length(poss) == (tot + 2)
13 poss = poss(2:end-1);
14 elseif length(poss) == (tot + 1)
15 poss = poss(1:end-1);
16 end
10 17
11 % Use first tot grid points 18 % Use first tot grid points
12 if length(poss)<tot && eta(end)>eta(1) 19 if length(poss)<tot && eta(end)>eta(1)
13 index=1:tot; 20 index=1:tot;
14 pol=(x(1:tot)-x(1))/(x(tot)-x(1)); 21 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
47 for j = 1:s_order 54 for j = 1:s_order
48 S(j,i) = (-1)^(i-1)*pol(i)^(j-1); 55 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
49 end 56 end
50 end 57 end
51 58
52
53 A = [M;S]; 59 A = [M;S];
54 60
55 d = A\b; 61 d = A\b;
56 ret = x*0; 62 ret = x*0;
57 ret(index) = d/h*h_pol; 63 ret(index) = d/h*h_pol;