1% Testing divand in 2 dimensions with independent verification.
2
3% grid of background field
4[xi,yi] = ndgrid(linspace(0,1,10));
5
6mask = ones(size(xi));
7pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
8pn = ones(size(xi)) / (yi(1,2)-yi(1,1));
9
10% grid of observations
11[x,y] = ndgrid(linspace(eps,1-eps,10));
12x = x(:);
13y = y(:);
14v = sin( x*6 ) .* cos( y*6);
15
16
17lenx = .15;
18leny = .15;
19
20lambda = 20;
21
22[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,{lenx,leny},lambda,'diagnostics',1,'primal',1);
23
24
25
26iR = inv(full(s.R));
27iB = full(s.iB);
28H = full(s.H);
29sv = s.sv;
30
31iP = iB + H'*iR*H;
32
33P = inv(iP);
34
35xa2 = P* (H'*iR*v(:));
36
37[fi2] = statevector_unpack(sv,xa2);
38fi2(~s.mask) = NaN;
39
40rms_diff = [];
41rms_diff(end+1) = divand_rms(va,fi2);
42rms_diff(end+1) = divand_rms(diag(P),err(:));
43
44
45if any(rms_diff > 1e-6)
46  error('unexpected large difference');
47end
48
49fprintf('(max rms=%g) ',max(rms_diff));
50
51% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
52%
53% This program is free software; you can redistribute it and/or modify it under
54% the terms of the GNU General Public License as published by the Free Software
55% Foundation; either version 2 of the License, or (at your option) any later
56% version.
57%
58% This program is distributed in the hope that it will be useful, but WITHOUT
59% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
60% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
61% details.
62%
63% You should have received a copy of the GNU General Public License along with
64% this program; if not, see <http://www.gnu.org/licenses/>.
65