Re: divergence function

[Javier Enciso]

Soren Hauberg wrote:
> tir, 29 09 2009 kl. 16:11 +0200, skrev Javier Enciso:
> > I already made the corrections and I hope the code is ready to be 
> > shipped. If you (or someone else) have any other concern, please let me 
> > know.
>
> A few minor things:
>
>   * In your cases you have, e.g.
>
>       case 2
>         [u, v] = deal (varargin{1:2});
>
>     since you know that 'varargin' has two elements, you can just write
>
>       case 2
>         [u, v] = deal (varargin{:});
>
>   * When calling 'error', the message should start with the name of the
>     function raising the error, i.e.
>
>       error ("divergence: size of matrices must be equal");
>
>   * Instead of
>
>       function [] = check_dim (varargin)
>
>     you should probably just write
>
>       function check_dim (varargin)
>
> Otherwise I think it looks good :-)
>
> Soren

Hi Soren,

Thanks again for your suggestions. Here it goes the correction and if 
any other improvement is spotted, please let me know.

Btw, is there a formal guideline for coding style? I guess I'm missing 
that document. ;-)

Cheers,
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-dimensional 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-dimensional 
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, v] = deal (varargin{:});
      check_dim (u, v);
      [du, dump] = gradient (u);
      [dump, dv] = gradient (v);
      div = du + dv;
    case 3
      [u, v, w] = deal (varargin{:});
      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, y, u, v] = deal (varargin{:});
      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, y, z, u, v, w] = deal (varargin{:});
      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)
  if (!size_equal (varargin{:}))
    error ("divergence: size of matrices must be equal");
  endif

  switch (nargin)
  case {2, 4}
    if (any (cellfun (@ndims, varargin) != 2))
      error ("divergence: matrices must have 2 dimensions");
    endif
  case {3, 6}
    if (any (cellfun (@ndims, varargin) != 3))
      error ("divergence: matrices must have 3 dimensions");
    endif
  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);