Mercurial > repos > public > sbplib

annotate diracDiscr.m @ 755:14f0058356f2 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/grids
author Martin Almquist <malmquist@stanford.edu>
date Fri, 15 Jun 2018 18:10:26 -0700
parents 5264ce57b573
children c70131daaa6e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
754
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H, h)
645
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
2
651
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
3 m = length(x);
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
4
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
5 % Return zeros if x0 is outside grid
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
6 if(x_0in < x(1) || x_0in > x(end) )
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
7
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
8 ret = zeros(size(x));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
9
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
10 else
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
11
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
12 fnorm = diag(H);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
13 eta = abs(x-x_0in);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
14 tot = m_order+s_order;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
15 S = [];
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
16 M = [];
754
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
17
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
18 % Get interior grid spacing
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
19 middle = floor(m/2);
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
20 h = x(middle+1) - x(middle);
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
21
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
22 poss = find(tot*h/2 >= eta);
754
5264ce57b573 Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
Martin Almquist <malmquist@stanford.edu>
parents: 651
diff changeset
23
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
24 % Ensure that poss is not too long
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
25 if length(poss) == (tot + 2)
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
26 poss = poss(2:end-1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
27 elseif length(poss) == (tot + 1)
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
28 poss = poss(1:end-1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
29 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
30
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
31 % Use first tot grid points
651
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
32 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
33 index=1:tot;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
34 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
35 x_0=(x_0in-x(1))/(x(tot)-x(1));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
36 norm=fnorm(1:tot)/h;
646
0990765e3e4d Fix bug in diracDiscr
Martin Almquist <malmquist@stanford.edu>
parents: 645
diff changeset
37
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
38 % Use last tot grid points
651
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
39 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
40 index = length(x)-tot+1:length(x);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
41 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
42 norm = fnorm(end-tot+1:end)/h;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
43 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
44
651
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
45 % Interior, compensate for round-off errors.
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
46 elseif length(poss) < tot
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
47 if poss(end)<m
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
48 poss = [poss; poss(end)+1];
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
49 else
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
50 poss = [poss(1)-1; poss];
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
51 end
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
52 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
53 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
54 norm = fnorm(poss)/h;
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
55 index = poss;
4ee7d15bd8e6 Fix roundoff bug in diracDiscr and improve test.
Martin Almquist <malmquist@stanford.edu>
parents: 649
diff changeset
56
649
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
57 % Interior
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
58 else
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
59 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
60 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
61 norm = fnorm(poss)/h;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
62 index = poss;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
63 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
64
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
65 h_pol = pol(2)-pol(1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
66 b = zeros(m_order+s_order,1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
67
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
68 for i = 1:m_order
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
69 b(i,1) = x_0^(i-1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
70 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
71
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
72 for i = 1:(m_order+s_order)
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
73 for j = 1:m_order
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
74 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
75 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
76 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
77
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
78 for i = 1:(m_order+s_order)
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
79 for j = 1:s_order
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
80 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
81 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
82 end
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
83
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
84 A = [M;S];
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
85
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
86 d = A\b;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
87 ret = x*0;
1bdbe026abbc Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents: 648
diff changeset
88 ret(index) = d/h*h_pol;
646
0990765e3e4d Fix bug in diracDiscr
Martin Almquist <malmquist@stanford.edu>
parents: 645
diff changeset
89 end
645
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
90
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
91 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
92
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
93
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
94
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
95
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
96
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
97