1% Create the advection constrain. 2% 3% c = divand_constr_advec(s,velocity) 4% 5% Create the advection constrain using the specified velocity. 6% 7% Input: 8% s: structure created by divand_background 9% velocity: cell array of velocity vectors 10% 11% Output: 12% c: structure to be used by divand_addc with the following fields: R (a 13% covariance matrix), H (extraction operator) and yo (specified value for 14% the constrain). 15 16function c = divand_constr_advec(s,velocity) 17 18velocity = cat_cell_array(velocity); 19 20% check for NaNs 21nancount = sum(isnan(velocity(:))); 22if nancount > 0 23 error('%d velocity values are equal to NaN',nancount); 24end 25 26mask = s.mask; 27n = s.n; 28iscyclic = s.iscyclic; 29 30sz = size(mask); 31 32A = sparse(0); 33 34vel = reshape(velocity,numel(mask),n); 35 36 37 38for i=1:n 39 S = sparse_stagger(sz,i,iscyclic(i)); 40 m = (S * mask(:) == 1); 41 42 d = vel(:,i); 43 44 A = A + sparse_diag(d(mask==1)) * sparse_pack(mask) * S' * sparse_pack(m)' * s.Dx{i}; 45end 46 47l = size(A,1); 48 49c.H = A; 50c.yo = sparse(l,1); 51c.R = speye(l); 52 53 54% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be> 55% 56% This program is free software; you can redistribute it and/or modify it under 57% the terms of the GNU General Public License as published by the Free Software 58% Foundation; either version 2 of the License, or (at your option) any later 59% version. 60% 61% This program is distributed in the hope that it will be useful, but WITHOUT 62% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 63% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 64% details. 65% 66% You should have received a copy of the GNU General Public License along with 67% this program; if not, see <http://www.gnu.org/licenses/>. 68