1% Sparse operator for a gradient.
2%
3% [Dx1,Dx2,...,Dxn] = sparse_gradient(mask,pmn)
4%
5% Form the gradient using finite differences in all n-dimensions
6%
7% Input:
8%   mask: binary mask delimiting the domain. 1 is inside and 0 outside.
9%         For oceanographic application, this is the land-sea mask.
10%
11%   pmn: scale factor of the grid.
12%
13% Output:
14%   Dx1,Dx2,...,Dxn: sparce matrix represeting a gradient along
15%     different dimensions
16
17function varargout = sparse_gradient(mask,pmn,iscyclic)
18
19H = sparse_pack(mask);
20
21sz = size(mask);
22n = size(pmn,ndims(pmn));
23Pmn = reshape(pmn,[prod(sz) n]);
24
25if nargin == 2
26  iscyclic = zeros(n,1);
27end
28
29for i=1:n
30  % staggering operator
31  S = sparse_stagger(sz,i,iscyclic(i));
32
33  % mask for staggered variable
34  m = (S * mask(:) == 1);
35
36  d = m .* (S * Pmn(:,i));
37
38  varargout{i} = sparse_pack(m) * sparse_diag(d) * sparse_diff(sz,i,iscyclic(i)) * H';
39%  varargout{i} = sparse_diag(d) * sparse_diff(sz,i) * H';
40%  size(varargout{i})
41end
42
43
44
45% Copyright (C) 2009 Alexander Barth <a.barth@ulg.ac.be>
46%
47% This program is free software; you can redistribute it and/or modify
48% it under the terms of the GNU General Public License as published by
49% the Free Software Foundation; either version 2 of the License, or
50% (at your option) any later version.
51%
52% This program is distributed in the hope that it will be useful,
53% but WITHOUT ANY WARRANTY; without even the implied warranty of
54% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
55% GNU General Public License for more details.
56%
57% You should have received a copy of the GNU General Public License
58% along with this program; If not, see <http://www.gnu.org/licenses/>.
59
60