Mercurial > octave-nkf

changeset 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 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	Thu Dec 30 20:18:15 2010 -0500
+++ b/ChangeLog	Fri Dec 31 14:00:12 2010 +0100
@@ -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	Thu Dec 30 20:18:15 2010 -0500
+++ b/NEWS	Fri Dec 31 14:00:12 2010 +0100
@@ -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	Thu Dec 30 20:18:15 2010 -0500
+++ b/doc/ChangeLog	Fri Dec 31 14:00:12 2010 +0100
@@ -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	Thu Dec 30 20:18:15 2010 -0500
+++ b/doc/interpreter/arith.txi	Fri Dec 31 14:00:12 2010 +0100
@@ -194,10 +194,14 @@
 
 @DOCSTRING(ceil)
 
+@DOCSTRING(curl)
+
 @DOCSTRING(cross)
 
 @DOCSTRING(del2)
 
+@DOCSTRING(divergence)
+
 @DOCSTRING(factor)
 
 @DOCSTRING(factorial)
--- a/scripts/ChangeLog	Thu Dec 30 20:18:15 2010 -0500
+++ b/scripts/ChangeLog	Fri Dec 31 14:00:12 2010 +0100
@@ -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.
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/curl.m	Fri Dec 31 14:00:12 2010 +0100
@@ -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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/divergence.m	Fri Dec 31 14:00:12 2010 +0100
@@ -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	Thu Dec 30 20:18:15 2010 -0500
+++ b/scripts/general/module.mk	Fri Dec 31 14:00:12 2010 +0100
@@ -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	Thu Dec 30 20:18:15 2010 -0500
+++ b/scripts/miscellaneous/unimplemented.m	Fri Dec 31 14:00:12 2010 +0100
@@ -144,7 +144,6 @@
   "copyobj", 
   "createClassFromWsdl", 
   "createSoapMessage", 
-  "curl", 
   "customverctrl", 
   "daqread", 
   "datacursormode", 
@@ -160,7 +159,6 @@
   "deval", 
   "dialog", 
   "dither", 
-  "divergence", 
   "docopt", 
   "docsearch", 
   "dragrect",