# HG changeset patch
# User HG200
# Date 1574506036 -3600
# Sat Nov 23 11:47:16 2019 +0100
# Node ID 57f0da8c90daa96224ebadbb30506f977bafe22e
# Parent b70da79e4ed9d62248a737f0f933970bc297e59b
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.
diff -r b70da79e4ed9 -r 57f0da8c90da doc/interpreter/plot.txi
--- a/doc/interpreter/plot.txi Fri Nov 29 10:41:03 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 b70da79e4ed9 -r 57f0da8c90da libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk Fri Nov 29 10:41:03 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 \
diff -r b70da79e4ed9 -r 57f0da8c90da libinterp/corefcn/stream-euler.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 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
+.
+
+*/
+
+/*
+
+ 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");
+}
+
diff -r b70da79e4ed9 -r 57f0da8c90da scripts/help/__unimplemented__.m
--- a/scripts/help/__unimplemented__.m Fri Nov 29 10:41:03 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",
diff -r b70da79e4ed9 -r 57f0da8c90da scripts/plot/draw/module.mk
--- a/scripts/plot/draw/module.mk Fri Nov 29 10:41:03 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 \
diff -r b70da79e4ed9 -r 57f0da8c90da 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
+## .
+
+## -*- 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
+
+ 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);
diff -r b70da79e4ed9 -r 57f0da8c90da 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
+## .
+
+## -*- 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");
+ 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
+
+ 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);
diff -r b70da79e4ed9 -r 57f0da8c90da 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,160 @@
+## 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
+## .
+
+## -*- 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 ()
diff -r b70da79e4ed9 -r 57f0da8c90da 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,374 @@
+## 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
+## .
+
+## -*- 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);
+ 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;
+ pmc = colscale * (vact - vmin) * ones (4, npoly);
+ pc(:, from:to) = pmc;
+
+ 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;
+%! 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");
+
+%!error streamtube ()
+