--- a/libinterp/corefcn/stream-euler.cc	Thu Dec 12 20:35:09 2019 +0100
+++ b/libinterp/corefcn/stream-euler.cc	Thu Dec 12 14:50:17 2019 -0800
@@ -57,7 +57,7 @@
   double x, y;
 } Vector2;
 
-// The integer- and the fractional value from a point in C-Space.
+// The integer value and the fractional value from a point in C-Space.
 // Equivalent to the cell index the point is located in and the local
 // coordinates of the point in the cell.
 
@@ -503,10 +503,11 @@
 DEFUN (__streameuler2d__, args, ,
        doc: /* -*- texinfo -*-
 @deftypefn {} {} __streameuler2d__ (@var{U}, @var{V}, @var{TX}, @var{TY}, @var{ZETA}, @var{XI}, @var{H}, @var{MAXNVERTS})
-Calculates the streamline in a vector field [@var{U}, @var{V}] starting from a
-seed point at position [@var{ZETA}, @var{XI}].  The integrator used is
-Heun's Scheme.  The step size can be controlled by @var{H}.  The Jacobian
-matrix can be defined for each grid cell by [@var{TX}, @var{TY}].
+Calculates the streamline in a vector field @code{[@var{U}, @var{V}]} starting
+from a seed point at position @code{[@var{ZETA}, @var{XI}]}.  The integrator
+used is Heun's Scheme.  The step size can be controlled by @var{H}.  The
+Jacobian matrix can be defined for each grid cell by
+@code{[@var{TX}, @var{TY}]}.
 
 @seealso{streamline, stream2, stream3, __streameuler3d__}
 @end deftypefn */)
@@ -517,11 +518,11 @@
 DEFUN (__streameuler3d__, args, ,
        doc: /* -*- texinfo -*-
 @deftypefn {} {} __streameuler3d__ (@var{U}, @var{V}, @var{W}, @var{TX}, @var{TY}, @var{TZ}, @var{ZETA}, @var{XI}, @var{RHO}, @var{H}, @var{MAXNVERTS})
-Calculates the streamline in a vector field [@var{U}, @var{V}, @var{W}]
-starting from a seed point at position [@var{ZETA}, @var{XI}, @var{RHO}].  The
-integrator used is Heun's Scheme.  The step size can be controlled by @var{H}.
- The Jacobian matrix can be defined for each grid cell by [@var{TX}, @var{TY},
-@var{TZ}].
+Calculates the streamline in a vector field @code{[@var{U}, @var{V}, @var{W}]}
+starting from a seed point at position
+@code{[@var{ZETA}, @var{XI}, @var{RHO}]}.  The integrator used is Heun's
+Scheme.  The step size can be controlled by @var{H}.  The Jacobian matrix can
+be defined for each grid cell by @code{[@var{TX}, @var{TY}, @var{TZ}]}.
 
 @seealso{streamline, stream2, stream3, __streameuler2d__}
 @end deftypefn */)

--- a/scripts/plot/draw/stream2.m	Thu Dec 12 20:35:09 2019 +0100
+++ b/scripts/plot/draw/stream2.m	Thu Dec 12 14:50:17 2019 -0800
@@ -19,7 +19,7 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {@var{xy} =} stream2 (@var{x}, @var{y}, @var{u}, @var{v}, @var{sx}, @var{sy})
 ## @deftypefnx {} {@var{xy} =} stream2 (@var{u}, @var{v}, @var{sx}, @var{sy})
-## @deftypefnx {} {@var{xy} =} stream2 (@dots{}, "@var{options}")
+## @deftypefnx {} {@var{xy} =} stream2 (@dots{}, @var{options})
 ## Compute 2-D streamline data.
 ##
 ## Calculates streamlines of a vector field given by @code{[@var{u}, @var{v}]}.
@@ -30,9 +30,9 @@
 ## @code{[]} is returned.
 ##
 ## The input parameter @var{options} is a 2-D vector of the form
-## @code{[@var{stepsize}, @var{maxnumbervertices}]}.  The first parameter
-## specifies the step size used for trajectory integration (default 0.1).  It
-## is allowed to set a negative value to control the direction of integration.
+## @code{[@var{stepsize}, @var{max_vertices}]}.  The first parameter
+## specifies the step size used for trajectory integration (default 0.1).  A
+## negative value is allowed which will reverse the direction of integration.
 ## The second parameter specifies the maximum number of segments used to
 ## create a streamline (default 10,000).
 ##
@@ -51,7 +51,6 @@
 ## @end example
 ##
 ## @seealso{streamline, stream3}
-##
 ## @end deftypefn
 
 ## References:
@@ -75,52 +74,53 @@
 function xy = stream2 (varargin)
 
   options = [];
-  switch (length (varargin))
-    case (0)
+  switch (numel (varargin))
+    case 0
       print_usage ();
     case {4,5}
-      if (length (varargin) == 4)
+      if (numel (varargin) == 4)
         [u, v, spx, spy] = varargin{:};
       else
         [u, v, spx, spy, options] = varargin{:};
       endif
       [m, n] = size (u);
       [x, y] = meshgrid (1:n, 1:m);
-    case (6)
+    case 6
       [x, y, u, v, spx, spy] = varargin{:};
-    case (7)
+    case 7
       [x, y, u, v, spx, spy, options] = varargin{:};
     otherwise
-      error ("stream2: unknown input parameter count");
+      error ("stream2: invalid number of inputs");
   endswitch
 
-  h = 0.1;
-  maxnverts = 10000;
+  stepsize = 0.1;
+  max_vertices = 10_000;
   if (! isempty (options))
-    switch (length (options))
-      case (1)
-        h = options(1);
-      case (2)
-        h = options(1);
-        maxnverts = options(2);
+    switch (numel (options))
+      case 1
+        stepsize = options(1);
+      case 2
+        stepsize = options(1);
+        max_vertices = options(2);
       otherwise
-        error ("stream2: wrong options length");
+        error ("stream2: invalid number of OPTIONS elements");
     endswitch
+
+    if (! isreal (stepsize) || stepsize == 0)
+      error ("stream2: STEPSIZE must be a real scalar != 0");
+    endif
+    if (! isreal (max_vertices) || max_vertices < 1)
+      error ("stream2: MAX_VERTICES must be an integer > 0");
+    endif
+    max_vertices = fix (max_vertices);
   endif
 
-  if ((! isnumeric (h)) || (h == 0))
-    error ("stream2: step size error");
-  endif
-  if ((! isnumeric (maxnverts)) || (maxnverts < 1))
-    error ("stream2: max num vertices error");
-  endif
-  if (! (isequal (size (u), size (v), size (x), size (y)) && ...
-         isequal (size (spx), size (spy))) )
+  if (! (size_equal (u, v, x, y) && size_equal (spx, spy)))
     error ("stream2: matrix dimensions must match");
   endif
-  if (iscomplex (u) || iscomplex (v) || iscomplex (x) || ...
-      iscomplex (y) || iscomplex (spx) || iscomplex (spy))
-    error ("stream2: input must be real valued");
+  if (iscomplex (u) || iscomplex (v) || iscomplex (x) || iscomplex (y)
+      || iscomplex (spx) || iscomplex (spy))
+    error ("stream2: all inputs must be real-valued");
   endif
 
   gx = x(1,:);
@@ -130,12 +130,12 @@
   dx = diff (gx);
   dy = diff (gy);
   ## "<" used to check if the mesh is ascending
-  if (any (dx <= 0) || any (dy <= 0) || ...
-      any (isnan (dx)) || any (isnan (dy)))
-    error ("stream2: ill shaped elements in mesh");
+  if (any (dx <= 0) || any (dy <= 0)
+      || any (isnan (dx)) || any (isnan (dy)))
+    error ("stream2: non-monotonically increasing or NaN values found in mesh");
   endif
-  tx = 1./dx;
-  ty = 1./dy;
+  tx = 1 ./ dx;
+  ty = 1 ./ dy;
   ## "Don't cares" used for handling points located on the border
   tx(end + 1) = 0;
   ty(end + 1) = 0;
@@ -145,15 +145,15 @@
   px = spx(:);
   py = spy(:);
 
-  for nseed = 1:length (px)
+  for nseed = 1 : numel (px)
 
     xp = px(nseed);
     yp = py(nseed);
-    idx = find (logical (diff (gx <= xp)), 1);
+    idx = find (diff (gx <= xp), 1);
     if (gx(end) == xp)
       idx = numel (gx);
     endif
-    idy = find (logical (diff (gy <= yp)), 1);
+    idy = find (diff (gy <= yp), 1);
     if (gy(end) == yp)
       idy = numel (gy);
     endif
@@ -165,7 +165,7 @@
       zeta = (idx - 1) + (xp - gx(idx)) * tx(idx);
       xi = (idy - 1) + (yp - gy(idy)) * ty(idy);
 
-      C = __streameuler2d__ (u, v, tx, ty, zeta, xi, h, maxnverts);
+      C = __streameuler2d__ (u, v, tx, ty, zeta, xi, stepsize, max_vertices);
 
       ## Transform from C coordinates to P coordinates
       idu = floor (C(:,1));
@@ -178,6 +178,7 @@
 
 endfunction
 
+
 %!demo
 %! clf;
 %! [x, y] = meshgrid (-5:5, -4:4);
@@ -195,5 +196,34 @@
 %! axis equal;
 
 %!test
-%! xy = stream2([1,1,1;2,2,2;3,3,3], [1,1,1;2,2,2;3,3,3], 1, 1, [0.01,5]);
+%! xy = stream2 ([1,1,1;2,2,2;3,3,3], [1,1,1;2,2,2;3,3,3], 1, 1, [0.01,5]);
 %! assert (numel (xy{:}), 10);
+
+## Test input validation
+%!error stream2 ()
+%!error <invalid number of inputs> stream2 (1)
+%!error <invalid number of inputs> stream2 (1,2)
+%!error <invalid number of inputs> stream2 (1,2,3)
+%!error <invalid number of OPTIONS> stream2 (1,2,3,4, [1,2,3])
+%!error <STEPSIZE must be a real scalar != 0> stream2 (1,2,3,4, [1i])
+%!error <STEPSIZE must be a real scalar != 0> stream2 (1,2,3,4, [0])
+%!error <MAX_VERTICES must be an integer> stream2 (1,2,3,4, [1, 1i])
+%!error <MAX_VERTICES must be an integer> stream2 (1,2,3,4, [1, 0])
+%!error <matrix dimensions must match> stream2 ([1 1],2,3,4)
+%!error <matrix dimensions must match> stream2 (1,[2 2],3,4)
+%!error <matrix dimensions must match> stream2 (1,2,[3 3],4)
+%!error <matrix dimensions must match> stream2 (1,2,3,[4 4])
+%!error <all inputs must be real-valued> stream2 (1i,2,3,4)
+%!error <all inputs must be real-valued> stream2 (1,2i,3,4)
+%!error <all inputs must be real-valued> stream2 (1,2,3i,4)
+%!error <all inputs must be real-valued> stream2 (1,2,3,4i)
+%!error <non-monotonically increasing or NaN values found in mesh>
+%! stream2 ([2 1], [1 2], [1 1], [2 2], [3 3], [4 4]);
+%!error <non-monotonically increasing or NaN values found in mesh>
+%! stream2 ([1 NaN], [1 2], [1 1], [2 2], [3 3], [4 4]);
+## FIXME: vectors representing x, y mesh are not accepted.
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream2 ([1 2], [2 1], [1 1], [2 2], [3 3], [4 4]);
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream2 ([1 2], [1 NaN], [1 1], [2 2], [3 3], [4 4]);
+

--- a/scripts/plot/draw/stream3.m	Thu Dec 12 20:35:09 2019 +0100
+++ b/scripts/plot/draw/stream3.m	Thu Dec 12 14:50:17 2019 -0800
@@ -19,7 +19,7 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {@var{xyz} =} stream3 (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
 ## @deftypefnx {} {@var{xyz} =} stream3 (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
-## @deftypefnx {} {@var{xyz} =} stream3 (@dots{}, "@var{options}")
+## @deftypefnx {} {@var{xyz} =} stream3 (@dots{}, @var{options})
 ## Compute 3-D streamline data.
 ##
 ## Calculate streamlines of a vector field given by @code{[@var{u}, @var{v},
@@ -30,9 +30,9 @@
 ## the vector field, @code{[]} is returned.
 ##
 ## The input parameter @var{options} is a 2-D vector of the form
-## @code{[@var{stepsize}, @var{maxnumbervertices}]}.  The first parameter
-## specifies the step size used for trajectory integration (default 0.1).  It
-## is allowed to set a negative value to control the direction of integration.
+## @code{[@var{stepsize}, @var{max_vertices}]}.  The first parameter
+## specifies the step size used for trajectory integration (default 0.1).  A
+## negative value is allowed which will reverse the direction of integration.
 ## The second parameter specifies the maximum number of segments used to
 ## create a streamline (default 10,000).
 ##
@@ -52,7 +52,6 @@
 ## @end example
 ##
 ## @seealso{stream2, streamline, streamtube}
-##
 ## @end deftypefn
 
 ## References:
@@ -76,53 +75,54 @@
 function xyz = stream3 (varargin)
 
   options = [];
-  switch (length (varargin))
-    case (0)
+  switch (numel (varargin))
+    case 0
       print_usage ();
     case {6,7}
-      if (length (varargin) == 6)
+      if (numel (varargin) == 6)
         [u, v, w, spx, spy, spz] = varargin{:};
       else
         [u, v, w, spx, spy, spz, options] = varargin{:};
       endif
       [m, n, p] = size (u);
       [x, y, z] = meshgrid (1:n, 1:m, 1:p);
-    case (9)
+    case 9
       [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
-    case (10)
+    case 10
       [x, y, z, u, v, w, spx, spy, spz, options] = varargin{:};
     otherwise
-      error ("stream3: unknown input parameter count");
+      error ("stream3: invalid number of inputs");
   endswitch
 
-  h = 0.1;
-  maxnverts = 10000;
+  stepsize = 0.1;
+  max_vertices = 10_000;
   if (! isempty (options))
-    switch (length (options))
-      case (1)
-        h = options(1);
-      case (2)
-        h = options(1);
-        maxnverts = options(2);
+    switch (numel (options))
+      case 1
+        stepsize = options(1);
+      case 2
+        stepsize = options(1);
+        max_vertices = options(2);
       otherwise
-        error ("stream3: wrong options length");
+        error ("stream3: invalid number of OPTIONS elements");
     endswitch
+
+    if (! isreal (stepsize) || stepsize == 0)
+      error ("stream2: STEPSIZE must be a real scalar != 0");
+    endif
+    if (! isreal (max_vertices) || max_vertices < 1)
+      error ("stream2: MAX_VERTICES must be an integer > 0");
+    endif
+    max_vertices = fix (max_vertices);
   endif
 
-  if ((! isnumeric (h)) || (h == 0))
-    error ("stream3: step size error");
-  endif
-  if ((! isnumeric (maxnverts)) || (maxnverts < 1))
-    error ("stream3: max num vertices error");
-  endif
-  if (! (isequal (size (u), size (v), size (w), size (x), size (y), size (z)) && ...
-         isequal (size (spx), size (spy))) )
+  if (! (size_equal (u, v, w, x, y, z) && size_equal (spx, spy, spz)))
     error ("stream3: matrix dimensions must match");
   endif
-  if (iscomplex (u) || iscomplex (v) || iscomplex (w) || ...
-      iscomplex (x) || iscomplex (y) || iscomplex (z) || ...
-      iscomplex (spx) || iscomplex (spy) || iscomplex (spz))
-    error ("stream3: input must be real valued");
+  if (iscomplex (u) || iscomplex (v) || iscomplex (w)
+      || iscomplex (x) || iscomplex (y) || iscomplex (z)
+      || iscomplex (spx) || iscomplex (spy) || iscomplex (spz))
+    error ("stream3: all inputs must be real-valued");
   endif
 
   gx = x(1, :, 1);
@@ -135,13 +135,13 @@
   dy = diff (gy);
   dz = diff (gz);
   ## "<" used to check if the mesh is ascending
-  if (any (dx <= 0) || any (dy <= 0) || any (dz <= 0) || ...
-      any (isnan (dx)) || any (isnan (dy)) || any (isnan (dz)))
-    error ("stream3: ill shaped elements in mesh");
+  if (any (dx <= 0) || any (dy <= 0) || any (dz <= 0)
+      || any (isnan (dx)) || any (isnan (dy)) || any (isnan (dz)))
+    error ("stream3: non-monotonically increasing or NaN values found in mesh");
   endif
-  tx = 1./dx;
-  ty = 1./dy;
-  tz = 1./dz;
+  tx = 1 ./ dx;
+  ty = 1 ./ dy;
+  tz = 1 ./ dz;
   ## "Don't cares" used for handling points located on the border
   tx(end + 1) = 0;
   ty(end + 1) = 0;
@@ -154,20 +154,20 @@
   py = spy(:);
   pz = spz(:);
 
-  for nseed = 1:length (px)
+  for nseed = 1 : numel (px)
 
     xp = px(nseed);
     yp = py(nseed);
     zp = pz(nseed);
-    idx = find (logical (diff (gx <= xp)), 1);
+    idx = find (diff (gx <= xp), 1);
     if (gx(end) == xp)
       idx = numel (gx);
     endif
-    idy = find (logical (diff (gy <= yp)), 1);
+    idy = find (diff (gy <= yp), 1);
     if (gy(end) == yp)
       idy = numel (gy);
     endif
-    idz = find (logical (diff (gz <= zp)), 1);
+    idz = find (diff (gz <= zp), 1);
     if (gz(end) == zp)
       idz = numel (gz);
     endif
@@ -180,7 +180,8 @@
       xi = (idy - 1) + (yp - gy(idy)) * ty(idy);
       rho = (idz - 1) + (zp - gz(idz)) * tz(idz);
 
-      C = __streameuler3d__ (u, v, w, tx, ty, tz, zeta, xi, rho, h, maxnverts);
+      C = __streameuler3d__ (u, v, w, tx, ty, tz, zeta, xi, rho, ...
+                             stepsize, max_vertices);
 
       ## Transform from C coordinates to P coordinates
       idu = floor (C(:, 1));
@@ -195,6 +196,7 @@
 
 endfunction
 
+
 %!demo
 %! clf;
 %! [x, y, z] = meshgrid (-30:1:30, -30:1:30, 0:1:50);
@@ -220,3 +222,42 @@
 %! [u, v, w] = meshgrid (0:3, 0:3, 0:3);
 %! xyz = stream3 (u, v, w, 2, 2, 2, [0.01,5]);
 %! assert (numel (xyz{:}), 15);
+
+## Test input validation
+%!error stream3 ()
+%!error <invalid number of inputs> stream3 (1)
+%!error <invalid number of inputs> stream3 (1,2)
+%!error <invalid number of inputs> stream3 (1,2,3)
+%!error <invalid number of inputs> stream3 (1,2,3,4)
+%!error <invalid number of inputs> stream3 (1,2,3,4,5)
+%!error <invalid number of inputs> stream3 (1,2,3,4,5,6,7,8)
+%!error <invalid number of OPTIONS> stream3 (1,2,3,4,5,6, [1,2,3])
+%!error <STEPSIZE must be a real scalar != 0> stream3 (1,2,3,4,5,6, [1i])
+%!error <STEPSIZE must be a real scalar != 0> stream3 (1,2,3,4,5,6, [0])
+%!error <MAX_VERTICES must be an integer> stream3 (1,2,3,4,5,6, [1, 1i])
+%!error <MAX_VERTICES must be an integer> stream3 (1,2,3,4,5,6, [1, 0])
+%!error <matrix dimensions must match> stream3 ([1 1],2,3,4,5,6)
+%!error <matrix dimensions must match> stream3 (1,[2 2],3,4,5,6)
+%!error <matrix dimensions must match> stream3 (1,2,[3 3],4,5,6)
+%!error <matrix dimensions must match> stream3 (1,2,3,[4 4],5,6)
+%!error <matrix dimensions must match> stream3 (1,2,3,4,[5 5],6)
+%!error <matrix dimensions must match> stream3 (1,2,3,4,5,[6 6])
+%!error <all inputs must be real-valued> stream3 (1i,2,3,4,5,6)
+%!error <all inputs must be real-valued> stream3 (1,2i,3,4,5,6)
+%!error <all inputs must be real-valued> stream3 (1,2,3i,4,5,6)
+%!error <all inputs must be real-valued> stream3 (1,2,3,4i,5,6)
+%!error <all inputs must be real-valued> stream3 (1,2,3,4,5i,6)
+%!error <all inputs must be real-valued> stream3 (1,2,3,4,5,6i)
+%!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([2 1], [1 2], [3 3], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);
+%!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([1 NaN], [1 2], [3 3], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);
+## FIXME: vectors representing x, y, z mesh are not accepted.
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([1 2], [2 1], [3 3], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([1 2], [1 NaN], [3 3], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([1 2], [1 2], [2 1], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);
+%#!error <non-monotonically increasing or NaN values found in mesh>
+%! stream3 ([1 2], [1 2], [1 NaN], [4 4], [5 5], [6 6], [7 7], [8 8], [9 9]);

--- a/scripts/plot/draw/streamline.m	Thu Dec 12 20:35:09 2019 +0100
+++ b/scripts/plot/draw/streamline.m	Thu Dec 12 14:50:17 2019 -0800
@@ -19,7 +19,7 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {} streamline (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
 ## @deftypefnx {} {} streamline (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
-## @deftypefnx {} {} streamline (@dots{}, "@var{options}")
+## @deftypefnx {} {} streamline (@dots{}, @var{options})
 ## @deftypefnx {} {} streamline (@var{hax}, @dots{})
 ## @deftypefnx {} {@var{h} =} streamline (@dots{})
 ## Plot streamlines of 2-D or 3-D vector fields.
@@ -31,9 +31,9 @@
 ## points @code{[@var{sx}, @var{sy}]} or @code{[@var{sx}, @var{sy}, @var{sz}]}.
 ##
 ## The input parameter @var{options} is a 2-D vector of the form
-## @code{[@var{stepsize}, @var{maxnumbervertices}]}.  The first parameter
-## specifies the step size used for trajectory integration (default 0.1).  It
-## is allowed to set a negative value to control the direction of integration.
+## @code{[@var{stepsize}, @var{max_vertices}]}.  The first parameter
+## specifies the step size used for trajectory integration (default 0.1).  A
+## negative value is allowed which will reverse the direction of integration.
 ## The second parameter specifies the maximum number of segments used to
 ## create a streamline (default 10,000).
 ##
@@ -55,16 +55,20 @@
 ## @end example
 ##
 ## @seealso{stream2, stream3, streamtube}
-##
 ## @end deftypefn
 
 function h = streamline (varargin)
 
+  [hax, varargin, nargin] = __plt_get_axis_arg__ ("streamline", varargin{:});
+
   if (nargin == 0)
     print_usage ();
   endif
 
-  [hax, varargin] = __plt_get_axis_arg__ ("streamline", varargin{:});
+  nd = ndims (varargin{1});
+  if (nd > 3)
+    error ("streamline: input data must be 2-D or 3-D");
+  endif
 
   if (isempty (hax))
     hax = gca ();
@@ -73,33 +77,31 @@
   endif
 
   h = [];
-  argval = varargin(1);
-  switch (numel (size (argval{:})))
-    case (2)
-      xy = stream2 (varargin{:});
-      for i = 1:length (xy)
-        sl = xy{i};
-        if (! isempty (sl))
-          tmp = line (hax, "xdata", sl(:, 1), "ydata", sl(:, 2), "color", "b");
-          h = [h; tmp];
-        endif
-      endfor
-    case (3)
-      xyz = stream3 (varargin{:});
-      for i = 1:length (xyz)
-        sl = xyz{i};
-        if (! isempty (sl))
-          tmp = line (hax, "xdata", sl(:, 1), "ydata", sl(:, 2), "zdata", sl(:, 3), ...
-                      "color", "b");
-          h = [h; tmp];
-        endif
-      endfor
-    otherwise
-      error ("streamline: 2D or 3D input data only");
-  endswitch
+  if (nd == 2)
+    xy = stream2 (varargin{:});
+    for i = 1 : length (xy)
+      sl = xy{i};
+      if (! isempty (sl))
+        htmp = line (hax, "xdata", sl(:, 1), "ydata", sl(:, 2), "color", "b");
+        h = [h; htmp];
+      endif
+    endfor
+  else
+    xyz = stream3 (varargin{:});
+    for i = 1 : length (xyz)
+      sl = xyz{i};
+      if (! isempty (sl))
+        htmp = line (hax,
+                     "xdata", sl(:, 1), "ydata", sl(:, 2), "zdata", sl(:, 3),
+                     "color", "b");
+        h = [h; htmp];
+      endif
+    endfor
+  endif
 
 endfunction
 
+
 %!demo
 %! clf;
 %! [x, y] = meshgrid (-2:0.5:2);
@@ -157,4 +159,14 @@
 %! grid on;
 %! axis equal tight;
 
+## Test input validation
 %!error streamline ()
+%!error <Invalid call to streamline>
+%! hf = figure ("visible", "off");
+%! unwind_protect
+%!   hax = axes ();
+%!   streamline (hax);
+%! unwind_protect_cleanup
+%!   close (hf);
+%! end_unwind_protect
+%!error <input data must be 2-D or 3-D> streamline (ones (2,2,2,2))

--- a/scripts/plot/draw/streamtube.m	Thu Dec 12 20:35:09 2019 +0100
+++ b/scripts/plot/draw/streamtube.m	Thu Dec 12 14:50:17 2019 -0800
@@ -20,7 +20,7 @@
 ## @deftypefn  {} {} streamtube (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
 ## @deftypefnx {} {} streamtube (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
 ## @deftypefnx {} {} streamtube (@var{vertices}, @var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w})
-## @deftypefnx {} {} streamtube (@dots{}, "@var{options}")
+## @deftypefnx {} {} streamtube (@dots{}, @var{options})
 ## @deftypefnx {} {} streamtube (@var{hax}, @dots{})
 ## @deftypefnx {} {@var{h} =} streamtube (@dots{})
 ## Calculate and display streamtubes.
@@ -34,14 +34,14 @@
 ## The streamtubes start at the seed points
 ## @code{[@var{sx}, @var{sy}, @var{sz}]}.
 ##
-## The tubes are colored depending on the local vector field strength.
+## The tubes are colored based on the local vector field strength.
 ##
 ## The input parameter @var{options} is a 2-D vector of the form
 ## @code{[@var{scale}, @var{n}]}.  The first parameter scales the start radius
 ## of the streamtubes (default 1).  The second parameter specifies the number
 ## of patches used for the streamtube circumference (default 20).
 ##
-## @code{streamtube} can be called with a cell array containing precomputed
+## @code{streamtube} can be called with a cell array containing pre-computed
 ## streamline data.  To do this, @var{vertices} must be created with the
 ## @code{stream3} function.  This option is useful if you need to alter the
 ## integrator step size or the maximum number of vertices of the streamline.
@@ -65,7 +65,6 @@
 ## @end example
 ##
 ## @seealso{stream3, streamline}
-##
 ## @end deftypefn
 
 ## References:
@@ -78,28 +77,18 @@
 
 function h = streamtube (varargin)
 
-  if (nargin == 0)
-    print_usage ();
-  endif
-
-  [hax, varargin] = __plt_get_axis_arg__ ("streamtube", varargin{:});
-
-  if (isempty (hax))
-    hax = gca ();
-  else
-    hax = hax(1);
-  endif
+  [hax, varargin, nargin] = __plt_get_axis_arg__ ("streamtube", varargin{:});
 
   options = [];
   xyz = [];
-  switch (length (varargin))
-    case (0)
+  switch (nargin)
+    case 0
       print_usage ();
-    case (6)
+    case 6
       [u, v, w, spx, spy, spz] = varargin{:};
       [m, n, p] = size (u);
       [x, y, z] = meshgrid (1:n, 1:m, 1:p);
-    case (7)
+    case 7
       if (iscell (varargin{1}))
         [xyz, x, y, z, u, v, w] = varargin{:};
       else
@@ -107,57 +96,64 @@
         [m, n, p] = size (u);
         [x, y, z] = meshgrid (1:n, 1:m, 1:p);
       endif
-    case (8)
+    case 8
       [xyz, x, y, z, u, v, w, options] = varargin{:};
-    case (9)
+    case 9
       [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
-    case (10)
+    case 10
       [x, y, z, u, v, w, spx, spy, spz, options] = varargin{:};
     otherwise
-      error ("streamtube: unknown input parameter count");
+      error ("streamtube: invalid number of inputs");
   endswitch
 
   scale = 1;
   n = 20;
   if (! isempty (options))
-    switch (length (options))
-      case (1)
+    switch (numel (options))
+      case 1
         scale = options(1);
-      case (2)
+      case 2
         scale = options(1);
         n = options(2);
       otherwise
-        error ("streamtube: wrong options length");
+        error ("streamtube: invalid number of OPTIONS elements");
     endswitch
+
+    if (! isreal (scale) || scale <= 0)
+      error ("streamtube: SCALE must be a real scalar > 0");
+    endif
+    if (! isreal (n) || n < 3)
+      error ("streamtube: number of polygons N must be greater than 2");
+    endif
+    n = fix (n);
   endif
 
-  if ((! isnumeric (scale)) || (scale <= 0))
-    error ("streamtube: scale error");
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
   endif
-  if ((! isnumeric (n)) || (n < 3))
-    error ("streamtube: number of polygons for tube circumference too small");
-  endif
+
   if isempty (xyz)
     xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.5);
   endif
 
   div = divergence (x, y, z, u, v, w);
   vn = sqrt (u.*u + v.*v + w.*w);
-  vmax = max (max (max (vn)));
-  vmin = min (min (min (vn)));
+  vmax = max (vn(:));
+  vmin = min (vn(:));
 
   ## Radius estimator
   [ny, nx, nz] = size (x);
-  dx = (max (max (max (x))) - min (min (min (x)))) / nx;
-  dy = (max (max (max (y))) - min (min (min (y)))) / ny;
-  dz = (max (max (max (z))) - min (min (min (z)))) / nz;
-  r0 = scale * sqrt (dx * dx + dy * dy + dz * dz);
+  dx = (max (x(:)) - min (x(:))) / nx;
+  dy = (max (y(:)) - min (y(:))) / ny;
+  dz = (max (z(:)) - min (z(:))) / nz;
+  r0 = scale * sqrt (dx*dx + dy*dy + dz*dz);
 
   h = [];
-  for i = 1:length (xyz)
-
+  for i = 1 : length (xyz)
     sl = xyz{i};
-    [nverts, ~] = size(sl);
+    nverts = rows (sl);
     if (! isempty (sl)) && (nverts > 2)
 
       divsl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
@@ -166,11 +162,10 @@
       wsl = interp3 (x, y, z, w, sl(:, 1), sl(:, 2), sl(:, 3));
       vv = sqrt (usl.*usl + vsl.*vsl + wsl.*wsl);
 
-      tmp = plottube (hax, sl, divsl, vv, vmax, vmin, r0, n);
-      h = [h, tmp];
+      htmp = plottube (hax, sl, divsl, vv, vmax, vmin, r0, n);
+      h = [h; htmp];
 
     endif
-
   endfor
 
 endfunction
@@ -184,7 +179,7 @@
     maxnverts = length (sl);
   endif
   if (maxnverts < 3)
-    error ("streamtube: too less data to show");
+    error ("streamtube: too little data to plot");
   endif
 
   if (vmax == vmin)
@@ -193,7 +188,7 @@
     colscale = 1.0 / (vmax - vmin);
   endif
 
-  phi = linspace (0, 2*pi (), npoly);
+  phi = linspace (0, 2*pi, npoly);
 
   X0 = sl(1, :);
   X1 = sl(2, :);
@@ -265,7 +260,7 @@
     N = [S(3), S(3), - S(1) - S(2)];
   endif
 
-  N = N / norm(N);
+  N /= norm (N);
 
 endfunction
 
@@ -273,7 +268,7 @@
 ## from starting point XS to ending point XE
 function [px, py, pz] = segment_patch_data (XS, XE)
 
-  [~, npoly] = size (XS);
+  npoly = columns (XS);
 
   px = zeros (4, npoly);
   py = zeros (4, npoly);
@@ -324,6 +319,7 @@
 
 endfunction
 
+
 %!demo
 %! clf;
 %! [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3.5:0.1:0);
@@ -368,5 +364,15 @@
 %! set (gca, "cameraviewanglemode", "manual");
 %! title ("Integration Towards Sink");
 
+## Test input validation
 %!error streamtube ()
-
+%!error <invalid number of inputs> streamtube (1)
+%!error <invalid number of inputs> streamtube (1,2)
+%!error <invalid number of inputs> streamtube (1,2,3)
+%!error <invalid number of inputs> streamtube (1,2,3,4)
+%!error <invalid number of inputs> streamtube (1,2,3,4,5)
+%!error <invalid number of OPTIONS> streamtube (1,2,3,4,5,6, [1,2,3])
+%!error <SCALE must be a real scalar . 0> streamtube (1,2,3,4,5,6, [1i])
+%!error <SCALE must be a real scalar . 0> streamtube (1,2,3,4,5,6, [0])
+%!error <N must be greater than 2> streamtube (1,2,3,4,5,6, [1, 1i])
+%!error <N must be greater than 2> streamtube (1,2,3,4,5,6, [1, 2])