--- a/doc/interpreter/plot.txi	Sun Feb 23 14:23:09 2020 +0100
+++ b/doc/interpreter/plot.txi	Tue Feb 25 19:26:33 2020 +0100
@@ -247,6 +247,8 @@
 
 @DOCSTRING(streamtube)
 
+@DOCSTRING(ostreamtube)
+
 @DOCSTRING(streamline)
 
 @DOCSTRING(stream2)

--- a/etc/NEWS.6	Sun Feb 23 14:23:09 2020 +0100
+++ b/etc/NEWS.6	Tue Feb 25 19:26:33 2020 +0100
@@ -206,6 +206,7 @@
 * `mustBeReal`
 * `namedargs2cell`
 * `newline`
+* `ostreamtube`
 * `rescale`
 * `rotx`
 * `roty`

--- a/libgui/src/gui-preferences-sc.h	Sun Feb 23 14:23:09 2020 +0100
+++ b/libgui/src/gui-preferences-sc.h	Tue Feb 25 19:26:33 2020 +0100
@@ -112,7 +112,7 @@
 const sc_pref sc_main_window_editor (sc_main_window + ":editor", PRE + CTRL + Qt::Key_4);
 const sc_pref sc_main_window_doc (sc_main_window + ":doc", PRE + CTRL + Qt::Key_5);
 const sc_pref sc_main_window_variable_editor (sc_main_window + ":variable_editor", PRE + CTRL + Qt::Key_6);
-const sc_pref sc_main_window_previous_dock (sc_main_window + ":previous_widget", PRE + CTRL_ALT + Qt::Key_0);
+const sc_pref sc_main_window_previous_dock (sc_main_window + ":previous_widget", PRE + CTRL_ALT + Qt::Key_P);
 const sc_pref sc_main_window_reset (sc_main_window + ":reset", QKeySequence::UnknownKey);
 
 // help

--- a/scripts/plot/draw/module.mk	Sun Feb 23 14:23:09 2020 +0100
+++ b/scripts/plot/draw/module.mk	Tue Feb 25 19:26:33 2020 +0100
@@ -61,6 +61,7 @@
   %reldir%/mesh.m \
   %reldir%/meshc.m \
   %reldir%/meshz.m \
+  %reldir%/ostreamtube.m \
   %reldir%/pareto.m \
   %reldir%/patch.m \
   %reldir%/pcolor.m \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/ostreamtube.m	Tue Feb 25 19:26:33 2020 +0100
@@ -0,0 +1,366 @@
+########################################################################
+##
+## Copyright (C) 2019-2020 The Octave Project Developers
+##
+## See the file COPYRIGHT.md in the top-level directory of this
+## distribution or <https://octave.org/copyright/>.
+##
+## 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
+## <https://www.gnu.org/licenses/>.
+##
+########################################################################
+
+## -*- texinfo -*-
+## @deftypefn  {} {} ostreamtube (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
+## @deftypefnx {} {} ostreamtube (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
+## @deftypefnx {} {} ostreamtube (@var{xyz}, @var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w})
+## @deftypefnx {} {} ostreamtube (@dots{}, @var{options})
+## @deftypefnx {} {} ostreamtube (@var{hax}, @dots{})
+## @deftypefnx {} {@var{h} =} ostreamtube (@dots{})
+## Calculate and display streamtubes.
+##
+## Streamtubes are approximated by connecting circular crossflow areas
+## along a streamline.  The expansion of the flow is determined by the local
+## crossflow divergence.
+##
+## The vector field is given by @code{[@var{u}, @var{v}, @var{w}]} and is
+## defined over a rectangular grid given by @code{[@var{x}, @var{y}, @var{z}]}.
+## The streamtubes start at the seed points
+## @code{[@var{sx}, @var{sy}, @var{sz}]}.
+##
+## 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 vertices that are used to construct the tube circumference (default 20).
+##
+## @code{ostreamtube} can be called with a cell array containing pre-computed
+## streamline data.  To do this, @var{xyz} 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.
+##
+## If the first argument @var{hax} is an axes handle, then plot into this axes,
+## rather than the current axes returned by @code{gca}.
+##
+## The optional return value @var{h} is a graphics handle to the plot
+## objects created for each streamtube.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3:0.1:0);
+## u = -x / 10 - y;
+## v = x - y / 10;
+## w = - ones (size (x)) / 10;
+## ostreamtube (x, y, z, u, v, w, 1, 0, 0);
+## @end group
+## @end example
+##
+## @seealso{stream3, streamline, streamtube}
+## @end deftypefn
+
+## References:
+##
+## @inproceedings{
+##    title = {Visualization of 3-D vector fields - Variations on a stream},
+##    author = {Dave Darmofal and Robert Haimes},
+##    year = {1992}
+## }
+##
+## @article{
+##    title = {Efficient streamline, streamribbon, and streamtube constructions on unstructured grids},
+##    author = {Ueng, Shyh-Kuang and Sikorski, C. and Ma, Kwan-Liu},
+##    year = {1996},
+##    month = {June},
+##    publisher = {IEEE Transactions on Visualization and Computer Graphics},
+## }
+
+function h = ostreamtube (varargin)
+
+  [hax, varargin, nargin] = __plt_get_axis_arg__ ("ostreamtube", varargin{:});
+
+  options = [];
+  xyz = [];
+  switch (nargin)
+    case 0
+      print_usage ();
+    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
+      if (iscell (varargin{1}))
+        [xyz, x, y, z, u, v, w] = varargin{:};
+      else
+        [u, v, w, spx, spy, spz, options] = varargin{:};
+        [m, n, p] = size (u);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 8
+      [xyz, x, y, z, u, v, w, options] = varargin{:};
+    case 9
+      [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
+    case 10
+      [x, y, z, u, v, w, spx, spy, spz, options] = varargin{:};
+    otherwise
+      error ("ostreamtube: invalid number of inputs");
+  endswitch
+
+  scale = 1;
+  num_circum = 20;
+  if (! isempty (options))
+    switch (numel (options))
+      case 1
+        scale = options(1);
+      case 2
+        scale = options(1);
+        num_circum = options(2);
+      otherwise
+        error ("ostreamtube: invalid number of OPTIONS elements");
+    endswitch
+
+    if (! isreal (scale) || scale <= 0)
+      error ("ostreamtube: SCALE must be a real scalar > 0");
+    endif
+    if (! isreal (num_circum) || num_circum < 3)
+      error ("ostreamtube: number of tube vertices N must be greater than 2");
+    endif
+    num_circum = fix (num_circum);
+  endif
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  endif
+
+  if (isempty (xyz))
+    xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.2);
+  endif
+
+  div = divergence (x, y, z, u, v, w);
+
+  ## Use the bounding box diagonal to determine the starting radius
+  mxx = mnx = mxy = mny = mxz = mnz = [];
+  j = 1;
+  for i = 1 : length (xyz)
+    sl = xyz{i};
+    if (! isempty (sl))
+      slx = sl(:, 1); sly = sl(:, 2); slz = sl(:, 3);
+      mxx(j) = max (slx); mnx(j) = min (slx);
+      mxy(j) = max (sly); mny(j) = min (sly);
+      mxz(j) = max (slz); mnz(j) = min (slz);
+      j += 1;
+    endif
+  endfor
+  dx = max (mxx) - min (mnx);
+  dy = max (mxy) - min (mny);
+  dz = max (mxz) - min (mnz);
+  rstart = scale * sqrt (dx*dx + dy*dy + dz*dz) / 25;
+
+  h = [];
+  for i = 1 : length (xyz)
+    sl = xyz{i};
+    num_vertices = rows (sl);
+    if (! isempty (sl) && num_vertices > 2)
+
+      usl = interp3 (x, y, z, u, sl(:, 1), sl(:, 2), sl(:, 3));
+      vsl = interp3 (x, y, z, v, sl(:, 1), sl(:, 2), sl(:, 3));
+      wsl = interp3 (x, y, z, w, sl(:, 1), sl(:, 2), sl(:, 3));
+      vv = sqrt (usl.*usl + vsl.*vsl + wsl.*wsl);
+
+      div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
+      is_singular_div = find (isnan (div_sl), 1, "first");
+
+      if (! isempty (is_singular_div))
+        max_vertices = is_singular_div - 1;
+      else
+        max_vertices = num_vertices;
+      endif
+
+      if (max_vertices > 2)
+
+        htmp = plottube (hax, sl, div_sl, vv, max_vertices, ...
+                         rstart, num_circum);
+        h = [h; htmp];
+
+      endif
+    endif
+  endfor
+
+endfunction
+
+function h = plottube (hax, sl, div_sl, vv, max_vertices, rstart, num_circum)
+
+  phi = linspace (0, 2*pi, num_circum);
+  cp = cos (phi);
+  sp = sin (phi);
+
+  X0 = sl(1, :);
+  X1 = sl(2, :);
+
+  ## 1st rotation axis
+  R = X1 - X0;
+  RE = R / norm (R);
+
+  ## Initial radius
+  vold = vv(1);
+  vact = vv(2);
+  ract = rstart * exp (0.5 * div_sl(2) * norm (R) / vact) * sqrt (vold / vact);
+  vold = vact;
+  rold = ract;
+
+  ## Guide point and its rotation to create a segment
+  N = get_normal1 (R);
+  K = ract * N;
+  XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
+
+  px = zeros (num_circum, max_vertices - 1);
+  py = zeros (num_circum, max_vertices - 1);
+  pz = zeros (num_circum, max_vertices - 1);
+  pc = zeros (num_circum, max_vertices - 1);
+
+  px(:, 1) = XS(1, :).';
+  py(:, 1) = XS(2, :).';
+  pz(:, 1) = XS(3, :).';
+  pc(:, 1) = vact * ones (num_circum, 1);
+
+  for i = 3 : max_vertices
+
+    KK = K;
+    X0 = X1;
+    X1 = sl(i, :);
+    R = X1 - X0;
+    RE = R / norm (R);
+
+    ## Tube radius
+    vact = vv(i);
+    ract = rold * exp (0.5 * div_sl(i) * norm (R) / vact) * sqrt (vold / vact);
+    vold = vact;
+    rold = ract;
+
+    ## Project KK onto RE and get the difference in order to calculate the next
+    ## guiding point
+    Kp = KK - RE * dot (KK, RE);
+    K = ract * Kp / norm (Kp);
+
+    ## Rotate around RE and collect surface patches
+    XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
+
+    px(:, i - 1) = XS(1, :).';
+    py(:, i - 1) = XS(2, :).';
+    pz(:, i - 1) = XS(3, :).';
+    pc(:, i - 1) = vact * ones (num_circum, 1);
+
+  endfor
+
+  h = surface (hax, px, py, pz, pc);
+
+endfunction
+
+## Arbitrary N normal to X
+function N = get_normal1 (X)
+
+  if ((X(3) == 0) && (X(1) == -X(2)))
+    N = [- X(2) - X(3), X(1), X(1)];
+  else
+    N = [X(3), X(3), - X(1) - X(2)];
+  endif
+
+  N /= norm (N);
+
+endfunction
+
+## Rotate X around U where |U| = 1
+## cp = cos (angle), sp = sin (angle)
+function Y = rotation (X, U, cp, sp)
+
+  ux = U(1);
+  uy = U(2);
+  uz = U(3);
+
+  Y(1, :) = X(1) * (cp + ux * ux * (1 - cp)) + ...
+            X(2) * (ux * uy * (1 - cp) - uz * sp) + ...
+            X(3) * (ux * uz * (1 - cp) + uy * sp);
+
+  Y(2, :) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ...
+            X(2) * (cp + uy * uy * (1 - cp)) + ...
+            X(3) * (uy * uz * (1 - cp) - ux * sp);
+
+  Y(3, :) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ...
+            X(2) * (uz * uy * (1 - cp) + ux * sp) + ...
+            X(3) * (cp + uz * uz * (1 - cp));
+
+endfunction
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3.5:0.1:0);
+%! a = 0.1;
+%! b = 0.1;
+%! u = - a * x - y;
+%! v = x - a * y;
+%! w = - b * ones (size (x));
+%! sx = 1.0;
+%! sy = 0.0;
+%! sz = 0.0;
+%! ostreamtube (x, y, z, u, v, w, sx, sy, sz, [1.2, 30]);
+%! colormap (jet);
+%! shading interp;
+%! view ([-47, 24]);
+%! camlight ();
+%! lighting gouraud;
+%! grid on;
+%! view (3);
+%! axis equal;
+%! set (gca, "cameraviewanglemode", "manual");
+%! title ("Spiral Sink");
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-2:0.5:2);
+%! t = sqrt (1.0./(x.^2 + y.^2 + z.^2)).^3;
+%! u = - x.*t;
+%! v = - y.*t;
+%! w = - z.*t;
+%! [sx, sy, sz] = meshgrid (-2:4:2);
+%! xyz = stream3 (x, y, z, u, v, w, sx, sy, sz, [0.1, 60]);
+%! ostreamtube (xyz, x, y, z, u, v, w, [2, 50]);
+%! colormap (jet);
+%! shading interp;
+%! view ([-47, 24]);
+%! camlight ();
+%! lighting gouraud;
+%! grid on;
+%! view (3);
+%! axis equal;
+%! set (gca, "cameraviewanglemode", "manual");
+%! title ("Integration Towards Sink");
+
+## Test input validation
+%!error ostreamtube ()
+%!error <invalid number of inputs> ostreamtube (1)
+%!error <invalid number of inputs> ostreamtube (1,2)
+%!error <invalid number of inputs> ostreamtube (1,2,3)
+%!error <invalid number of inputs> ostreamtube (1,2,3,4)
+%!error <invalid number of inputs> ostreamtube (1,2,3,4,5)
+%!error <invalid number of OPTIONS> ostreamtube (1,2,3,4,5,6, [1,2,3])
+%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [1i])
+%!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [0])
+%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 1i])
+%!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 2])

--- a/scripts/plot/draw/stream3.m	Sun Feb 23 14:23:09 2020 +0100
+++ b/scripts/plot/draw/stream3.m	Tue Feb 25 19:26:33 2020 +0100
@@ -58,7 +58,7 @@
 ## @end group
 ## @end example
 ##
-## @seealso{stream2, streamline, streamtube}
+## @seealso{stream2, streamline, streamtube, ostreamtube}
 ## @end deftypefn
 
 ## References:

--- a/scripts/plot/draw/streamline.m	Sun Feb 23 14:23:09 2020 +0100
+++ b/scripts/plot/draw/streamline.m	Tue Feb 25 19:26:33 2020 +0100
@@ -61,7 +61,7 @@
 ## @end group
 ## @end example
 ##
-## @seealso{stream2, stream3, streamtube}
+## @seealso{stream2, stream3, streamtube, ostreamtube}
 ## @end deftypefn
 
 function h = streamline (varargin)

--- a/scripts/plot/draw/streamtube.m	Sun Feb 23 14:23:09 2020 +0100
+++ b/scripts/plot/draw/streamtube.m	Tue Feb 25 19:26:33 2020 +0100
@@ -26,85 +26,91 @@
 ## -*- texinfo -*-
 ## @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 (@var{xyz}, @var{x}, @var{y}, @var{z}, @var{div})
+## @deftypefnx {} {} streamtube (@var{xyz}, @var{div})
+## @deftypefnx {} {} streamtube (@var{xyz}, @var{dia})
 ## @deftypefnx {} {} streamtube (@dots{}, @var{options})
 ## @deftypefnx {} {} streamtube (@var{hax}, @dots{})
 ## @deftypefnx {} {@var{h} =} streamtube (@dots{})
-## Calculate and display streamtubes.
+## Plot tubes scaled by the divergence along streamlines.
 ##
-## Streamtubes are approximated by connecting circular crossflow areas
-## along a streamline.  The expansion of the flow is determined by the local
-## crossflow divergence.
+## @code{streamtube} draws tubes whose diameter is scaled by the divergence of
+## a vector field.  The vector field is given by
+## @code{[@var{u}, @var{v}, @var{w}]} and is defined over a rectangular grid
+## given by @code{[@var{x}, @var{y}, @var{z}]}.  The tubes start at the
+## seed points @code{[@var{sx}, @var{sy}, @var{sz}]} and are plot along
+## streamlines.
 ##
-## The vector field is given by @code{[@var{u}, @var{v}, @var{w}]} and is
-## defined over a rectangular grid given by @code{[@var{x}, @var{y}, @var{z}]}.
-## The streamtubes start at the seed points
-## @code{[@var{sx}, @var{sy}, @var{sz}]}.
+## @code{streamtube} can also be called with a cell array containing
+## pre-computed streamline data.  To do this, @var{xyz} must be created with
+## the @code{stream3} command.  @var{div} is used to scale the tubes.
+## In order to plot tubes scaled by the vector field divergence, @var{div}
+## must be calculated with the @code{divergence} command.
 ##
-## The tubes are colored based on the local vector field strength.
+## A tube diameter of zero corresponds to the smallest scaling value along the
+## streamline and the largest tube diameter corresponds to the largest scaling
+## value.
+##
+## It is also possible to draw a tube along an arbitrary array of vertices
+## @var{xyz}.  The tube diameter can be specified by the vertex array @var{dia}
+## or by a constant.
 ##
 ## 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 polygon points used for the streamtube circumference (default 20).
-##
-## @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.
+## @code{[@var{scale}, @var{n}]}.  The first parameter scales the tube
+## diameter (default 1).  The second parameter specifies the number of vertices
+## that are used to construct the tube circumference (default 20).
 ##
 ## If the first argument @var{hax} is an axes handle, then plot into this axes,
 ## rather than the current axes returned by @code{gca}.
 ##
-## The optional return value @var{h} is a graphics handle to the patch plot
-## objects created for each streamtube.
-##
-## Example:
+## The optional return value @var{h} is a graphics handle to the plot objects
+## created for each tube.
 ##
-## @example
-## @group
-## [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3:0.1:0);
-## u = -x / 10 - y;
-## v = x - y / 10;
-## w = - ones (size (x)) / 10;
-## streamtube (x, y, z, u, v, w, 1, 0, 0);
-## @end group
-## @end example
-##
-## @seealso{stream3, streamline}
+## @seealso{stream3, streamline, ostreamtube}
 ## @end deftypefn
 
-## Reference:
-##
-## @inproceedings{
-##    title = {Visualization of 3-D vector fields - Variations on a stream},
-##    author = {Dave Darmofal and Robert Haimes},
-##    year = {1992}
-## }
-
 function h = streamtube (varargin)
 
   [hax, varargin, nargin] = __plt_get_axis_arg__ ("streamtube", varargin{:});
 
   options = [];
   xyz = [];
+  div = [];
+  dia = [];
   switch (nargin)
     case 0
       print_usage ();
+    case 2
+      ## "dia" can be a cell array or a constant
+      if (iscell (varargin{2}) || numel (varargin{2}) == 1)
+        [xyz, dia] = varargin{:};
+      else
+        [xyz, div] = varargin{:};
+        [m, n, p] = size (div);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 3
+      if (iscell (varargin{2}))
+        [xyz, dia, options] = varargin{:};
+      else
+        [xyz, div, options] = varargin{:};
+        [m, n, p] = size (div);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 5
+        [xyz, x, y, z, div] = varargin{:};
     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
       if (iscell (varargin{1}))
-        [xyz, x, y, z, u, v, w] = varargin{:};
+        [xyz, x, y, z, div, options] = varargin{:};
       else
-        [u, v, w, spx, spy, spz, options] = varargin{:};
+        [u, v, w, spx, spy, spz] = varargin{:};
         [m, n, p] = size (u);
         [x, y, z] = meshgrid (1:n, 1:m, 1:p);
       endif
-    case 8
-      [xyz, x, y, z, u, v, w, options] = varargin{:};
+    case 7
+      [u, v, w, spx, spy, spz, options] = varargin{:};
+      [m, n, p] = size (u);
+      [x, y, z] = meshgrid (1:n, 1:m, 1:p);
     case 9
       [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
     case 10
@@ -114,14 +120,14 @@
   endswitch
 
   scale = 1;
-  num_poly = 20;
+  num_circum = 20;
   if (! isempty (options))
     switch (numel (options))
       case 1
         scale = options(1);
       case 2
         scale = options(1);
-        num_poly = options(2);
+        num_circum = options(2);
       otherwise
         error ("streamtube: invalid number of OPTIONS elements");
     endswitch
@@ -129,10 +135,26 @@
     if (! isreal (scale) || scale <= 0)
       error ("streamtube: SCALE must be a real scalar > 0");
     endif
-    if (! isreal (num_poly) || num_poly < 3)
-      error ("streamtube: number of polygons N must be greater than 2");
+    if (! isreal (num_circum) || num_circum < 3)
+      error ("streamtube: number of tube vertices N must be greater than 2");
     endif
-    num_poly = fix (num_poly);
+    num_circum = fix (num_circum);
+  endif
+
+  if (isempty (xyz))
+    xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.2);
+  endif
+
+  if (isempty (div) && isempty (dia))
+    div = divergence (x, y, z, u, v, w);
+  endif
+
+  if (! isempty (dia) && iscell (dia))
+    for i = 1 : length (xyz)
+      if (rows (dia{i}) != rows (xyz{i}))
+        error ("streamtube: DIA must have same length then XYZ");
+      endif
+    endfor
   endif
 
   if (isempty (hax))
@@ -141,54 +163,55 @@
     hax = hax(1);
   endif
 
-  if (isempty (xyz))
-    xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.2);
-  endif
-
-  div = divergence (x, y, z, u, v, w);
-
-  ## Determine start radius
+  ## Derive final scale factor from the bounding box diagonal
   mxx = mnx = mxy = mny = mxz = mnz = [];
   j = 1;
   for i = 1 : length (xyz)
     sl = xyz{i};
     if (! isempty (sl))
       slx = sl(:, 1); sly = sl(:, 2); slz = sl(:, 3);
-      mxx(j) = max (slx(:)); mnx(j) = min (slx(:));
-      mxy(j) = max (sly(:)); mny(j) = min (sly(:));
-      mxz(j) = max (slz(:)); mnz(j) = min (slz(:));
+      mxx(j) = max (slx); mnx(j) = min (slx);
+      mxy(j) = max (sly); mny(j) = min (sly);
+      mxz(j) = max (slz); mnz(j) = min (slz);
       j += 1;
     endif
   endfor
   dx = max (mxx) - min (mnx);
   dy = max (mxy) - min (mny);
   dz = max (mxz) - min (mnz);
-  rstart = scale * sqrt (dx*dx + dy*dy + dz*dz) / 25;
+  clen = scale * sqrt (dx*dx + dy*dy + dz*dz) / 40;
 
   h = [];
   for i = 1 : length (xyz)
     sl = xyz{i};
     num_vertices = rows (sl);
-    if (! isempty (sl) && num_vertices > 2)
-
-      usl = interp3 (x, y, z, u, sl(:, 1), sl(:, 2), sl(:, 3));
-      vsl = interp3 (x, y, z, v, sl(:, 1), sl(:, 2), sl(:, 3));
-      wsl = interp3 (x, y, z, w, sl(:, 1), sl(:, 2), sl(:, 3));
-      vv = sqrt (usl.*usl + vsl.*vsl + wsl.*wsl);
-
-      div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
-      is_singular_div = find (isnan (div_sl), 1, "first");
+    if (! isempty (sl) && num_vertices > 1)
+      if (isempty (dia))
 
-      if (! isempty (is_singular_div))
-        max_vertices = is_singular_div - 1;
+        ## Plot a tube based on normalized divergence
+        [div_sl, max_vertices] = interp_sl (x, y, z, div, sl);
+        if (max_vertices > 1)
+          ## Nomalize the divergence along the streamline
+          mn = min (div_sl);
+          mx = max (div_sl);
+          if (mn == mx)
+            radius_sl = clen * ones (max_vertices, 1);
+          else
+            radius_sl = clen * (div_sl - mn) / (mx - mn);
+          endif
+          htmp = plottube (hax, sl, radius_sl, max_vertices, num_circum);
+          h = [h; htmp];
+        endif
+
       else
-        max_vertices = num_vertices;
-      endif
 
-      if (max_vertices > 2)
-
-        htmp = plottube (hax, sl, div_sl, vv, max_vertices, ...
-                         rstart, num_poly);
+        ## Plot a tube from external data (vertex array or constant)
+        if (iscell (dia))
+          radius_sl = 0.5 * scale * dia{i};
+        else
+          radius_sl = 0.5 * scale * dia * ones (1, num_vertices);
+        endif
+        htmp = plottube (hax, sl, radius_sl, num_vertices, num_circum);
         h = [h; htmp];
 
       endif
@@ -197,9 +220,9 @@
 
 endfunction
 
-function h = plottube (hax, sl, div_sl, vv, max_vertices, rstart, num_poly)
+function h = plottube (hax, sl, radius_sl, max_vertices, num_circum)
 
-  phi = linspace (0, 2*pi, num_poly);
+  phi = linspace (0, 2*pi, num_circum);
   cp = cos (phi);
   sp = sin (phi);
 
@@ -210,58 +233,71 @@
   R = X1 - X0;
   RE = R / norm (R);
 
-  ## Initial radius
-  vold = vv(1);
-  vact = vv(2);
-  ract = rstart * exp (0.5 * div_sl(2) * norm (R) / vact) * sqrt (vold / vact);
-  vold = vact;
-  rold = ract;
+  ## Guide point and its rotation to create a segment
+  KE = get_normal1 (RE);
+  K = radius_sl(1) * KE;
+  XS0 = rotation (K, RE, cp, sp) + repmat (X0.', 1, num_circum);
+  K = radius_sl(2) * KE;
+  XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
 
-  ## Guide point and its rotation to create a segment
-  N = get_normal1 (R);
-  K = ract * N;
-  XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_poly);
+  px = zeros (num_circum, max_vertices);
+  py = zeros (num_circum, max_vertices);
+  pz = zeros (num_circum, max_vertices);
 
-  px = zeros (num_poly, max_vertices - 1);
-  py = zeros (num_poly, max_vertices - 1);
-  pz = zeros (num_poly, max_vertices - 1);
-  pc = zeros (num_poly, max_vertices - 1);
+  px(:, 1) = XS0(1, :).';
+  py(:, 1) = XS0(2, :).';
+  pz(:, 1) = XS0(3, :).';
 
-  px(:, 1) = XS(1, :).';
-  py(:, 1) = XS(2, :).';
-  pz(:, 1) = XS(3, :).';
-  pc(:, 1) = vact * ones (num_poly, 1);
+  px(:, 2) = XS(1, :).';
+  py(:, 2) = XS(2, :).';
+  pz(:, 2) = XS(3, :).';
 
   for i = 3 : max_vertices
 
-    KK = K;
+    KEold = KE;
     X0 = X1;
+
     X1 = sl(i, :);
     R = X1 - X0;
     RE = R / norm (R);
 
-    ## Tube radius
-    vact = vv(i);
-    ract = rold * exp (0.5 * div_sl(i) * norm (R) / vact) * sqrt (vold / vact);
-    vold = vact;
-    rold = ract;
-
-    ## Project K onto RE and get the difference in order to calculate the next
+    ## Project KE onto RE and get the difference in order to calculate the next
     ## guiding point
-    Kp = KK - RE * dot (KK, RE);
-    K = ract * Kp / norm (Kp);
+    Kp = KEold - RE * dot (KEold, RE);
+    KE = Kp / norm (Kp);
+    K = radius_sl(i) * KE;
 
-    ## Rotate the guiding point around R and collect patch vertices
-    XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_poly);
+    ## Rotate around RE and collect surface patches
+    XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
 
-    px(:, i - 1) = XS(1, :).';
-    py(:, i - 1) = XS(2, :).';
-    pz(:, i - 1) = XS(3, :).';
-    pc(:, i - 1) = vact * ones (num_poly, 1);
+    px(:, i) = XS(1, :).';
+    py(:, i) = XS(2, :).';
+    pz(:, i) = XS(3, :).';
 
   endfor
 
-  h = surface (hax, px, py, pz, pc);
+  h = surface (hax, px, py, pz);
+
+endfunction
+
+## Interpolate onto the streamline vertices and return the first chunck of
+## valid samples until a singularity is hit (NaN or +-Inf) or
+## the streamline vertex array "sl" ends
+function [div_sl_crop, max_vertices] = interp_sl (x, y, z, div, sl)
+
+  div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
+  is_nan = find (isnan (div_sl), 1, "first");
+  is_inf = find (isinf (div_sl), 1, "first");
+
+  max_vertices = rows (sl);
+  if (! isempty (is_nan))
+    max_vertices = min (max_vertices, is_nan - 1);
+  endif
+  if (! isempty (is_inf))
+    max_vertices = min (max_vertices, is_inf - 1);
+  endif
+
+  div_sl_crop = div_sl(1 : max_vertices);
 
 endfunction
 
@@ -300,60 +336,48 @@
 
 endfunction
 
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-3:0.15:3, -1:0.1:1, -1:0.1:1);
+%! u = 2 + 8 * exp (-2.0*x.*x);
+%! v = zeros (size (x));
+%! w = zeros (size (x));
+%! h = streamtube (x, y, z, u, v, w, -3, 0, 0, [5, 60]);
+%! set (h, "facecolor", "r", "edgecolor", "none");
+%! hold on;
+%! camlight ();
+%! lighting gouraud;
+%! view (3);
+%! grid on;
+%! quiver3 (x, y, z, u, v, w);
+%! axis tight equal;
+%! title ("Divergence Plot");
 
 %!demo
 %! clf;
-%! [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3.5:0.1:0);
-%! a = 0.1;
-%! b = 0.1;
-%! u = - a * x - y;
-%! v = x - a * y;
-%! w = - b * ones (size (x));
-%! sx = 1.0;
-%! sy = 0.0;
-%! sz = 0.0;
-%! streamtube (x, y, z, u, v, w, sx, sy, sz, [1.2, 30]);
+%! t = 0:.15:15;
+%! xyz{1} = [cos(t)' sin(t)' (t/3)'];
+%! dia{1} = cos(t)';
+%! streamtube (xyz, dia);
+%! grid on;
+%! axis tight equal;
 %! colormap (jet);
 %! shading interp;
-%! view ([-47, 24]);
 %! camlight ();
 %! lighting gouraud;
-%! grid on;
 %! view (3);
-%! axis equal;
-%! set (gca, "cameraviewanglemode", "manual");
-%! title ("Spiral Sink");
-
-%!demo
-%! clf;
-%! [x, y, z] = meshgrid (-2:0.5:2);
-%! t = sqrt (1.0./(x.^2 + y.^2 + z.^2)).^3;
-%! u = - x.*t;
-%! v = - y.*t;
-%! w = - z.*t;
-%! [sx, sy, sz] = meshgrid (-2:4:2);
-%! xyz = stream3 (x, y, z, u, v, w, sx, sy, sz, [0.1, 60]);
-%! streamtube (xyz, x, y, z, u, v, w, [2, 50]);
-%! colormap (jet);
-%! shading interp;
-%! view ([-47, 24]);
-%! camlight ();
-%! lighting gouraud;
-%! grid on;
-%! view (3);
-%! axis equal;
-%! set (gca, "cameraviewanglemode", "manual");
-%! title ("Integration Towards Sink");
+%! title ("Plot Arbitrary Tube");
 
 ## 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])
+%!error <invalid number of inputs> streamtube (1,2,3,4,5,6,7,8)
+%!error <invalid number of inputs> streamtube (1,2,3,4,5,6,7,8,9,10,11)
+%!error <invalid number of OPTIONS elements> streamtube (1,2,[1,2,3])
+%!error <SCALE must be a real scalar . 0> streamtube (1,2,[1i])
+%!error <SCALE must be a real scalar . 0> streamtube (1,2,[0])
+%!error <SCALE must be a real scalar . 0> streamtube (1,2,[-1])
+%!error <N must be greater than 2> streamtube (1,2,[1,1i])
+%!error <N must be greater than 2> streamtube (1,2,[1,2])
+%!error <DIA must have same length then XYZ> streamtube ({[1,1,1;2,2,2]},{[1,1,1]})

--- a/src/mkoctfile.in.cc	Sun Feb 23 14:23:09 2020 +0100
+++ b/src/mkoctfile.in.cc	Tue Feb 25 19:26:33 2020 +0100
@@ -792,7 +792,7 @@
       else if (starts_with (arg, "-Wl,") || starts_with (arg, "-l")
                || starts_with (arg, "-L") || starts_with (arg, "-R"))
         {
-          ldflags += (' ' + arg);
+          ldflags += (' ' + quote_path (arg));
         }
 #if ! defined (OCTAVE_USE_WINDOWS_API)
       else if (arg == "-pthread")