Mercurial > repos > public > sbplib
annotate diracDiscr.m @ 649:1bdbe026abbc feature/d1_staggered
Make diracDiscr return zeros if x0 is outside grid.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Tue, 14 Nov 2017 15:40:06 -0800 |
parents | 9e5dd0d3cf60 |
children | 4ee7d15bd8e6 |
rev | line source |
---|---|
645
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H) |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
2 |
649
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
3 % 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
|
4 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
|
5 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
6 ret = zeros(size(x)); |
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 else |
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 fnorm = diag(H); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
11 eta = abs(x-x_0in); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
12 h = x(2)-x(1); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
13 tot = m_order+s_order; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
14 S = []; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
15 M = []; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
16 poss = find(tot*h/2 >= eta); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
17 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
18 % 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
|
19 if length(poss) == (tot + 2) |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
20 poss = poss(2:end-1); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
21 elseif length(poss) == (tot + 1) |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
22 poss = poss(1:end-1); |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
23 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
24 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
25 % Use first tot grid points |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
26 if length(poss)<tot && eta(end)>eta(1) |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
27 index=1:tot; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
28 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
|
29 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
|
30 norm=fnorm(1:tot)/h; |
646
0990765e3e4d
Fix bug in diracDiscr
Martin Almquist <malmquist@stanford.edu>
parents:
645
diff
changeset
|
31 |
649
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
32 % Use last tot grid points |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
33 elseif length(poss)<tot && eta(end)<eta(1) |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
34 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
|
35 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
|
36 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
|
37 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
|
38 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
39 % Interior |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
40 else |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
41 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
|
42 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
|
43 norm = fnorm(poss)/h; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
44 index = poss; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
45 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
46 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
47 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
|
48 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
|
49 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
50 for i = 1:m_order |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
51 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
|
52 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
53 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
54 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
|
55 for j = 1:m_order |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
56 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
|
57 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
58 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
59 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
60 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
|
61 for j = 1:s_order |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
62 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
|
63 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
64 end |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
65 |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
66 A = [M;S]; |
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 d = A\b; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
69 ret = x*0; |
1bdbe026abbc
Make diracDiscr return zeros if x0 is outside grid.
Martin Almquist <malmquist@stanford.edu>
parents:
648
diff
changeset
|
70 ret(index) = d/h*h_pol; |
646
0990765e3e4d
Fix bug in diracDiscr
Martin Almquist <malmquist@stanford.edu>
parents:
645
diff
changeset
|
71 end |
645
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
72 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
73 end |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
74 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
75 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
76 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
77 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
78 |
a5307adc6177
Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff
changeset
|
79 |