Mercurial > repos > public > sbplib

annotate diracDiscr.m @ 645:a5307adc6177 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add diracDiscr.m. Noticed bug if source is on grid point.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 14 Nov 2017 14:28:35 -0800
parents
children 0990765e3e4d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
3 fnorm = diag(H);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
4 eta = abs(x-x_0in);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
5 h = x(2)-x(1);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
6 tot = m_order+s_order;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
7 S = [];
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
8 M = [];
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
9 poss = find(tot*h/2>eta);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
10
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
11 % Use first tot grid points
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
12 if length(poss)<tot && eta(end)>eta(1)
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
13 index=1:tot;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
14 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
15 x_0=(x_0in-x(1))/(x(tot)-x(1));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
16 norm=fnorm(1:tot)/h;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
17
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
18 % Use last tot grid points
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
19 elseif length(poss)<tot && eta(end)<eta(1)
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
20 index = length(x)-tot:tot;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
21 pol = (x(end-tot:end)-x(end-tot))/(x(end)-x(end-tot));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
22 norm = fnorm(end-tot:end)/h;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
23 x_0 = (x_0in-x(end-tot))/(x(end)-x(end-tot));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
24
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
25 % Interior
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
26 else
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
27 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
28 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
29 norm = fnorm(poss)/h;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
30 index = poss;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
31 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
32
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
33 h_pol = pol(2)-pol(1);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
34 b = zeros(m_order+s_order,1);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
35
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
36 for i = 1:m_order
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
37 b(i,1) = x_0^(i-1);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
38 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
39
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
40 for i = 1:(m_order+s_order)
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
41 for j = 1:m_order
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
42 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
43 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
44 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
45
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
46 for i = 1:(m_order+s_order)
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
47 for j = 1:s_order
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
48 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
49 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
50 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
51
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
52
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
53 A = [M;S];
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
54
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
55 d = A\b;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
56 ret = x*0;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
57 ret(index) = d/h*h_pol;
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
58 end
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
59
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
60
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
61
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
62
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
63
a5307adc6177 Add diracDiscr.m. Noticed bug if source is on grid point.
Martin Almquist <malmquist@stanford.edu>
parents:
diff changeset
64