divergence function

[Javier Enciso]

Hi All,

I' not sure whether this email already reach you guys, I wasn't member 
of the mailing list, anyway I just resend it.

I'm glad to contribute with the 'divergence' function. This function
is used to compute divergence of vector field and is part of Matlab's
core functions.

In doubt, please feel free to contact me.

Regards,
Javier

## Copyright (C) 2009   Javier Enciso   <address@hidden>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see 
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function file} {} divergence (@var{X}, @var{Y}, @var{Z}, 
@var{U}, @var{V}, @var{W})
## @deftypefnx {Function file} {} divergence (@var{U}, @var{V}, @var{W})
## @deftypefnx {Function file} {} divergence (@var{X}, @var{Y}, @var{U}, 
@var{V})
## @deftypefnx {Function file} {} divergence (@var{U}, @var{V})
## Compute divergence of vector field. 
## 
## @code{divergence (X, Y, Z, U, V, W)} computes the divergence of a 3-D vector 
field 
## @var{U}, @var{V}, @var{W} at the points @var{X}, @var{Y}, @var{Z}. 
Coordinate points 
## must preserve the given order and define a 3-dimensional grid.
## 
## If @var{X}, @var{Y}, and @var{Z} are omitted, @code{divergence (U, V, W)} 
assumes 
## that the coordinates are given by the expression @code{[X, Y, Z] = meshgrid 
(1:n, 1:m, 1:p)}
## where @code{[m, n, p] = size (U)}.
## 
## @code{divergence (X, Y, U, V)} computes the divergence of a 2-D vector field 
@var{U}, 
## @var{V} at the points @var{X}, @var{Y}. Coordinate points must preserve the 
given order 
## and define a 2-dimensional grid.
## 
## If @var{X} and @var{Y} are omitted, @code{divergence (U, V)} assumes that 
the 
## coordinate points @var{X} and @var{Y} are given by the expression @code{[X, 
Y] = 
## meshgrid (1:n, 1:m)} where @code{[m, n] = size (U)}.
## @end deftypefn

function [div] = divergence (varargin)
        switch (nargin)
                case 2
                        u = varargin {1};
                        v = varargin {2};
                        check_dim (u, v);
                        [du, dump] = gradient (u);
                        [dump, dv] = gradient (v);
                        div = du + dv;                  
                case 3
                        u = varargin {1};
                        v = varargin {2};
                        w = varargin {3};
                        check_dim (u, v, w);
                        [du, dump, dump] = gradient (u);
                        [dump, dv, dump] = gradient (v);
                        [dump, dump, dw] = gradient (w);
                        div = du + dv + dw;                     
                case 4
                        x = varargin {1};
                        y = varargin {2};
                        u = varargin {3};
                        v = varargin {4};
                        check_dim (x, y, u, v);
                        [du, dump] = gradient (u);
                        [dump, dv] = gradient (v);
                        [dx, dump] = gradient (x);
                        [dump, dy] = gradient (y);
                        div = du./dx + dv./dy;
                case 6
                        x = varargin {1};
                        y = varargin {2};
                        z = varargin {3};
                        u = varargin {4};
                        v = varargin {5};
                        w = varargin {6};
                        check_dim (x, y, z, u, v, w);
                        [du, dump, dump] = gradient (u);
                        [dump, dv, dump] = gradient (v);
                        [dump, dump, dw] = gradient (w);
                        [dx, dump, dump] = gradient (x);
                        [dump, dy, dump] = gradient (y);
                        [dump, dump, dz] = gradient (z);
                        div = du./dx + dv./dy + dw./dz;
                otherwise
                        print_usage ();
        endswitch
endfunction

function [] = check_dim (varargin)
        for m = 1:nargin-1
                if (!size_equal (varargin{m}, varargin{m+1}))
                        error ("Size of matrices must be equal");
                endif
        endfor
        
        switch (nargin)
                case {2, 4}
                        for m = 1:nargin
                                if (ndims (varargin{m} ) != 2)
                                        error ("Matrices must have 2 
dimensions");
                                endif
                        endfor
                case {3, 6}                     
                        for m = 1:nargin                        
                                if (ndims (varargin{m}) != 3)
                                        error ("Matrices must have 3 
dimensions");
                                endif                           
                        endfor          
        endswitch
endfunction

%!test
%! a = [1, 3; 1, 2]; 
%! b = [3, 3; 1, 4];
%! expect = [0, 3; -1, 2];
%! get = divergence (a, b);
%! if (any(size (expect) != size (get)))
%!    error ("wrong size: expected %d,%d but got %d,%d",
%!           size (expect), size (get));
%! elseif (any (any (expect!=get)))
%!    error ("didn't get what was expected.");
%! endif

%!test
%! a(:,:,1) = [1, 3; 1, 2];
%! a(:,:,2) = [2, 1; 2, 2];
%! b(:,:,1) = [5, 5; 0, 4];
%! b(:,:,2) = [6, 2; 4, 2];
%! c(:,:,1) = [0, 1; 2, 3];
%! c(:,:,2) = [4, 5; 5, 4];
%! expect(:,:,1) = [1, 5; -1, 1];
%! expect(:,:,2) = [1, 3; 1, 1];
%! get = divergence (a, b, c);
%! if (any(size (expect) != size (get)))
%!    error ("wrong size: expected %d,%d but got %d,%d",
%!           size (expect), size (get));
%! elseif (any (any (expect!=get)))
%!    error ("didn't get what was expected.");
%! endif

%!test
%! a = rand (3, 3, 3);
%! b = rand (3, 3, 3);
%! c = rand (3, 3, 4);
%!error divergence (a, b, c);

%!test
%! a = rand (2, 3);
%! b = rand (2, 3);
%! c = rand (2, 3);
%!error divergence (a, b, c);