1 ## Copyright (C) 2009-2012 Kai Habel
3 ## This file is part of Octave.
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.
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.
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/>.
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.
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)$$
35 ## curl F(x,y,z) = | -- Fz - -- Fy, -- Fx - -- Fz, -- Fy - -- Fx |
36 ## \ dy dz dz dx dx dy /
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
47 ## @seealso{divergence, gradient, del2, cross}
50 ## Author: Kai Habel <kai.habel@gmx.de>
52 function varargout = curl (varargin)
56 sz = size (varargin{fidx});
60 sz = size (varargin{fidx});
66 dx = varargin{1}(1,:);
67 dy = varargin{2}(:,1);
70 dx = varargin{1}(1,:,1)(:);
71 dy = varargin{2}(:,1,1)(:);
72 dz = varargin{3}(1,1,:)(:);
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");
87 dFx_dy = gradient (varargin{fidx}.', dy, dx).';
88 dFy_dx = gradient (varargin{fidx + 1}, dx, dy);
89 rot_z = dFy_dx - dFx_dy;
91 if (nargout == 0 || nargout == 1)
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");
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);
120 if (nargout == 0 || nargout == 1)
123 varargout{1} = rot_x;
124 varargout{2} = rot_y;
125 varargout{3} = rot_z;
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));