Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 754:5264ce57b573 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 15 Jun 2018 14:01:13 -0700
parents 4ee7d15bd8e6
children c70131daaa6e
comparison
equal deleted inserted replaced
753:44c46bd6913a 754:5264ce57b573
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H) 1 function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
2 2
3 m = length(x); 3 m = length(x);
4 4
5 % Return zeros if x0 is outside grid 5 % Return zeros if x0 is outside grid
6 if(x_0in < x(1) || x_0in > x(end) ) 6 if(x_0in < x(1) || x_0in > x(end) )
9 9
10 else 10 else
11 11
12 fnorm = diag(H); 12 fnorm = diag(H);
13 eta = abs(x-x_0in); 13 eta = abs(x-x_0in);
14 h = x(2)-x(1);
15 tot = m_order+s_order; 14 tot = m_order+s_order;
16 S = []; 15 S = [];
17 M = []; 16 M = [];
17
18 % Get interior grid spacing
19 middle = floor(m/2);
20 h = x(middle+1) - x(middle);
21
18 poss = find(tot*h/2 >= eta); 22 poss = find(tot*h/2 >= eta);
19 23
20 % Ensure that poss is not too long 24 % Ensure that poss is not too long
21 if length(poss) == (tot + 2) 25 if length(poss) == (tot + 2)
22 poss = poss(2:end-1); 26 poss = poss(2:end-1);
23 elseif length(poss) == (tot + 1) 27 elseif length(poss) == (tot + 1)
24 poss = poss(1:end-1); 28 poss = poss(1:end-1);