Mercurial > hg > octave-lojdl
changeset 11428:13f128bd6a6f
Add curl and divergence functions
author | Kai Habel <kai.habel@gmx.de> |
---|---|
date | Fri, 31 Dec 2010 14:00:12 +0100 |
parents | dc983f92e774 |
children | 4286b8b999ab |
files | ChangeLog NEWS doc/ChangeLog doc/interpreter/arith.txi scripts/ChangeLog scripts/general/curl.m scripts/general/divergence.m scripts/general/module.mk scripts/miscellaneous/unimplemented.m |
diffstat | 9 files changed, 298 insertions(+), 12 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2010-12-31 Kai Habel <kai.habel@gmx.de> + + * NEWS: Add curl and divergence to new functions list. + 2010-12-20 Ben Abbott <bpabbott@mac.com> * README.MacOS, Makefile.am: Update build insructions for MacOS and
--- a/NEWS +++ b/NEWS @@ -313,15 +313,15 @@ ** The following functions are new in Octave 3.4: - accumdim fileread nfields ppjumps strread - bitpack fminbnd nth_element pqpnonneg textread - bitunpack fskipl onCleanup randi uigetdir - blkmm ifelse pbaspect repelems uigetfile - cbrt ishermitian pie3 reset uiputfile - chop isindex powerset rsf2csf uimenu - daspect luupdate ppder saveas whitebg - erfcx merge ppint saveas - + accumdim erfcx nfields pqpnonneg uigetfile + bitpack fileread nth_element randi uiputfile + bitunpack fminbnd onCleanup repelems uimenu + blkmm fskipl pbaspect reset whitebg + cbrt ifelse pie3 rsf2csf + curl ishermitian powerset saveas + chop isindex ppder strread + daspect luupdate ppint textread + divergence merge ppjumps uigetdir ** Using the image function to view images with external programs such as display, xv, and xloadimage is no longer supported. The
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,7 @@ +2010-12-31 Kai Habel <kai.habel@gmx.de> + + * interpreter/arith.txi: Add curl and divergence function. + 2010-12-22 Tatsuro Matsuoka <tmacchant@yahoo.co.jp> * nonlin.txi: Correct docs for using Jacobian with fsolve.
--- a/doc/interpreter/arith.txi +++ b/doc/interpreter/arith.txi @@ -194,10 +194,14 @@ @DOCSTRING(ceil) +@DOCSTRING(curl) + @DOCSTRING(cross) @DOCSTRING(del2) +@DOCSTRING(divergence) + @DOCSTRING(factor) @DOCSTRING(factorial)
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,9 @@ +2010-12-31 Kai Habel <kai.habel@gmx.de> + + * general/curl.m, general/divergence.m: New functions. + * general/module.mk: Add new file to list. + * miscellaneous/unimplemented.m: Remove curl and divergence from list. + 2010-12-30 Ben Abbott <bpabbott@mac.com> * plot/contour3.m, plot/mesh.m, plot/meshc.m, plot/surfc.m: Grid on @@ -104,7 +110,7 @@ 2010-12-16 Ben Abbott <bpabbott@mac.com> - * miscellaneous/unimplemented.m: Remove isprop() form list of + * miscellaneous/unimplemented.m: Remove isprop() from list of unimplemented functions. * plot/isprop.m: New function.
new file mode 100644 --- /dev/null +++ b/scripts/general/curl.m @@ -0,0 +1,157 @@ +## Copyright (C) 2009, 2010 Kai Habel +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{cx}, @var{cy}, @var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz}) +## @deftypefnx {Function File} {[@var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{fx}, @var{fy}) +## @deftypefnx {Function File} {[@dots{}] =} curl (@var{fx}, @var{fy}, @var{fz}) +## @deftypefnx {Function File} {[@dots{}] =} curl (@var{fx}, @var{fy}) +## @deftypefnx {Function File} {@var{v} =} curl (@dots{}) +## Calculate curl of vector field given by the arrays @var{fx}, @var{fy}, and @var{fz} +## or @var{fx}, @var{fy} respectively. +## @iftex +## @tex +## $$ 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)$$ +## @end tex +## @end iftex +## @ifnottex +## @example +## @group +## / d d d d d d \ +## curl F(x,y,z) = | -- Fz - -- Fy, -- Fx - -- Fz, -- Fy - -- Fx | +## \ dy dz dz dx dx dy / +## @end group +## @end example +## @end ifnottex +## The coordinates of the vector field can be given by the arguments @var{x}, @var{y}, @var{z} +## or @var{x}, @var{y} respectively. +## @var{v} calculates the scalar component of the angular velocity vector in direction +## of the z-axis for two-dimensional input. For three-dimensional input the scalar +## rotation is calculated at each grid point in direction of the vector field +## at that point. +## @seealso{divergence, gradient, del2} +## @end deftypefn + +## Author: Kai Habel <kai.habel@gmx.de> + +function varargout = curl (varargin) + + fidx = 1; + if (nargin == 2) + sz = size (varargin{fidx}); + dx = (1:sz(2))(:); + dy = (1:sz(1))(:); + elseif (nargin == 3) + sz = size (varargin{fidx}); + dx = (1:sz(2))(:); + dy = (1:sz(1))(:); + dz = (1:sz(3))(:); + elseif (nargin == 4) + fidx = 3; + dx = varargin{1}(1,:); + dy = varargin{2}(:,1); + elseif (nargin == 6) + fidx = 4; + dx = varargin{1}(1,:,1)(:); + dy = varargin{2}(:,1,1)(:); + dz = varargin{3}(1,1,:)(:); + else + print_usage(); + endif + + if ((nargin == 4) || (nargin == 2)) + if (!size_equal (varargin{fidx}, varargin{fidx + 1})) + error ("curl: size of x and y must match."); + elseif (ndims (varargin{fidx}) != 2) + error ("curl: expected two-dimensional matrices x and y."); + elseif ((length (dx) != columns (varargin{fidx})) + || (length (dy) != rows (varargin{fidx}))) + error ("curl: size of dx and dy must match the respective dimension of x and y"); + endif + + dFx_dy = gradient (varargin{fidx}.', dy, dx).'; + dFy_dx = gradient (varargin{fidx + 1}, dx, dy); + rot_z = dFy_dx - dFx_dy; + av = rot_z / 2; + if (nargout == 0) + av + elseif (nargout == 1) + varargout{1} = av; + elseif (nargout == 2) + varargout{1} = rot_z; + varargout{2} = av; + else + error ("curl: number of output arguments must be 1 or 2.") + endif + + elseif ((nargin == 6) || (nargin == 3)) + if (!size_equal (varargin{fidx}, varargin{fidx + 1}, varargin{fidx + 2})) + error ("curl: size of x, y, and z must match") + elseif (ndims (varargin{fidx}) != 3) + error ("curl: expected two-dimensional matrices x, y, and z."); + elseif ((length (dx) != size (varargin{fidx}, 2)) + || (length (dy) != size (varargin{fidx}, 1)) + || (length (dz) != size (varargin{fidx}, 3))) + error ("curl: size of dx, dy, and dz must match the respective dimesion of x, y, and z."); + endif + + [~, dFx_dy, dFx_dz] = gradient (varargin{fidx}, dx, dy, dz); + [dFy_dx, ~, dFy_dz] = gradient (varargin{fidx + 1}, dx, dy, dz); + [dFz_dx, dFz_dy] = gradient (varargin{fidx + 2}, dx, dy, dz); + rot_x = dFz_dy - dFy_dz; + rot_y = dFx_dz - dFz_dx; + rot_z = dFy_dx - dFx_dy; + l = sqrt(varargin{fidx}.^2 + varargin{fidx + 1}.^2 + varargin{fidx + 2}.^2); + av = (rot_x .* varargin{fidx} + + rot_y .* varargin{fidx + 1} + + rot_z .* varargin{fidx + 2}) ./ (2 * l); + + if (nargout == 0) + varargout{1} = av; + elseif (nargout == 1) + varargout{1} = av; + elseif (nargout == 2) + varargout{1} = rot_x; + varargout{2} = rot_y; + elseif (nargout == 3) + varargout{1} = rot_x; + varargout{2} = rot_y; + varargout{3} = rot_z; + elseif (nargout == 4) + varargout{1} = rot_x; + varargout{2} = rot_y; + varargout{3} = rot_z; + varargout{4} = av; + else + error ("curl: number of output arguments must be smaller than 5"); + endif + endif + +endfunction + +%!test +%! [X,Y]=meshgrid(-20:20,-22:22); +%! av = curl(2*(X-Y),Y); +%! assert(all(av(:)==1)); +%! [cz,av] = curl(2*(X-Y),Y); +%! assert(all(cz(:)==2)); +%! assert(all(av(:)==1)); +%! [cz,av] = curl(X/2,Y/2,2*(X-Y),Y); +%! assert(all(cz(:)==4)); +%! assert(all(av(:)==2)); +%! assert(size_equal(X,Y,cz,av)); \ No newline at end of file
new file mode 100644 --- /dev/null +++ b/scripts/general/divergence.m @@ -0,0 +1,111 @@ +## Copyright (C) 2009, 2010 Kai Habel +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{div} =} divergence (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz}) +## @deftypefnx {Function File} {@var{div} =} divergence (@var{fx}, @var{fy}, @var{fz}) +## @deftypefnx {Function File} {@var{div} =} divergence (@var{x}, @var{y}, @var{fx}, @var{fy}) +## @deftypefnx {Function File} {@var{div} =} divergence (@var{fx}, @var{fy}) +## Calculate divegence of vector field given by the arrays @var{fx}, @var{fy}, and @var{fz} or @var{fx}, @var{fy} respectively. +## @iftex +## @tex +## $$ +## div F(x,y,z) = \partial_x{F} + \partial_y{F} + \partial_z{F} +## $$ +## @end tex +## @end iftex +## @ifnottex +## @example +## @group +## d d d +## div F(x,y,z) = -- F(x,y,z) + -- F(x,y,z) + -- F(x,y,z) +## dx dy dz +## @end group +## @end example +## @end ifnottex +## The coordinates of the vector field can be given by the arguments @var{x}, @var{y}, @var{z} +## or @var{x}, @var{y} respectively. +## +## @seealso{curl, gradient, del2} +## @end deftypefn + +## Author: Kai Habel <kai.habel@gmx.de> + +function retval = divergence (varargin) + + fidx = 1; + if (nargin == 2) + sz = size (varargin{fidx}); + dx = (1:sz(2))(:); + dy = (1:sz(1))(:); + elseif (nargin == 3) + sz = size (varargin{fidx}); + dx = 1:sz(2); + dy = 1:sz(1); + dz = 1:sz(3); + elseif (nargin == 4) + fidx = 3; + dx = varargin{1}(1,:); + dy = varargin{2}(:,1); + elseif (nargin == 6) + fidx = 4; + dx = varargin{1}(1,:,1)(:); + dy = varargin{2}(:,1,1)(:); + dz = varargin{3}(1,1,:)(:); + else + print_usage(); + endif + + if ((nargin == 4) || (nargin == 2)) + if (!size_equal (varargin{fidx},varargin{fidx + 1})) + error ("divergence: size of x and y must match."); + elseif (ndims (varargin{fidx}) != 2) + error ("divergence: expected two-dimensional matrices x and y."); + elseif ((length (dx) != columns (varargin{fidx})) || (length (dy) != rows (varargin{fidx}))) + error ("divergence: size of dx and dy must match the respective dimension of x and y"); + endif + + retval = gradient(varargin{fidx}, dx, dy); + retval += gradient(varargin{fidx + 1}.', dy, dx).'; + + elseif ((nargin == 6) || (nargin == 3)) + if (!size_equal (varargin{fidx},varargin{fidx + 1},varargin{fidx + 2})) + error ("divergence: size of x, y, and z must match") + elseif (ndims (varargin{fidx}) != 3) + error ("divergence: expected two-dimensional matrices x, y, and z."); + elseif ((length (dx) != size (varargin{fidx}, 2)) + || (length (dy) != size (varargin{fidx}, 1)) + || (length (dz) != size (varargin{fidx}, 3))) + error ("divergence: size of dx, dy, and dz must match the respective dimesion of x, y, and z."); + endif + + ## x-direction + retval = gradient (varargin{fidx}, dx, dy, dz); + ## y-direction + retval += shiftdim (gradient (shiftdim (varargin{fidx + 1}, 2), dy), 1); + ## z-direction + retval += shiftdim (gradient (shiftdim (varargin{fidx + 2}, 1), dz), 2); + endif + +endfunction + +%!test +%! [X,Y]=meshgrid(-20:20,-22:22); +%! div = divergence(X-Y,Y); +%! assert(all(div(:)==2)); +%! assert(size_equal(X,Y,div));
--- a/scripts/general/module.mk +++ b/scripts/general/module.mk @@ -23,10 +23,12 @@ general/common_size.m \ general/cplxpair.m \ general/cumtrapz.m \ + general/curl.m \ general/dblquad.m \ general/deal.m \ general/del2.m \ general/display.m \ + general/divergence.m \ general/flipdim.m \ general/fliplr.m \ general/flipud.m \
--- a/scripts/miscellaneous/unimplemented.m +++ b/scripts/miscellaneous/unimplemented.m @@ -144,7 +144,6 @@ "copyobj", "createClassFromWsdl", "createSoapMessage", - "curl", "customverctrl", "daqread", "datacursormode", @@ -160,7 +159,6 @@ "deval", "dialog", "dither", - "divergence", "docopt", "docsearch", "dragrect",