[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
divergence function
From: |
Javier Enciso |
Subject: |
divergence function |
Date: |
Tue, 29 Sep 2009 11:16:22 +0200 |
User-agent: |
Thunderbird 2.0.0.23 (Windows/20090812) |
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);