Mercurial > hg > octave-lyh

comparison scripts/general/cumtrapz.m @ 10690:35adf2a71f3f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use common code block to find first non-singleton dimension.
author Rik <octave@nomad.inbox5.com>
date Tue, 08 Jun 2010 22:09:25 -0700
parents 95c3e38098bf
children fd0a3ac60b0e
comparison
equal deleted inserted replaced
10689:6622772a0add 10690:35adf2a71f3f
15 ## You should have received a copy of the GNU General Public License 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {@var{z} =} cumtrapz (@var{y}) 20 ## @deftypefn {Function File} {@var{z} =} cumtrapz (@var{y})
21 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@var{x}, @var{y}) 21 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@var{x}, @var{y})
22 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@dots{}, @var{dim}) 22 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@dots{}, @var{dim})
23 ## 23 ##
24 ## Cumulative numerical integration using trapezoidal method. 24 ## Cumulative numerical integration using trapezoidal method.
25 ## @code{cumtrapz (@var{y})} computes the cumulative integral of the 25 ## @code{cumtrapz (@var{y})} computes the cumulative integral of the
26 ## @var{y} along the first non-singleton dimension. If the argument 26 ## @var{y} along the first non-singleton dimension. If the argument
27 ## @var{x} is omitted a equally spaced vector is assumed. @code{cumtrapz 27 ## @var{x} is omitted an equally spaced vector is assumed. @code{cumtrapz
28 ## (@var{x}, @var{y})} evaluates the cumulative integral with respect 28 ## (@var{x}, @var{y})} evaluates the cumulative integral with respect
29 ## to @var{x}. 29 ## to @var{x}.
30 ## 30 ##
31 ## @seealso{trapz,cumsum} 31 ## @seealso{trapz,cumsum}
32 ## @end deftypefn 32 ## @end deftypefn
48 have_x = false; 48 have_x = false;
49 have_dim = false; 49 have_dim = false;
50 if (nargin == 3) 50 if (nargin == 3)
51 have_x = true; 51 have_x = true;
52 have_dim = true; 52 have_dim = true;
53 endif 53 elseif (nargin == 2)
54 if (nargin == 2)
55 if (! size_equal (x, y) && isscalar (y)) 54 if (! size_equal (x, y) && isscalar (y))
56 dim = y; 55 dim = y;
57 have_dim = true; 56 have_dim = true;
58 else 57 else
59 have_x = true; 58 have_x = true;
60 endif 59 endif
61 endif 60 endif
62 61
63 if (! have_dim) 62 if (! have_dim)
64 ## Find the first singleton dimension. 63 ## Find the first non-singleton dimension.
65 dim = 0; 64 dim = find (sz > 1, 1);
66 while (dim < nd && sz(dim+1) == 1) 65 if (isempty (dim))
67 dim++;
68 endwhile
69 dim++;
70 if (dim > nd)
71 dim = 1; 66 dim = 1;
72 endif 67 endif
73 else 68 else
74 dim = floor (dim); 69 dim = floor (dim);
75 if (dim < 1 || dim > nd) 70 if (! (isscalar (dim) && 1 <= dim && dim <= nd))
76 error ("cumtrapz: invalid dimension along which to sort"); 71 error ("cumtrapz: invalid dimension DIM");
77 endif 72 endif
78 endif 73 endif
79 74
80 n = sz(dim); 75 n = sz(dim);
81 idx1 = cell (); 76 idx1 = cell ();
88 83
89 if (! have_x) 84 if (! have_x)
90 z = 0.5 * cumsum (x(idx1{:}) + x(idx2{:}), dim); 85 z = 0.5 * cumsum (x(idx1{:}) + x(idx2{:}), dim);
91 else 86 else
92 if (! size_equal (x, y)) 87 if (! size_equal (x, y))
93 error ("cumtrapz: x and y must have same shape"); 88 error ("cumtrapz: X and Y must have the same shape");
94 endif 89 endif
95 z = 0.5 * cumsum ((x(idx1{:}) - x(idx2{:})) .* 90 z = 0.5 * cumsum ((x(idx1{:}) - x(idx2{:})) .*
96 (y(idx1{:}) + y(idx2{:})), dim); 91 (y(idx1{:}) + y(idx2{:})), dim);
97 endif 92 endif
98 93