octave_packages/m/general/curl.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{cx}, @var{cy}, @var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz})
  21 ## @deftypefnx {Function File} {[@var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{fx}, @var{fy})
  22 ## @deftypefnx {Function File} {[@dots{}] =} curl (@var{fx}, @var{fy}, @var{fz})
  23 ## @deftypefnx {Function File} {[@dots{}] =} curl (@var{fx}, @var{fy})
  24 ## @deftypefnx {Function File} {@var{v} =} curl (@dots{})
  25 ## Calculate curl of vector field given by the arrays @var{fx}, @var{fy}, and
  26 ## @var{fz} or @var{fx}, @var{fy} respectively.
  27 ## @tex
  28 ## $$ curl F(x,y,z) = \left( {\partial{d} \over \partial{y}} F_z - {\partial{d} \over \partial{z}} F_y, {\partial{d} \over \partial{z}} F_x - {\partial{d} \over \partial{x}} F_z, {\partial{d} \over \partial{x}} F_y - {\partial{d} \over \partial{y}} F_x \right)$$
  29 ## @end tex
  30 ## @ifnottex
  31 ##
  32 ## @example
  33 ## @group
  34 ##                   / d         d       d         d       d         d     \
  35 ## curl F(x,y,z)  =  | -- Fz  -  -- Fy,  -- Fx  -  -- Fz,  -- Fy  -  -- Fx |
  36 ##                   \ dy        dz      dz        dx      dx        dy    /
  37 ## @end group
  38 ## @end example
  39 ##
  40 ## @end ifnottex
  41 ## The coordinates of the vector field can be given by the arguments @var{x},
  42 ## @var{y}, @var{z} or @var{x}, @var{y} respectively.  @var{v} calculates the
  43 ## scalar component of the angular velocity vector in direction of the z-axis
  44 ## for two-dimensional input.  For three-dimensional input the scalar
  45 ## rotation is calculated at each grid point in direction of the vector field
  46 ## at that point.
  47 ## @seealso{divergence, gradient, del2, cross}
  48 ## @end deftypefn
  49
  50 ## Author: Kai Habel <kai.habel@gmx.de>
  51
  52 function varargout = curl (varargin)
  53
  54   fidx = 1;
  55   if (nargin == 2)
  56     sz = size (varargin{fidx});
  57     dx = (1:sz(2))(:);
  58     dy = (1:sz(1))(:);
  59   elseif (nargin == 3)
  60     sz = size (varargin{fidx});
  61     dx = (1:sz(2))(:);
  62     dy = (1:sz(1))(:);
  63     dz = (1:sz(3))(:);
  64   elseif (nargin == 4)
  65     fidx = 3;
  66     dx = varargin{1}(1,:);
  67     dy = varargin{2}(:,1);
  68   elseif (nargin == 6)
  69     fidx = 4;
  70     dx = varargin{1}(1,:,1)(:);
  71     dy = varargin{2}(:,1,1)(:);
  72     dz = varargin{3}(1,1,:)(:);
  73   else
  74     print_usage();
  75   endif
  76
  77   if ((nargin == 4) || (nargin == 2))
  78     if (!size_equal (varargin{fidx}, varargin{fidx + 1}))
  79       error ("curl: size of X and Y must match");
  80     elseif (ndims (varargin{fidx}) != 2)
  81       error ("curl: expected two-dimensional matrices X and Y");
  82     elseif ((length (dx) != columns (varargin{fidx}))
  83          || (length (dy) != rows (varargin{fidx})))
  84       error ("curl: size of dx and dy must match the respective dimension of X and Y");
  85     endif
  86
  87     dFx_dy = gradient (varargin{fidx}.', dy, dx).';
  88     dFy_dx = gradient (varargin{fidx + 1}, dx, dy);
  89     rot_z = dFy_dx - dFx_dy;
  90     av = rot_z / 2;
  91     if (nargout == 0 || nargout == 1)
  92       varargout{1} = av;
  93     else
  94       varargout{1} = rot_z;
  95       varargout{2} = av;
  96     endif
  97
  98   elseif ((nargin == 6) || (nargin == 3))
  99     if (!size_equal (varargin{fidx}, varargin{fidx + 1}, varargin{fidx + 2}))
 100       error ("curl: size of X, Y, and Z must match");
 101     elseif (ndims (varargin{fidx}) != 3)
 102       error ("curl: expected two-dimensional matrices X, Y, and Z");
 103     elseif ((length (dx) != size (varargin{fidx}, 2))
 104          || (length (dy) != size (varargin{fidx}, 1))
 105          || (length (dz) != size (varargin{fidx}, 3)))
 106       error ("curl: size of dx, dy, and dz must match the respective dimesion of X, Y, and Z");
 107     endif
 108
 109     [~, dFx_dy, dFx_dz] = gradient (varargin{fidx}, dx, dy, dz);
 110     [dFy_dx, ~, dFy_dz] = gradient (varargin{fidx + 1}, dx, dy, dz);
 111     [dFz_dx, dFz_dy] = gradient (varargin{fidx + 2}, dx, dy, dz);
 112     rot_x = dFz_dy - dFy_dz;
 113     rot_y = dFx_dz - dFz_dx;
 114     rot_z = dFy_dx - dFx_dy;
 115     l = sqrt(varargin{fidx}.^2 + varargin{fidx + 1}.^2 + varargin{fidx + 2}.^2);
 116     av = (rot_x .* varargin{fidx} +
 117           rot_y .* varargin{fidx + 1} +
 118           rot_z .* varargin{fidx + 2}) ./ (2 * l);
 119
 120     if (nargout == 0 || nargout == 1)
 121       varargout{1} = av;
 122     else
 123       varargout{1} = rot_x;
 124       varargout{2} = rot_y;
 125       varargout{3} = rot_z;
 126       varargout{4} = av;
 127     endif
 128   endif
 129
 130 endfunction
 131
 132 %!test
 133 %! [X,Y]=meshgrid(-20:20,-22:22);
 134 %! av = curl(2*(X-Y),Y);
 135 %! assert(all(av(:)==1));
 136 %! [cz,av] = curl(2*(X-Y),Y);
 137 %! assert(all(cz(:)==2));
 138 %! assert(all(av(:)==1));
 139 %! [cz,av] = curl(X/2,Y/2,2*(X-Y),Y);
 140 %! assert(all(cz(:)==4));
 141 %! assert(all(av(:)==2));
 142 %! assert(size_equal(X,Y,cz,av));