Mercurial > octave

changeset 27759:992e702ef0d7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add stream* functions (patch #9859). * libinterp/corefcn/stream-euler.cc: New file. * libinterp/corefcn/module.mk: Add new file to list of source files. * stream2.m, stream3.m, streamline.m, streamtube.m: New functions. * plot.txi: Add documentation of new functions to manual. * __unimplemented__.m: Remove implemented functions from list. * NEWS: Announce new functions.
author Markus Meisinger <chloros2@gmx.de>
date Sat, 23 Nov 2019 11:47:16 +0100
parents fc8668e1cb26
children 16e83787f970
files NEWS doc/interpreter/plot.txi libinterp/corefcn/module.mk libinterp/corefcn/stream-euler.cc scripts/help/__unimplemented__.m scripts/plot/draw/module.mk scripts/plot/draw/stream2.m scripts/plot/draw/stream3.m scripts/plot/draw/streamline.m scripts/plot/draw/streamtube.m
diffstat 10 files changed, 1500 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS	Sun Nov 17 10:26:50 2019 +0100
+++ b/NEWS	Sat Nov 23 11:47:16 2019 +0100
@@ -140,6 +140,10 @@
 * `rotx`
 * `roty`
 * `rotz`
+* `stream2`
+* `stream3`
+* `streamline`
+* `streamtube`
 * `uisetfont`
 * `verLessThan`
 * `web`
--- a/doc/interpreter/plot.txi	Sun Nov 17 10:26:50 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)
--- a/libinterp/corefcn/module.mk	Sun Nov 17 10:26:50 2019 +0100
+++ b/libinterp/corefcn/module.mk	Sat Nov 23 11:47:16 2019 +0100
@@ -238,6 +238,7 @@
   %reldir%/spparms.cc \
   %reldir%/sqrtm.cc \
   %reldir%/stack-frame.cc \
+  %reldir%/stream-euler.cc \
   %reldir%/strfind.cc \
   %reldir%/strfns.cc \
   %reldir%/sub2ind.cc \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/stream-euler.cc	Sat Nov 23 11:47:16 2019 +0100
@@ -0,0 +1,530 @@
+/*
+
+Copyright (C) 2019 Markus Meisinger
+
+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.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;
+    }
+}
+
+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 (true)
+    {
+      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 (true)
+    {
+      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 = buffer.extract (0, 0, nverts-1, 1);
+
+  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 = buffer.extract (0, 0, nverts-1, 2);
+
+  return octave_value (xyz);
+}
+
+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}].
+
+@seealso{streamline, stream2, stream3, streameuler3d}
+@end deftypefn */)
+{
+  return streameuler2d_internal (args);
+}
+
+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}].
+
+@seealso{streamline, stream2, stream3, streameuler2d}
+@end deftypefn */)
+{
+  return streameuler3d_internal (args, "streameuler3d");
+}
+
--- a/scripts/help/__unimplemented__.m	Sun Nov 17 10:26:50 2019 +0100
+++ b/scripts/help/__unimplemented__.m	Sat Nov 23 11:47:16 2019 +0100
@@ -1207,13 +1207,9 @@
   "step",
   "stopasync",
   "str2mat",
-  "stream2",
-  "stream3",
-  "streamline",
   "streamparticles",
   "streamribbon",
   "streamslice",
-  "streamtube",
   "string",
   "strings",
   "strip",
--- a/scripts/plot/draw/module.mk	Sun Nov 17 10:26:50 2019 +0100
+++ b/scripts/plot/draw/module.mk	Sat Nov 23 11:47:16 2019 +0100
@@ -94,6 +94,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 \
--- /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,199 @@
+## Copyright (C) 2019 Markus Meisinger
+##
+## 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{:};
+      endif
+      [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");
+  endswitch
+
+  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");
+    endswitch
+  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))) )
+    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");
+  endif
+
+  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 ("stream2: ill shaped elements in mesh");
+  endif
+  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);
+    endif
+    idy = find (logical (diff (gy <= yp)), 1);
+    if (gy(end) == yp)
+      idy = numel (gy);
+    endif
+
+    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).')];
+    endif
+
+  endfor
+
+endfunction
+
+%!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);
--- /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,223 @@
+## Copyright (C) 2019 Markus Meisinger
+##
+## 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{:};
+      endif
+      [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");
+  endswitch
+
+  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");
+    endswitch
+  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))) )
+    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");
+  endif
+
+  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 ("stream3: ill shaped elements in mesh");
+  endif
+  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);
+    endif
+    idy = find (logical (diff (gy <= yp)), 1);
+    if (gy(end) == yp)
+      idy = numel (gy);
+    endif
+    idz = find (logical (diff (gz <= zp)), 1);
+    if (gz(end) == zp)
+      idz = numel (gz);
+    endif
+
+    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).')];
+    endif
+
+  endfor
+
+endfunction
+
+%!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);
--- /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,160 @@
+## Copyright (C) 2019 Markus Meisinger
+##
+## 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 ();
+  endif
+
+  [hax, varargin] = __plt_get_axis_arg__ ("streamline", varargin{:});
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  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
+
+endfunction
+
+%!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;
+
+%!error streamline ()
--- /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,371 @@
+## Copyright (C) 2019 Markus Meisinger
+##
+## 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 tubes are colored depending on the local vector field strength.
+##
+## 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 ();
+  endif
+
+  [hax, varargin] = __plt_get_axis_arg__ ("streamtube", varargin{:});
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  endif
+
+  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);
+      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 ("streamtube: unknown input parameter count");
+  endswitch
+
+  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");
+    endswitch
+  endif
+
+  if ((! isnumeric (scale)) || (scale <= 0))
+    error ("streamtube: scale error");
+  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)));
+
+  ## 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];
+
+    endif
+
+  endfor
+
+endfunction
+
+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);
+  endif
+  if (maxnverts < 3)
+    error ("streamtube: too less data to show");
+  endif
+
+  if (vmax == vmin)
+    colscale = 0.0;
+  else
+    colscale = 1.0 / (vmax - vmin);
+  endif
+
+  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;
+    pc(:, from:to) = colscale * (vact - vmin) * ones (4, npoly);
+
+  endfor
+
+  h = patch (hax, px, py, pz, pc);
+
+endfunction
+
+## 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)];
+  endif
+
+  N = N / norm(N);
+
+endfunction
+
+## 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)];
+
+endfunction
+
+## 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));
+
+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;
+%! streamtube (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]);
+%! 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");
+
+%!error streamtube ()
+