octave_packages/m/general/divergence.m

   1 ## Copyright (C) 2009-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{div} =} divergence (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz})
  21 ## @deftypefnx {Function File} {@var{div} =} divergence (@var{fx}, @var{fy}, @var{fz})
  22 ## @deftypefnx {Function File} {@var{div} =} divergence (@var{x}, @var{y}, @var{fx}, @var{fy})
  23 ## @deftypefnx {Function File} {@var{div} =} divergence (@var{fx}, @var{fy})
  24 ## Calculate divergence of a vector field given by the arrays @var{fx},
  25 ## @var{fy}, and @var{fz} or @var{fx}, @var{fy} respectively.
  26 ## @tex
  27 ## $$
  28 ## div F(x,y,z) = \partial_x{F} + \partial_y{F} + \partial_z{F}
  29 ## $$
  30 ## @end tex
  31 ## @ifnottex
  32 ##
  33 ## @example
  34 ## @group
  35 ##                   d               d               d
  36 ## div F(x,y,z)  =   -- F(x,y,z)  +  -- F(x,y,z)  +  -- F(x,y,z)
  37 ##                   dx              dy              dz
  38 ## @end group
  39 ## @end example
  40 ##
  41 ## @end ifnottex
  42 ## The coordinates of the vector field can be given by the arguments @var{x},
  43 ## @var{y}, @var{z} or @var{x}, @var{y} respectively.
  44 ##
  45 ## @seealso{curl, gradient, del2, dot}
  46 ## @end deftypefn
  47
  48 ## Author: Kai Habel <kai.habel@gmx.de>
  49
  50 function retval = divergence (varargin)
  51
  52   fidx = 1;
  53   if (nargin == 2)
  54     sz = size (varargin{fidx});
  55     dx = (1:sz(2))(:);
  56     dy = (1:sz(1))(:);
  57   elseif (nargin == 3)
  58     sz = size (varargin{fidx});
  59     dx = 1:sz(2);
  60     dy = 1:sz(1);
  61     dz = 1:sz(3);
  62   elseif (nargin == 4)
  63     fidx = 3;
  64     dx = varargin{1}(1,:);
  65     dy = varargin{2}(:,1);
  66   elseif (nargin == 6)
  67     fidx = 4;
  68     dx = varargin{1}(1,:,1)(:);
  69     dy = varargin{2}(:,1,1)(:);
  70     dz = varargin{3}(1,1,:)(:);
  71   else
  72     print_usage();
  73   endif
  74
  75   if ((nargin == 4) || (nargin == 2))
  76     if (!size_equal (varargin{fidx},varargin{fidx + 1}))
  77       error ("divergence: size of X and Y must match");
  78     elseif (ndims (varargin{fidx}) != 2)
  79       error ("divergence: expected two-dimensional matrices X and Y");
  80     elseif ((length (dx) != columns (varargin{fidx})) || (length (dy) != rows (varargin{fidx})))
  81       error ("divergence: size of dx and dy must match the respective dimension of X and Y");
  82     endif
  83
  84     retval = gradient(varargin{fidx}, dx, dy);
  85     retval += gradient(varargin{fidx + 1}.', dy, dx).';
  86
  87   elseif ((nargin == 6) || (nargin == 3))
  88     if (!size_equal (varargin{fidx},varargin{fidx + 1},varargin{fidx + 2}))
  89       error ("divergence: size of X, Y, and Z must match");
  90     elseif (ndims (varargin{fidx}) != 3)
  91       error ("divergence: expected two-dimensional matrices X, Y, and Z");
  92     elseif ((length (dx) != size (varargin{fidx}, 2))
  93          || (length (dy) != size (varargin{fidx}, 1))
  94          || (length (dz) != size (varargin{fidx}, 3)))
  95       error ("divergence: size of dx, dy, and dz must match the respective dimesion of X, Y, and Z");
  96     endif
  97
  98     ## x-direction
  99     retval = gradient (varargin{fidx}, dx, dy, dz);
 100     ## y-direction
 101     retval += shiftdim (gradient (shiftdim (varargin{fidx + 1}, 2), dy), 1);
 102     ## z-direction
 103     retval += shiftdim (gradient (shiftdim (varargin{fidx + 2}, 1), dz), 2);
 104   endif
 105
 106 endfunction
 107
 108 %!test
 109 %! [X,Y]=meshgrid(-20:20,-22:22);
 110 %! div = divergence(X-Y,Y);
 111 %! assert(all(div(:)==2));
 112 %! assert(size_equal(X,Y,div));