Mercurial > hg > octave-nkf

comparison scripts/general/curl.m @ 11428:13f128bd6a6f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add curl and divergence functions
author Kai Habel <kai.habel@gmx.de>
date Fri, 31 Dec 2010 14:00:12 +0100
parents
children 1740012184f9
comparison
equal deleted inserted replaced
11427:dc983f92e774 11428:13f128bd6a6f
1 ## Copyright (C) 2009, 2010 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 @var{fz}
26 ## or @var{fx}, @var{fy} respectively.
27 ## @iftex
28 ## @tex
29 ## $$ 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)$$
30 ## @end tex
31 ## @end iftex
32 ## @ifnottex
33 ## @example
34 ## @group
35 ## / d d d d d d \
36 ## curl F(x,y,z) = | -- Fz - -- Fy, -- Fx - -- Fz, -- Fy - -- Fx |
37 ## \ dy dz dz dx dx dy /
38 ## @end group
39 ## @end example
40 ## @end ifnottex
41 ## The coordinates of the vector field can be given by the arguments @var{x}, @var{y}, @var{z}
42 ## or @var{x}, @var{y} respectively.
43 ## @var{v} calculates the scalar component of the angular velocity vector in direction
44 ## of the z-axis 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}
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)
92 av
93 elseif (nargout == 1)
94 varargout{1} = av;
95 elseif (nargout == 2)
96 varargout{1} = rot_z;
97 varargout{2} = av;
98 else
99 error ("curl: number of output arguments must be 1 or 2.")
100 endif
101
102 elseif ((nargin == 6) || (nargin == 3))
103 if (!size_equal (varargin{fidx}, varargin{fidx + 1}, varargin{fidx + 2}))
104 error ("curl: size of x, y, and z must match")
105 elseif (ndims (varargin{fidx}) != 3)
106 error ("curl: expected two-dimensional matrices x, y, and z.");
107 elseif ((length (dx) != size (varargin{fidx}, 2))
108 || (length (dy) != size (varargin{fidx}, 1))
109 || (length (dz) != size (varargin{fidx}, 3)))
110 error ("curl: size of dx, dy, and dz must match the respective dimesion of x, y, and z.");
111 endif
112
113 [~, dFx_dy, dFx_dz] = gradient (varargin{fidx}, dx, dy, dz);
114 [dFy_dx, ~, dFy_dz] = gradient (varargin{fidx + 1}, dx, dy, dz);
115 [dFz_dx, dFz_dy] = gradient (varargin{fidx + 2}, dx, dy, dz);
116 rot_x = dFz_dy - dFy_dz;
117 rot_y = dFx_dz - dFz_dx;
118 rot_z = dFy_dx - dFx_dy;
119 l = sqrt(varargin{fidx}.^2 + varargin{fidx + 1}.^2 + varargin{fidx + 2}.^2);
120 av = (rot_x .* varargin{fidx} +
121 rot_y .* varargin{fidx + 1} +
122 rot_z .* varargin{fidx + 2}) ./ (2 * l);
123
124 if (nargout == 0)
125 varargout{1} = av;
126 elseif (nargout == 1)
127 varargout{1} = av;
128 elseif (nargout == 2)
129 varargout{1} = rot_x;
130 varargout{2} = rot_y;
131 elseif (nargout == 3)
132 varargout{1} = rot_x;
133 varargout{2} = rot_y;
134 varargout{3} = rot_z;
135 elseif (nargout == 4)
136 varargout{1} = rot_x;
137 varargout{2} = rot_y;
138 varargout{3} = rot_z;
139 varargout{4} = av;
140 else
141 error ("curl: number of output arguments must be smaller than 5");
142 endif
143 endif
144
145 endfunction
146
147 %!test
148 %! [X,Y]=meshgrid(-20:20,-22:22);
149 %! av = curl(2*(X-Y),Y);
150 %! assert(all(av(:)==1));
151 %! [cz,av] = curl(2*(X-Y),Y);
152 %! assert(all(cz(:)==2));
153 %! assert(all(av(:)==1));
154 %! [cz,av] = curl(X/2,Y/2,2*(X-Y),Y);
155 %! assert(all(cz(:)==4));
156 %! assert(all(av(:)==2));
157 %! assert(size_equal(X,Y,cz,av));