1% Form the different components of the background error covariance matrix. 2% 3% [iB_,iB] = divand_background_components(s,alpha) 4% 5% Compute the components of the background error covariance matrix iB_ and 6% their sum based on alpha (the a-dimensional coefficients for norm, gradient, 7% laplacian,...). 8 9function [iB_,iB] = divand_background_components(s,alpha) 10 11WE = s.WE; 12D = s.D; 13coeff = s.coeff; 14n = s.n; 15 16iB_ = cell(length(alpha),1); 17 18% constrain of total norm 19 20iB_{1} = (1/coeff) * (WE'*WE); 21 22 23% loop over all derivatives 24 25for j=2:length(alpha) 26 % exponent of laplacian 27 k = floor((j-2)/2); 28 29 if mod(j,2) == 0 30 % constrain of derivative with uneven order (j-1) 31 % (gradient, gradient*laplacian,...) 32 % normalized by surface 33 34 iB_{j} = sparse(size(D,1),size(D,1)); 35 36 for i=1:n 37 Dx = s.WEss{i} * s.Dx{i} * D^k; 38 %Dx = s.WEs{i} * s.Dxs{i} * D^k; 39 iB_{j} = iB_{j} + Dx'*Dx; 40 end 41 else 42 % constrain of derivative with even order (j-1) 43 % (laplacian, biharmonic,...) 44 45 % normalize by surface of each cell 46 % such that inner produces (i.e. WE'*WE) 47 % become integrals 48 % WD: units length^(n/2) 49 50 WD = WE * D^(k+1); 51 iB_{j} = WD'*WD; 52 end 53 54 iB_{j} = iB_{j}/coeff; 55end 56 57 58% sum all terms of iB 59% iB is adimentional 60iB = alpha(1) * iB_{1}; 61for j=2:length(alpha) 62 iB = iB + alpha(j) * iB_{j}; 63end 64 65% LocalWords: iB divand 66 67% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be> 68% 69% This program is free software; you can redistribute it and/or modify it under 70% the terms of the GNU General Public License as published by the Free Software 71% Foundation; either version 2 of the License, or (at your option) any later 72% version. 73% 74% This program is distributed in the hope that it will be useful, but WITHOUT 75% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 76% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 77% details. 78% 79% You should have received a copy of the GNU General Public License along with 80% this program; if not, see <http://www.gnu.org/licenses/>. 81