Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1245:0a1c64d3c717 feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Avoid indentation of whole function
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 20 Nov 2019 20:30:45 +0100
parents ff613067dec6
children 25efceb0c392
comparison
equal deleted inserted replaced
1244:745dc1f1a069 1245:0a1c64d3c717
34 34
35 % Return zeros if x0 is outside grid 35 % Return zeros if x0 is outside grid
36 if x_s < x(1) || x_s > x(end) 36 if x_s < x(1) || x_s > x(end)
37 ret = zeros(size(x)); 37 ret = zeros(size(x));
38 return 38 return
39 else 39 end
40 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for 40
41 S = []; 41 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
42 M = []; 42 S = [];
43 M = [];
43 44
44 % Get interior grid spacing 45 % Get interior grid spacing
45 middle = floor(length(x)/2); 46 middle = floor(length(x)/2);
46 h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids. 47 h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids.
47 48
48 index = sourceIndices(x_s, x, tot_order, h); 49 index = sourceIndices(x_s, x, tot_order, h);
49 50
50 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1))); 51 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
51 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1))); 52 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
52 53
53 quadrature = diag(H); 54 quadrature = diag(H);
54 quadrature_weights = quadrature(index)/h; 55 quadrature_weights = quadrature(index)/h;
55 56
56 h_polynomial = polynomial(2)-polynomial(1); 57 h_polynomial = polynomial(2)-polynomial(1);
57 b = zeros(tot_order,1); 58 b = zeros(tot_order,1);
58 59
59 for i = 1:m_order 60 for i = 1:m_order
60 b(i,1) = x_0^(i-1); 61 b(i,1) = x_0^(i-1);
61 end
62
63 for i = 1:tot_order
64 for j = 1:m_order
65 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
66 end
67 end
68
69 for i = 1:tot_order
70 for j = 1:s_order
71 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
72 end
73 end
74
75 A = [M;S];
76
77 d = A\b;
78 ret = x*0;
79 ret(index) = d/h*h_polynomial;
80 end 62 end
81 63
64 for i = 1:tot_order
65 for j = 1:m_order
66 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
67 end
68 end
69
70 for i = 1:tot_order
71 for j = 1:s_order
72 S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1);
73 end
74 end
75
76 A = [M;S];
77
78 d = A\b;
79 ret = x*0;
80 ret(index) = d/h*h_polynomial;
82 end 81 end
83 82
84 83
85 function I = sourceIndices(x_s, x, tot_order, h) 84 function I = sourceIndices(x_s, x, tot_order, h)
86 % Find the indices that are within range of of the point source location 85 % Find the indices that are within range of of the point source location