# HG changeset patch
# User HG200
# Date 1574506036 -3600
#      Sat Nov 23 11:47:16 2019 +0100
# Branch stable
# Node ID fcae01084f8fc667f2d621ee3a6c84e7fd3142c8
# Parent  f998e243fa785af25b8b2ab63eea5cb018690ac3
branch 'stable'
added libinterp/dldfcn/stream_euler.cc
added scripts/plot/draw/stream2.m
added scripts/plot/draw/stream3.m
added scripts/plot/draw/streamline.m
added scripts/plot/draw/streamtube.m
changed doc/interpreter/plot.txi
changed libinterp/dldfcn/module-files
changed scripts/plot/draw/module.mk

diff -r f998e243fa78 -r fcae01084f8f doc/interpreter/plot.txi
--- a/doc/interpreter/plot.txi	Wed Nov 13 11:11:49 2019 +0100
+++ b/doc/interpreter/plot.txi	Sat Nov 23 11:47:16 2019 +0100
@@ -245,6 +245,14 @@
 
 @DOCSTRING(quiver3)
 
+@DOCSTRING(streamtube)
+
+@DOCSTRING(streamline)
+
+@DOCSTRING(stream2)
+
+@DOCSTRING(stream3)
+
 @DOCSTRING(compass)
 
 @DOCSTRING(feather)
diff -r f998e243fa78 -r fcae01084f8f libinterp/dldfcn/module-files
--- a/libinterp/dldfcn/module-files	Wed Nov 13 11:11:49 2019 +0100
+++ b/libinterp/dldfcn/module-files	Sat Nov 23 11:47:16 2019 +0100
@@ -18,5 +18,6 @@
 fftw.cc|$(FFTW_XCPPFLAGS)|$(FFTW_XLDFLAGS)|$(FFTW_XLIBS)
 gzip.cc|$(Z_CPPFLAGS) $(BZ2_CPPFLAGS)|$(Z_LDFLAGS) $(BZ2_LDFLAGS)|$(Z_LIBS) $(BZ2_LIBS)
 qr.cc|$(QRUPDATE_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(QRUPDATE_LDFLAGS) $(SPARSE_XLDFLAGS)|$(QRUPDATE_LIBS) $(SPARSE_XLIBS)
+stream_euler.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
 symbfact.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
 symrcm.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
diff -r f998e243fa78 -r fcae01084f8f libinterp/dldfcn/stream_euler.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/dldfcn/stream_euler.cc	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,549 @@
+/*
+
+Copyright (C) 2019 HG200
+
+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/>.
+
+*/
+
+/*
+
+ References:
+
+ @article{
+    title = {Particle Tracing Algorithms for 3D Curvilinear Grids},
+    year = {2000},
+    author = {Nielson, Gregory and Uller, H. and Sadarjoen, I. and Walsum, Theo and Hin, Andrea and Post, Frits}
+ }
+
+ @article{
+    title = {Sources of error in the graphical analysis of CFD results},
+    publisher = {Journal of Scientific Computing},
+    year = {1988},
+    volume = {3},
+    number = {2},
+    pages = {149--164},
+    author = {Buning, Pieter G.},
+ }
+
+*/
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include "defun-dld.h"
+#include "error.h"
+#include "ovl.h"
+
+// Coordinates of a point in C-Space (unit square mesh)
+
+typedef struct
+{
+  double x, y;
+} Vector2;
+
+// The integer- 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.
+
+typedef struct
+{
+  double fcx, fcy;
+  signed long idx, idy;
+} Cell2;
+
+typedef struct
+{
+  double x, y, z;
+} Vector3;
+
+typedef struct
+{
+  double fcx, fcy, fcz;
+  signed long idx, idy, idz;
+} Cell3;
+
+static inline void
+number_to_fractional (signed long *id, double *fc, const double u)
+{
+  *id = floor(u);
+  *fc = u - *id;
+}
+
+static inline octave_idx_type
+handle_border_index (const octave_idx_type id, const octave_idx_type N)
+{
+  return (id < N - 1 ? id : N - 2);
+}
+
+static inline void
+handle_border (octave_idx_type *id2, double *fc2, const octave_idx_type id1,
+               const double fc1, const octave_idx_type N)
+{
+  if (id1 < N - 1)
+    {
+      *id2 = id1;
+      *fc2 = fc1;
+    }
+  else
+    {
+      *id2 = N - 2;
+      *fc2 = 1.0;
+    }
+}
+
+// Is there an API function for the command "dst=src(1:nverts,:)" ?
+
+static Matrix
+crop_buffer (const Matrix src, const octave_idx_type nverts,
+             const octave_idx_type cols)
+{
+  Matrix dst (nverts, cols);
+
+  for (octave_idx_type j = 0; j < cols; j++)
+    for (octave_idx_type i = 0; i < nverts; i++)
+      {
+        dst(i, j) = src(i, j);
+      }
+
+  return (dst);
+}
+
+static inline double
+bilinear (const double u11, const double u21, const double u12,
+          const double u22, const double x, const double y)
+{
+  return (u11 * (1-x) * (1-y) +
+          u21 * x * (1-y) +
+          u12 * (1-x) * y +
+          u22 * x * y);
+}
+
+static inline Cell2
+vector_to_cell2d (const Vector2 X)
+{
+  Cell2 Z;
+
+  number_to_fractional (&Z.idx, &Z.fcx, X.x);
+  number_to_fractional (&Z.idy, &Z.fcy, X.y);
+
+  return (Z);
+}
+
+static inline bool
+is_in_definition_set2d (const Cell2 X, const octave_idx_type cols,
+                        const octave_idx_type rows)
+{
+ return ( (((X.idx >= 0) && (X.idx < cols-1)) ||
+           ((X.idx == cols-1) && (X.fcx == 0.0))) &&
+          (((X.idy >= 0) && (X.idy < rows-1)) ||
+           ((X.idy == rows-1) && (X.fcy == 0.0))) );
+}
+
+static inline Vector2
+add2d (const Cell2 X, const Vector2 Y)
+{
+  Vector2 Z = {X.idx + X.fcx + Y.x,
+               X.idy + X.fcy + Y.y};
+
+  return (Z);
+}
+
+static inline Vector2
+vector_interpolation2d (const Cell2 X, const Matrix& u, const Matrix& v,
+                        const octave_idx_type cols, const octave_idx_type rows)
+{
+  Vector2 V;
+  double fcx, fcy;
+  octave_idx_type idx, idy;
+
+  handle_border (&idx, &fcx, X.idx, X.fcx, cols);
+  handle_border (&idy, &fcy, X.idy, X.fcy, rows);
+
+  V.x = bilinear (u(idy, idx), u(idy, idx+1), u(idy+1, idx),
+                  u(idy+1, idx+1), fcx, fcy);
+  V.y = bilinear (v(idy, idx), v(idy, idx+1), v(idy+1, idx),
+                  v(idy+1, idx+1), fcx, fcy);
+
+  return (V);
+}
+
+// Apply the Jacobian matrix on the vector V.
+// The step vector length is set to h.
+
+static inline Vector2
+calculate_step_vector2d (const Cell2 X, const Vector2 V,
+                         const RowVector& tx, const RowVector& ty,
+                         const octave_idx_type cols, const octave_idx_type rows,
+                         const double h)
+{
+  Vector2 S;
+
+  const octave_idx_type idx = handle_border_index (X.idx, cols);
+  const octave_idx_type idy = handle_border_index (X.idy, rows);
+
+  const double x = V.x * tx(idx);
+  const double y = V.y * ty(idy);
+  const double n = 1.0 / sqrt (x*x + y*y);
+  S.x = h * n * x;
+  S.y = h * n * y;
+
+  return (S);
+}
+
+static inline bool
+is_singular2d (const Vector2 V)
+{
+  return ((octave::math::isnan (V.x) || octave::math::isnan (V.y)) ||
+          ((V.x == 0) && (V.y == 0)));
+}
+
+static void
+euler2d (const octave_idx_type cols, const octave_idx_type rows,
+         const Matrix& u, const Matrix& v,
+         const RowVector& tx, const RowVector& ty,
+         const double zeta, const double xi,
+         const double h, const octave_idx_type maxnverts,
+         Matrix& buffer, octave_idx_type *nverts)
+{
+  Vector2 V0, V1, S0, X1, Xnxt, S1;
+  const Vector2 X0 = {zeta, xi};
+  Cell2 X0f, X1f;
+
+  octave_idx_type i = 0;
+
+  buffer(i, 0) = X0.x;
+  buffer(i, 1) = X0.y;
+
+  X0f = vector_to_cell2d (X0);
+  while (1)
+    {
+      if (! is_in_definition_set2d (X0f, cols, rows))
+        break;
+
+      V0 = vector_interpolation2d (X0f, u, v, cols, rows);
+      if (is_singular2d (V0))
+        break;
+
+      S0 = calculate_step_vector2d (X0f, V0, tx, ty, cols, rows, h);
+
+      X1 = add2d (X0f, S0);
+      X1f = vector_to_cell2d (X1);
+      if (! is_in_definition_set2d (X1f, cols, rows))
+        break;
+
+      V1 = vector_interpolation2d (X1f, u, v, cols, rows);
+      if (is_singular2d (V1))
+        break;
+
+      S1 = calculate_step_vector2d (X1f, V1, tx, ty, cols, rows, h);
+
+      // Runge Kutta - Heun's Scheme
+      const Vector2 S = {0.5 * (S0.x + S1.x),
+                         0.5 * (S0.y + S1.y)};
+      Xnxt = add2d (X0f, S);
+
+      X0f = vector_to_cell2d (Xnxt);
+      if (! is_in_definition_set2d (X0f, cols, rows))
+        break;
+
+      i++;
+      buffer(i, 0) = Xnxt.x;
+      buffer(i, 1) = Xnxt.y;
+
+      if (i + 1 >= maxnverts)
+        break;
+    }
+
+  *nverts = i + 1;
+}
+
+static inline double
+trilinear (const double u111, const double u211, const double u121,
+           const double u221, const double u112, const double u212,
+           const double u122, const double u222,
+           const double x, const double y, const double z)
+{
+  return (u111 * (1-x) * (1-y) * (1-z) +
+          u211 * x * (1-y) * (1-z) +
+          u121 * (1-x) * y * (1-z) +
+          u221 * x * y * (1-z) +
+          u112 * (1-x) * (1-y) * z +
+          u212 * x * (1-y) * z +
+          u122 * (1-x) * y * z +
+          u222 * x * y * z);
+}
+
+static inline Cell3
+vector_to_cell3d (const Vector3 X)
+{
+  Cell3 Z;
+
+  number_to_fractional (&Z.idx, &Z.fcx, X.x);
+  number_to_fractional (&Z.idy, &Z.fcy, X.y);
+  number_to_fractional (&Z.idz, &Z.fcz, X.z);
+
+  return (Z);
+}
+
+static inline bool
+is_in_definition_set3d (const Cell3 X, const octave_idx_type nx,
+                        const octave_idx_type ny, const octave_idx_type nz)
+{
+  return ( (((X.idx >= 0) && (X.idx < nx-1)) ||
+            ((X.idx == nx-1) && (X.fcx == 0.0))) &&
+           (((X.idy >= 0) && (X.idy < ny-1)) ||
+            ((X.idy == ny-1) && (X.fcy == 0.0))) &&
+           (((X.idz >= 0) && (X.idz < nz-1)) ||
+            ((X.idz == nz-1) && (X.fcz == 0.0))) );
+}
+
+static inline Vector3
+add3d (const Cell3 X, const Vector3 Y)
+{
+  Vector3 Z = {X.idx + X.fcx + Y.x,
+               X.idy + X.fcy + Y.y,
+               X.idz + X.fcz + Y.z};
+
+  return (Z);
+}
+
+static inline Vector3
+vector_interpolation3d (const Cell3 X, const NDArray& u, const NDArray& v,
+                        const NDArray& w, const octave_idx_type nx,
+                        const octave_idx_type ny, const octave_idx_type nz)
+{
+  Vector3 V;
+  double fcx, fcy, fcz;
+  octave_idx_type idx, idy, idz;
+
+  handle_border (&idx, &fcx, X.idx, X.fcx, nx);
+  handle_border (&idy, &fcy, X.idy, X.fcy, ny);
+  handle_border (&idz, &fcz, X.idz, X.fcz, nz);
+
+  V.x = trilinear (u(idy, idx, idz), u(idy, idx+1, idz),
+                   u(idy+1, idx, idz), u(idy+1, idx+1, idz),
+                   u(idy, idx, idz+1), u(idy, idx+1, idz+1),
+                   u(idy+1, idx, idz+1), u(idy+1, idx+1, idz+1),
+                   fcx, fcy, fcz);
+  V.y = trilinear (v(idy, idx, idz), v(idy, idx+1, idz),
+                   v(idy+1, idx, idz), v(idy+1, idx+1, idz),
+                   v(idy, idx, idz+1), v(idy, idx+1, idz+1),
+                   v(idy+1, idx, idz+1), v(idy+1, idx+1, idz+1),
+                   fcx, fcy, fcz);
+  V.z = trilinear (w(idy, idx, idz), w(idy, idx+1, idz),
+                   w(idy+1, idx, idz), w(idy+1, idx+1, idz),
+                   w(idy, idx, idz+1), w(idy, idx+1, idz+1),
+                   w(idy+1, idx, idz+1), w(idy+1, idx+1, idz+1),
+                   fcx, fcy, fcz);
+
+  return (V);
+}
+
+static inline Vector3
+calculate_step_vector3d (const Cell3 X, const Vector3 V,
+                         const RowVector& tx, const RowVector& ty, const RowVector& tz,
+                         const octave_idx_type nx, const octave_idx_type ny,
+                         const octave_idx_type nz, const double h)
+{
+  Vector3 S;
+
+  const octave_idx_type idx = handle_border_index (X.idx, nx);
+  const octave_idx_type idy = handle_border_index (X.idy, ny);
+  const octave_idx_type idz = handle_border_index (X.idz, nz);
+
+  const double x = V.x * tx(idx);
+  const double y = V.y * ty(idy);
+  const double z = V.z * tz(idz);
+  const double n = 1.0 / sqrt (x*x + y*y + z*z);
+  S.x = h * n * x;
+  S.y = h * n * y;
+  S.z = h * n * z;
+
+  return (S);
+}
+
+static inline bool
+is_singular3d (const Vector3 V)
+{
+  return ((octave::math::isnan (V.x) || octave::math::isnan (V.y) ||
+           octave::math::isnan (V.z)) ||
+          ((V.x == 0) && (V.y == 0) && (V.z == 0)));
+}
+
+static void
+euler3d (const octave_idx_type nx, const octave_idx_type ny, const octave_idx_type nz,
+         const NDArray& u, const NDArray& v, const NDArray& w,
+         const RowVector& tx, const RowVector& ty, const RowVector& tz,
+         const double zeta, const double xi, const double rho,
+         const double h, const octave_idx_type maxnverts,
+         Matrix& buffer, octave_idx_type *nverts)
+{
+  Vector3 V0, V1, S0, X1, Xnxt, S1;
+  const Vector3 X0 = {zeta, xi, rho};
+  Cell3 X0f, X1f;
+
+  octave_idx_type i = 0;
+  buffer(i, 0) = X0.x;
+  buffer(i, 1) = X0.y;
+  buffer(i, 2) = X0.z;
+
+  X0f = vector_to_cell3d (X0);
+  while (1)
+    {
+      if (! is_in_definition_set3d (X0f, nx, ny, nz))
+        break;
+
+      V0 = vector_interpolation3d (X0f, u, v, w, nx, ny, nz);
+      if (is_singular3d (V0))
+        break;
+
+      S0 = calculate_step_vector3d (X0f, V0, tx, ty, tz, nx, ny, nz, h);
+
+      X1 = add3d (X0f, S0);
+
+      X1f = vector_to_cell3d (X1);
+      if (! is_in_definition_set3d (X1f, nx, ny, nz))
+        break;
+
+      V1 = vector_interpolation3d (X1f, u, v, w, nx, ny, nz);
+      if (is_singular3d (V1))
+        break;
+
+      S1 = calculate_step_vector3d (X1f, V1, tx, ty, tz, nx, ny, nz, h);
+
+      // Runge Kutta - Heun's Scheme
+      const Vector3 S = {0.5 * (S0.x + S1.x),
+                         0.5 * (S0.y + S1.y),
+                         0.5 * (S0.z + S1.z)};
+      Xnxt = add3d (X0f, S);
+
+      X0f = vector_to_cell3d (Xnxt);
+      if (! is_in_definition_set3d (X0f, nx, ny, nz))
+        break;
+
+      i++;
+      buffer(i, 0) = Xnxt.x;
+      buffer(i, 1) = Xnxt.y;
+      buffer(i, 2) = Xnxt.z;
+
+      if (i + 1 >= maxnverts)
+        break;
+
+    }
+
+  *nverts = i + 1;
+}
+
+static octave_value
+streameuler2d_internal (const octave_value_list& args)
+{
+
+  const int nargin = args.length ();
+  if (nargin != 8)
+    print_usage ();
+
+  const Matrix U = args(0).matrix_value ();
+  const Matrix V = args(1).matrix_value ();
+  const RowVector TX = args(2).row_vector_value ();
+  const RowVector TY = args(3).row_vector_value ();
+  const double zeta = args(4).double_value ();
+  const double xi = args(5).double_value ();
+  const double h = args(6).double_value ();
+  const octave_idx_type maxnverts = args(7).idx_type_value ();
+
+  const octave_idx_type rows = U.rows ();
+  const octave_idx_type cols = U.columns ();
+
+  octave_idx_type nverts;
+  Matrix buffer (maxnverts, 2);
+
+  euler2d (cols, rows, U, V, TX, TY, zeta, xi, h, maxnverts,
+           buffer, &nverts);
+
+  Matrix xy = crop_buffer (buffer, nverts, 2);
+
+  return  octave_value (xy);
+}
+
+static octave_value
+streameuler3d_internal (const octave_value_list& args, const char *fcn)
+{
+
+  const int nargin = args.length ();
+  if (nargin != 11)
+    print_usage ();
+
+  const NDArray U = args(0).array_value ();
+  const NDArray V = args(1).array_value ();
+  const NDArray W = args(2).array_value ();
+  const RowVector TX = args(3).row_vector_value ();
+  const RowVector TY = args(4).row_vector_value ();
+  const RowVector TZ = args(5).row_vector_value ();
+  const double zeta = args(6).double_value ();
+  const double xi = args(7).double_value ();
+  const double rho = args(8).double_value ();
+  const double h = args(9).double_value ();
+  const octave_idx_type maxnverts = args(10).idx_type_value ();
+
+  const dim_vector dims = args(0).dims ();
+  const int ndims = dims.ndims ();
+  if (ndims != 3)
+    {
+      error ("%s: dimension must be 3", fcn);
+    }
+
+  octave_idx_type nverts;
+  Matrix buffer (maxnverts, 3);
+
+  euler3d (dims(1), dims(0), dims(2), U, V, W, TX, TY, TZ, zeta, xi, rho,
+           h, maxnverts, buffer, &nverts);
+
+  Matrix xyz = crop_buffer (buffer, nverts, 3);
+
+  return  octave_value (xyz);
+}
+
+DEFUN_DLD (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}].
+
+@seealso{streamline, stream2, stream3, streameuler3d}
+@end deftypefn */)
+{
+  return streameuler2d_internal (args);
+}
+
+DEFUN_DLD (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}].
+
+@seealso{streamline, stream2, stream3, streameuler2d}
+@end deftypefn */)
+{
+  return streameuler3d_internal (args, "streameuler3d");
+}
+
diff -r f998e243fa78 -r fcae01084f8f scripts/plot/draw/module.mk
--- a/scripts/plot/draw/module.mk	Wed Nov 13 11:11:49 2019 +0100
+++ b/scripts/plot/draw/module.mk	Sat Nov 23 11:47:16 2019 +0100
@@ -93,6 +93,10 @@
   %reldir%/stem.m \
   %reldir%/stem3.m \
   %reldir%/stemleaf.m \
+  %reldir%/stream2.m \
+  %reldir%/stream3.m \
+  %reldir%/streamline.m \
+  %reldir%/streamtube.m \
   %reldir%/surf.m \
   %reldir%/surface.m \
   %reldir%/surfc.m \
diff -r f998e243fa78 -r fcae01084f8f scripts/plot/draw/stream2.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/stream2.m	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,195 @@
+## Copyright (C) 2019 HG200
+##
+## 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  {} {@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}")
+## Compute 2D streamline data.
+##
+## Calculates streamlines of a vector field given by [@var{u}, @var{v}].
+## The vector field is defined over a rectangular grid given by
+## [@var{x}, @var{y}]. The streamlines start at the seed points
+## [@var{sx}, @var{sy}]. The returned value @var{xy} contains a cell array
+## of vertex arrays. If the starting point is outside the vector field,
+## [] is returned.
+##
+## The input parameter @var{options} is a 2D vector of the form
+## [@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.
+## The second parameter specifies the maximum number of segments used to
+## create a streamline (default 10000).
+##
+## The return value @var{xy} is a @nospell{nverts x 2} matrix containing the
+## coordinates of the field line segments.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y] = meshgrid (0:3);
+## u = 2 * x;
+## v = y;
+## xy = stream2 (x, y, u, v, 1.0, 0.5);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream3}
+##
+## @end deftypefn
+
+## References:
+##
+## @article{
+##    title = {Particle Tracing Algorithms for 3D Curvilinear Grids},
+##    year = {2000},
+##    author = {Nielson, Gregory and Uller, H. and Sadarjoen, I. and Walsum, Theo and Hin, Andrea and Post, Frits}
+## }
+##
+## @article{
+##    title = {Sources of error in the graphical analysis of CFD results},
+##    publisher = {Journal of Scientific Computing},
+##    year = {1988},
+##    volume = {3},
+##    number = {2},
+##    pages = {149--164},
+##    author = {Buning, Pieter G.},
+## }
+
+function xy = stream2 (varargin)
+
+  options = [];
+  switch (length (varargin))
+    case (0)
+      print_usage ();
+    case {4,5}
+      if (length (varargin) == 4)
+        [u, v, spx, spy] = varargin{:};
+      else
+        [u, v, spx, spy, options] = varargin{:};
+      end
+      [m, n] = size (u);
+      [x, y] = meshgrid (1:n, 1:m);
+    case (6)
+      [x, y, u, v, spx, spy] = varargin{:};
+    case (7)
+      [x, y, u, v, spx, spy, options] = varargin{:};
+    otherwise
+      error ('stream2: unknown input parameter count');
+  end
+
+  h = 0.1;
+  maxnverts = 10000;
+  if (~ isempty (options))
+    switch (length (options))
+      case (1)
+        h = options(1);
+      case (2)
+        h = options(1);
+        maxnverts = options(2);
+      otherwise
+        error ('stream2: wrong options length');
+    end
+  end
+
+  if (~ isnumeric (h) || (h == 0))
+    error ('stream2: step size error');
+  end
+  if (~ isnumeric (maxnverts) || (maxnverts < 1))
+    error ('stream2: max num vertices error');
+  end
+  if (~ (isequal (size (u), size (v), size (x), size (y)) && ...
+         isequal (size (spx), size (spy))) )
+    error ('stream2: matrix dimensions must match');
+  end
+
+  gx = x(1, :);
+  gy = y(:, 1)';
+
+  ## Jacobian Matrix
+  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 ('sream2: ill shaped elements in mesh');
+  end
+  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;
+  dx(end + 1) = 0;
+  dy(end + 1) = 0;
+
+  px = spx(:);
+  py = spy(:);
+
+  for nseed = 1:length (px)
+
+    xp = px(nseed);
+    yp = py(nseed);
+    idx = find (logical (diff (gx <= xp)), 1);
+    if (gx(end) == xp)
+      idx = numel(gx);
+    end
+    idy = find (logical (diff (gy <= yp)), 1);
+    if (gy(end) == yp)
+      idy = numel(gy);
+    end
+
+    if (isempty (idx) || isempty (idy))
+      xy{nseed} = [];
+    else
+      ## Transform seed from P coordinates to C coordinates
+      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);
+
+      ## Transform from C coordinates to P coordinates
+      idu = floor (C(:, 1));
+      idv = floor (C(:, 2));
+      xy{nseed} = [gx(idu + 1)' + (C(:, 1) - idu).*(dx(idu + 1)'), ...
+                   gy(idv + 1)' + (C(:, 2) - idv).*(dy(idv + 1)')];
+    end
+
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y] = meshgrid (-5:5, -4:4);
+%! u = x - 2 * y;
+%! v = 2 * x - 3 * y;
+%! sx = [3, 0, -1, -2, -3, 0, 1, 2];
+%! sy = [3, 3, 3, 3, -3, -3, -3, -3];
+%! h = streamline (x, y, u, v, sx, sy, 0.05);
+%! set (h, 'color', 'r');
+%! hold on;
+%! quiver (x, y, u, v);
+%! scatter (sx(:), sy(:), 20, 'filled', 'o', 'markerfacecolor', 'r');
+%! grid on;
+%! title ('Asymptotically Stable Equilibrium');
+%! 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]);
+%! assert(numel(xy{:}),10);
diff -r f998e243fa78 -r fcae01084f8f scripts/plot/draw/stream3.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/stream3.m	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,218 @@
+## Copyright (C) 2019 HG200
+##
+## 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  {} {@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}")
+## Compute 3D streamline data.
+##
+## Calculates streamlines of a vector field given by [@var{u}, @var{v}, @var{w}].
+## The vector field is defined over a rectangular grid given by
+## [@var{x}, @var{y}, @var{z}]. The streamlines start at the seed points
+## [@var{sx}, @var{sy}, @var{sz}]. The returned value @var{xyz}
+## contains a cell array of vertex arrays. If the starting point is outside
+## the vector field, [] is returned.
+##
+## The input parameter @var{options} is a 2D vector of the form
+## [@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.
+## The second parameter specifies the maximum number of segments used to
+## create a streamline (default 10000).
+##
+## The return value @var{xyz} is a @nospell{nverts x 3} matrix containing the
+## coordinates of the field line segments.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y, z] = meshgrid (0:3);
+## u = 2 * x;
+## v = y;
+## w = 3 * z;
+## xyz = stream3 (x, y, z, u, v, w, 1.0, 0.5, 0.0);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream2}
+##
+## @end deftypefn
+
+## References:
+##
+## @article{
+##    title = {Particle Tracing Algorithms for 3D Curvilinear Grids},
+##    year = {2000},
+##    author = {Nielson, Gregory and Uller, H. and Sadarjoen, I. and Walsum, Theo and Hin, Andrea and Post, Frits}
+## }
+##
+## @article{
+##    title = {Sources of error in the graphical analysis of CFD results},
+##    publisher = {Journal of Scientific Computing},
+##    year = {1988},
+##    volume = {3},
+##    number = {2},
+##    pages = {149--164},
+##    author = {Buning, Pieter G.},
+## }
+
+function xyz = stream3 (varargin)
+
+  options = [];
+  switch (length (varargin))
+    case (0)
+      print_usage ();
+    case {6,7}
+      if (length (varargin) == 6)
+        [u, v, w, spx, spy, spz] = varargin{:};
+      else
+        [u, v, w, spx, spy, spz, options] = varargin{:};
+      end
+      [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)
+      [x, y, z, u, v, w, spx, spy, spz, options] = varargin{:};
+    otherwise
+      error ('stream3: unknown input parameter count');
+  end
+
+  h = 0.1;
+  maxnverts = 10000;
+  if (~ isempty (options))
+    switch (length (options))
+      case (1)
+        h = options(1);
+      case (2)
+        h = options(1);
+        maxnverts = options(2);
+      otherwise
+        error ('stream3: wrong options length');
+    end
+  end
+
+  if (~ isnumeric (h) || (h == 0))
+    error ('stream3: step size error');
+  end
+  if (~ isnumeric (maxnverts) || (maxnverts < 1))
+    error ('stream3: max num vertices error');
+  end
+  if (~ (isequal (size (u), size (v), size (w), size (x), size (y), size (z)) && ...
+         isequal (size (spx), size (spy))) )
+    error ('stream3: matrix dimensions must match');
+  end
+
+  gx = x(1, :, 1);
+  gy = y(:, 1, 1)';
+  tmp = z(1, 1, :);
+  gz = tmp(:)';
+
+  ## Jacobian Matrix
+  dx = diff (gx);
+  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 ('sream3: ill shaped elements in mesh');
+  end
+  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;
+  tz(end + 1) = 0;
+  dx(end + 1) = 0;
+  dy(end + 1) = 0;
+  dz(end + 1) = 0;
+
+  px = spx(:);
+  py = spy(:);
+  pz = spz(:);
+
+  for nseed = 1:length (px)
+
+    xp = px(nseed);
+    yp = py(nseed);
+    zp = pz(nseed);
+    idx = find (logical (diff (gx <= xp)), 1);
+    if (gx(end) == xp)
+      idx = numel(gx);
+    end
+    idy = find (logical (diff (gy <= yp)), 1);
+    if (gy(end) == yp)
+      idy = numel(gy);
+    end
+    idz = find (logical (diff (gz <= zp)), 1);
+    if (gz(end) == zp)
+      idz = numel(gz);
+    end
+
+    if (isempty (idx) || isempty (idy) || isempty (idz))
+      xyz{nseed} = [];
+    else
+      ## Transform seed from P coordinates to C coordinates
+      zeta = (idx - 1) + (xp - gx(idx)) * tx(idx);
+      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);
+
+      ## Transform from C coordinates to P coordinates
+      idu = floor (C(:, 1));
+      idv = floor (C(:, 2));
+      idw = floor (C(:, 3));
+      xyz{nseed} = [gx(idu + 1)' + (C(:, 1) - idu).*(dx(idu + 1)'), ...
+                    gy(idv + 1)' + (C(:, 2) - idv).*(dy(idv + 1)'), ...
+                    gz(idw + 1)' + (C(:, 3) - idw).*(dz(idw + 1)')];
+    end
+
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-30:1:30, -30:1:30, 0:1:50);
+%! s = 10;
+%! b = 8 / 3;
+%! r = 28;
+%! u = s * (y - x);
+%! v = r * x - y - x.*z;
+%! w = x.*y - b * z;
+%! hold on;
+%! sx = 0.1;
+%! sy = 0.1;
+%! sz = 0.1;
+%! plot3 (sx, sy, sz, '.r', 'markersize', 15);
+%! h = streamline (x, y, z, u, v, w, sx, sy, sz, [0.1, 50000]);
+%! set (h, 'color', 'r');
+%! view (3);
+%! title ('Lorenz System');
+%! grid on;
+%! axis equal;
+
+%!test
+%! [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);
diff -r f998e243fa78 -r fcae01084f8f scripts/plot/draw/streamline.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/streamline.m	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,161 @@
+## Copyright (C) 2019 HG200
+##
+## 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  {} {} 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 (@var{hax}, @dots{})
+## @deftypefnx {} {@var{h} =} streamline (@dots{})
+## Plot streamlines of 2D or 3D vector fields.
+##
+## Plot streamlines of a 2D or 3D vector field given by
+## [@var{u}, @var{v}] or [@var{u}, @var{v}, @var{w}]. The vector field
+## is defined over a rectangular grid given by [@var{x}, @var{y}]
+## or [@var{x}, @var{y}, @var{z}]. The streamlines start at the seed
+## points [@var{sx}, @var{sy}] or [@var{sx}, @var{sy}, @var{sz}].
+##
+## The input parameter @var{options} is a 2D vector of the form
+## [@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.
+## The second parameter specifies the maximum number of segments used to
+## create a streamline (default 10000).
+##
+## 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 hggroup
+## comprising the field lines.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y] = meshgrid (-1.5:0.2:2, -1:0.2:2);
+## u = - x / 4 - y;
+## v = x - y / 4;
+## streamline (x, y, u, v, 1.7, 1.5);
+## @end group
+## @end example
+##
+## @seealso{stream2, stream3}
+##
+## @end deftypefn
+
+function h = streamline (varargin)
+
+  if (nargin == 0)
+    print_usage ();
+  end
+
+  [hax, varargin] = __plt_get_axis_arg__ ("streamline", varargin{:});
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  end
+
+  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];
+        end
+      end
+    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];
+        end
+      end
+    otherwise
+      error ('streamline: 2D or 3D input data only');
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y] = meshgrid (-2:0.5:2);
+%! u = - y - x / 2;
+%! v = x - y / 2;
+%! [sx, sy] = meshgrid (-2:2:2);
+%! h = streamline (x, y, u, v, sx, sy);
+%! set (h, 'color', 'r');
+%! hold on;
+%! quiver (x, y, u, v);
+%! scatter (sx(:), sy(:), 20, 'filled', 'o', 'markerfacecolor', 'r');
+%! title ('Spiral Sink');
+%! grid on;
+%! axis equal;
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-3:3);
+%! u = - x / 2 - y;
+%! v = x - y / 2;
+%! w = - z;
+%! [sx, sy, sz] = meshgrid (3, 0:1.5:1.5, 0:1.5:3);
+%! h = streamline (x, y, z, u, v, w, sx, sy, sz);
+%! set (h, 'color', 'r');
+%! hold on;
+%! quiver3 (x, y, z, u, v, w);
+%! scatter3 (sx(:), sy(:), sz(:), 20, 'filled', 'o', 'markerfacecolor', 'r');
+%! view (3);
+%! title ('Spiral Sink');
+%! grid on;
+%! axis equal;
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-1:0.4:1, -1:0.4:1, -3:0.3:0);
+%! a = 0.08;
+%! b = 0.04;
+%! u = - a * x - y;
+%! v = x - a * y;
+%! w = - b * ones (size (x));
+%! hold on;
+%! sx = 1.0;
+%! sy = 0.0;
+%! sz = 0.0;
+%! plot3 (sx, sy, sz, '.r', 'markersize', 15);
+%! t = linspace (0, 12 * 2 * pi(), 500);
+%! tx = exp (-a * t).*cos (t);
+%! ty = exp (-a * t).*sin (t);
+%! tz = - b * t;
+%! plot3 (tx, ty, tz, '-b');
+%! h = streamline (x, y, z, u, v, w, sx, sy, sz);
+%! set (h, 'color', 'r');
+%! view (3);
+%! title ('Heuns Scheme (red) vs. Analytical Solution (blue)');
+%! grid on;
+%! axis equal tight;
+
+%!test
+%!error streamline ()
diff -r f998e243fa78 -r fcae01084f8f scripts/plot/draw/streamtube.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/streamtube.m	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,375 @@
+## Copyright (C) 2019 HG200
+##
+## 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  {} {} 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 (@var{hax}, @dots{})
+## @deftypefnx {} {@var{h} =} streamtube (@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 [@var{u}, @var{v}, @var{w}] and is defined over a
+## rectangular grid given by [@var{x}, @var{y}, @var{z}]. The streamtubes start
+## at the seed points [@var{sx}, @var{sy}, @var{sz}].
+##
+## The input parameter @var{options} is a 2D vector of the form
+## [@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).
+##
+## Streamtube can be called with a cell array containing precomputed streamline
+## data. To do this, @var{vertices} must be created with the 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 patch 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;
+## streamtube (x, y, z, u, v, w, 1, 0, 0);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream3}
+##
+## @end deftypefn
+
+## References:
+##
+## @inproceedings{
+##    title = {Visualization of 3-D vector fields - Variations on a stream},
+##    author = {Dave Darmofal and Robert Haimes},
+##    year = {1992}
+## }
+
+function h = streamtube (varargin)
+
+  if (nargin == 0)
+    print_usage ();
+  end
+
+  [hax, varargin] = __plt_get_axis_arg__ ("streamtube", varargin{:});
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  end
+
+  options = [];
+  xyz = [];
+  switch (length (varargin))
+    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);
+      end
+    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 ('streamtube: unknown input parameter count');
+  end
+
+  scale = 1;
+  n = 20;
+  if (~ isempty (options))
+    switch (length (options))
+      case (1)
+        scale = options(1);
+      case (2)
+        scale = options(1);
+        n = options(2);
+      otherwise
+        error ('streamtube: wrong options length');
+    end
+  end
+
+  if (~ isnumeric (scale) || (scale <= 0))
+    error ('streamtube: scale error');
+  end
+  if (~ isnumeric (n) || (n < 3))
+    error ('streamtube: number of polygons for tube circumference too small');
+  end
+  if isempty (xyz)
+    xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.5);
+  end
+
+  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)));
+
+  ## 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);
+
+  h = [];
+  for i = 1:length (xyz)
+
+    sl = xyz{i};
+    [nverts, ~] = size(sl);
+    if (~ isempty (sl)) && (nverts > 2)
+
+      divsl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
+      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);
+
+      tmp = plottube (hax, sl, divsl, vv, vmax, vmin, r0, n);
+      h = [h, tmp];
+
+    end
+
+  end
+
+end
+
+function h = plottube (hax, sl, divsl, vv, vmax, vmin, r0, npoly)
+
+  issingular = find (isnan (divsl), 1, 'first');
+  if ~ isempty (issingular)
+     maxnverts = issingular - 1;
+  else
+     maxnverts = length (sl);
+  end
+  if maxnverts < 3
+    error ('streamtube: too less data to show');
+  end
+
+  if (vmax == vmin)
+    colscale = 0.0;
+  else
+    colscale = 1.0 / (vmax - vmin);
+  end
+
+  phi = linspace (0, 2*pi (), npoly);
+
+  X0 = sl(1, :);
+  X1 = sl(2, :);
+
+  ## 1st rotation axis
+  R = X1 - X0;
+
+  ## Initial radius
+  vold = vv(1);
+  vact = vv(2);
+  ract = r0 * exp (0.5 * divsl(2) * norm (R) / vact) * sqrt (vold / vact);
+  vold = vact;
+  rold = ract;
+
+  ## Guide point and its rotation to create a segment
+  N = get_guide_point (X0, X1);
+  K = ract * N;
+  XS = rotation (R, K, phi) + repmat (X1', 1, npoly);
+
+  px = zeros (4, npoly * (maxnverts - 2));
+  py = zeros (4, npoly * (maxnverts - 2));
+  pz = zeros (4, npoly * (maxnverts - 2));
+  pc = zeros (4, npoly * (maxnverts - 2));
+
+  for j = 3:maxnverts
+
+    KK = K;
+    X0 = X1;
+    X1 = sl(j, :);
+    R = X1 - X0;
+
+    ## Tube radius
+    vact = vv(j);
+    ract = rold * exp (0.5 * divsl(j) * norm (R) / vact) * sqrt (vold / vact);
+    vold = vact;
+    rold = ract;
+
+    ## Project K onto R and get the difference in order to calculate the next
+    ## guiding point
+    Kp = KK - R * dot (KK, R) / (norm (R)^2);
+    K = ract * Kp / norm (Kp);
+
+    XSold = XS;
+    ## Rotate the guiding point around R and collect patch vertices
+    XS = rotation (R, K, phi) + repmat (X1', 1, npoly);
+    [tx, ty, tz] = segment_patch_data (XS, XSold);
+
+    from = (j - 3) * npoly + 1;
+    to = (j + 1 - 3) * npoly;
+    px(:, from:to) = tx;
+    py(:, from:to) = ty;
+    pz(:, from:to) = tz;
+    pmc = colscale * (vact - vmin) * ones (4, npoly);
+    pc(:, from:to) = pmc;
+
+  end
+
+  h = patch (hax, px, py, pz, pc);
+
+end
+
+## Find N orthogonal to (X1 - X0)
+function N = get_guide_point (X0, X1)
+
+  S = X1 - X0;
+
+  if ((S(3) == 0) && (S(1) == -S(2)))
+    N = [- S(2) - S(3), S(1), S(1)];
+  else
+    N = [S(3), S(3), - S(1) - S(2)];
+  end
+
+  N = N / norm(N);
+
+end
+
+## Create patch data to draw a segment
+## from starting point XS to ending point XE
+function [px, py, pz] = segment_patch_data (XS, XE)
+
+  [~, npoly] = size (XS);
+
+  px = zeros (4, npoly);
+  py = zeros (4, npoly);
+  pz = zeros (4, npoly);
+
+  px(1, :) = XS(1, :);
+  px(2, :) = XE(1, :);
+  px(3, :) = [XE(1, 2:end), XE(1, 1)];
+  px(4, :) = [XS(1, 2:end), XS(1, 1)];
+
+  py(1, :) = XS(2, :);
+  py(2, :) = XE(2, :);
+  py(3, :) = [XE(2, 2:end), XE(2, 1)];
+  py(4, :) = [XS(2, 2:end), XS(2, 1)];
+
+  pz(1, :) = XS(3, :);
+  pz(2, :) = XE(3, :);
+  pz(3, :) = [XE(3, 2:end), XE(3, 1)];
+  pz(4, :) = [XS(3, 2:end), XS(3, 1)];
+
+end
+
+## A: Axis of rotation
+## X: Guiding point
+## phi: Angles
+## Y: Rotated points
+function Y = rotation (A, X, phi)
+
+  U = A / norm (A);
+  cp = cos (phi);
+  sp = sin (phi);
+
+  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));
+
+end
+
+%!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;
+%! hax = axes ();
+%! streamtube (hax, 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;
+%! props = {'cameraviewangle', 'dataaspectratio', 'plotboxaspectratio'};
+%! set(hax, props, get (hax, props));
+%! 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);
+%! hax = axes ();
+%! xyz = stream3 (x, y, z, u, v, w, sx, sy, sz, [0.1, 60]);
+%! streamtube (hax, 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;
+%! props = {'cameraviewangle', 'dataaspectratio', 'plotboxaspectratio'};
+%! set(hax, props, get (hax, props));
+%! title('Integration Towards Sink');
+
+%!test
+%!error streamtube ()
+
