X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fm%2Fgeneral%2Ftrapz.m;fp=octave_packages%2Fm%2Fgeneral%2Ftrapz.m;h=e528df09a26ed04b8210a95bfd000277fc744703;hp=0000000000000000000000000000000000000000;hb=1c0469ada9531828709108a4882a751d2816994a;hpb=63de9f36673d49121015e3695f2c336ea92bc278

diff --git a/octave_packages/m/general/trapz.m b/octave_packages/m/general/trapz.m
new file mode 100644
index 0000000..e528df0
--- /dev/null
+++ b/octave_packages/m/general/trapz.m
@@ -0,0 +1,138 @@
+## Copyright (C) 2000-2012 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{q} =} trapz (@var{y})
+## @deftypefnx {Function File} {@var{q} =} trapz (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{q} =} trapz (@dots{}, @var{dim})
+##
+## Numerically evaluate the integral of points @var{y} using the trapezoidal
+## method.
+## @w{@code{trapz (@var{y})}} computes the integral of @var{y} along the first
+## non-singleton dimension.  When the argument @var{x} is omitted an
+## equally spaced @var{x} vector with unit spacing (1) is assumed.
+## @code{trapz (@var{x}, @var{y})} evaluates the integral with respect
+## to the spacing in @var{x} and the values in @var{y}.  This is useful if
+## the points in @var{y} have been sampled unevenly.
+## If the optional @var{dim} argument is given, operate along this dimension.
+##
+## If @var{x} is not specified then unit spacing will be used.  To scale
+## the integral to the correct value you must multiply by the actual spacing
+## value (deltaX).  As an example, the integral of @math{x^3} over the range
+## [0, 1] is @math{x^4/4} or 0.25.  The following code uses @code{trapz} to
+## calculate the integral in three different ways.
+##
+## @example
+## @group
+## x = 0:0.1:1;
+## y = x.^3;
+## q = trapz (y)
+##   @result{} q = 2.525   # No scaling
+## q * 0.1
+##   @result{} q = 0.2525  # Approximation to integral by scaling
+## trapz (x, y)
+##   @result{} q = 0.2525  # Same result by specifying @var{x}
+## @end group
+## @end example
+##
+## @seealso{cumtrapz}
+## @end deftypefn
+
+## Author:      Kai Habel <kai.habel@gmx.de>
+##
+## also: June 2000 - Paul Kienzle (fixes,suggestions)
+## 2006-05-12 David Bateman - Modified for NDArrays
+
+function z = trapz (x, y, dim)
+
+  if (nargin < 1) || (nargin > 3)
+    print_usage ();
+  endif
+
+  have_xy = have_dim = false;
+
+  if (nargin == 3)
+    have_xy = true;
+    have_dim = true;
+  elseif (nargin == 2)
+    if (! size_equal (x, y) && isscalar (y))
+      dim = y;
+      have_dim = true;
+    else
+      have_xy = true;
+    endif
+  endif
+
+  if (have_xy)
+    nd = ndims (y);
+    sz = size (y);
+  else
+    nd = ndims (x);
+    sz = size (x);
+  endif
+
+  if (! have_dim)
+    ## Find the first non-singleton dimension.
+    (dim = find (sz > 1, 1)) || (dim = 1);
+  else
+    if (!(isscalar (dim) && dim == fix (dim))
+        || !(1 <= dim && dim <= nd))
+      error ("trapz: DIM must be an integer and a valid dimension");
+    endif
+  endif
+
+  n = sz(dim);
+  idx1 = idx2 = repmat ({':'}, [nd, 1]);
+  idx1{dim} = 2 : n;
+  idx2{dim} = 1 : (n - 1);
+
+  if (! have_xy)
+    z = 0.5 * sum (x(idx1{:}) + x(idx2{:}), dim);
+  else
+    if (isvector (x) && !isvector (y))
+      if (length (x) != sz(dim))
+        error ("trapz: length of X and length of Y along DIM must match");
+      endif
+      ## Reshape vector to point along dimension DIM
+      shape = ones (nd, 1);
+      shape(dim) = sz(dim);
+      x = reshape (x, shape);
+      z = 0.5 * sum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim);
+    else
+      if (! size_equal (x, y))
+        error ("trapz: X and Y must have same shape");
+      endif
+      z = 0.5 * sum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim);
+    endif
+  endif
+endfunction
+
+
+%!assert (trapz(1:5), 12)
+%!assert (trapz(0:0.5:2,1:5), 6)
+%!assert (trapz([1:5;1:5].',1), [12,12])
+%!assert (trapz([1:5;1:5],2), [12;12])
+%!assert (trapz(repmat(reshape(1:5,1,1,5),2,2), 3), [12 12; 12 12])
+%!assert (trapz([0:0.5:2;0:0.5:2].',[1:5;1:5].',1), [6, 6])
+%!assert (trapz([0:0.5:2;0:0.5:2],[1:5;1:5],2), [6; 6])
+%!assert (trapz(repmat(reshape([0:0.5:2],1,1,5),2,2), ...
+%!              repmat(reshape(1:5,1,1,5),2,2), 3), [6 6; 6 6])
+%!assert (trapz(0:0.5:2,[(1:5)',(1:5)']), [6, 6])
+%!assert (trapz(0:0.5:2,[(1:5);(1:5)],2), [6; 6])
+%!assert (trapz(0:0.5:2,repmat(reshape(1:5,1,1,5),2,2),3), [6 6; 6 6])
+