# HG changeset patch
# User chloros2
# Date 1571039510 -7200
#      Mon Oct 14 09:51:50 2019 +0200
# Node ID f78d71f81e38820fde8c5b58df3001c7f357ba14
# Parent  3fe26656e73c73570a0422beab46a97f0439bf49
user: chloros2
added libinterp/corefcn/stream-euler.cc Forward euler method to calculate a streamline
added scripts/plot/draw/stream2.m Calculating 2D stream lines
added scripts/plot/draw/stream3.m Calculating 3D stream lines
added scripts/plot/draw/streamline.m Wrapper and plot function for stream2, stream3
changed libinterp/corefcn/module.mk Wire in stream-euler.cc
changed scripts/plot/draw/module.mk Wire in .m files
changed scripts/plot/draw/module.mk Wire in .m files
changed /doc/interpreter/plot.txi Wire in documentation

diff -r 3fe26656e73c -r f78d71f81e38 doc/interpreter/plot.txi
--- a/doc/interpreter/plot.txi	Fri Oct 11 17:00:28 2019 -0400
+++ b/doc/interpreter/plot.txi	Mon Oct 14 09:51:50 2019 +0200
@@ -245,6 +245,12 @@
 
 @DOCSTRING(quiver3)
 
+@DOCSTRING(streamline)
+
+@DOCSTRING(stream2)
+
+@DOCSTRING(stream3)
+
 @DOCSTRING(compass)
 
 @DOCSTRING(feather)
diff -r 3fe26656e73c -r f78d71f81e38 libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk	Fri Oct 11 17:00:28 2019 -0400
+++ b/libinterp/corefcn/module.mk	Mon Oct 14 09:51:50 2019 +0200
@@ -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 3fe26656e73c -r f78d71f81e38 libinterp/corefcn/stream-euler.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/stream-euler.cc	Mon Oct 14 09:51:50 2019 +0200
@@ -0,0 +1,340 @@
+/*
+
+Copyright (C) 2019 Chloros2
+
+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/>.
+
+*/
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "ovl.h"
+
+#define IPO4(u11,u21,u12,u22,x,y) (u11 * (1-x) * (1-y) + \
+                                   u21 * x * (1-y) +     \
+                                   u12 * (1-x) * y +     \
+                                   u22 * x * y)
+
+
+#define IPO8(u111,u211,u121,u221,u112,u212,u122,u222,x,y,z) (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)
+
+// Is there an API function for "dst = src(1:nverts, :)" ?
+Matrix
+cropM (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);
+}
+
+void
+euler2d (const RowVector& xx, const RowVector& yy,
+         const octave_idx_type cols, const octave_idx_type rows,
+         const Matrix& u, const Matrix& v,
+         double xp, double yp,
+         signed long idx, signed long idy,
+         double step, octave_idx_type maxverts,
+         Matrix& buffer, octave_idx_type *nverts)
+{
+  octave_idx_type i = 0;
+
+  buffer(i, 0) = xp;
+  buffer(i, 1) = yp;
+
+  while (1)
+    {
+      // exit if we leave the definition set
+      if ((idx < 0) || (idx >= (cols - 1)) || (idy < 0) || (idy >= (rows - 1)))
+        {
+          break;
+        }
+      if ((xx(idx+1) == xx(idx)) || (yy(idy+1) == yy(idy)))
+        {
+          break;
+        }
+
+      // canonical coordinates of current point position
+      double nnx = 1.0 / (xx(idx+1) - xx(idx));
+      double nny = 1.0 / (yy(idy+1) - yy(idy));
+      double zeta = (xp - xx(idx)) * nnx;
+      double xi = (yp - yy(idy)) * nny;
+
+      // linear interpolation onto the point (xp, yp)
+      double vpx = IPO4(u(idy, idx), u(idy, idx+1), u(idy+1, idx), u(idy+1, idx+1),
+                        zeta, xi);
+      double vpy = IPO4(v(idy, idx), v(idy, idx+1), v(idy+1, idx), v(idy+1, idx+1),
+                        zeta, xi);
+
+      // hit on singularity
+      if ((octave::math::isnan (vpx) || octave::math::isnan (vpy)) ||
+          ((vpx == 0) && (vpy == 0)))
+        {
+          break;
+        }
+
+      // take inner derivative
+      double npx = vpx * nnx;
+      double npy = vpy * nny;
+      // scale to step size
+      double norm = 1.0 / sqrt (npx * npx + npy * npy);
+      double dnpx = npx * step * norm;
+      double dnpy = npy * step * norm;
+
+      // next point
+      double xpnxt = xp + dnpx * (xx(idx+1) - xx(idx));
+      double ypnxt = yp + dnpy * (yy(idy+1) - yy(idy));
+      idx = idx + floor ((xpnxt - xx(idx)) * nnx);
+      idy = idy + floor ((ypnxt - yy(idy)) * nny);
+      xp = xpnxt;
+      yp = ypnxt;
+
+      // collect trajectory segments
+      i++;
+      buffer(i, 0) = xp;
+      buffer(i, 1) = yp;
+
+      // terminate at upper level of field line segments
+      if (i + 1 >= maxverts)
+        {
+          break;
+        }
+    }
+
+  *nverts = i + 1;
+}
+
+void
+euler3d (const RowVector& xx, const RowVector& yy, const RowVector& zz,
+         int dX, int dY, int dZ,
+         NDArray& u, NDArray& v, NDArray& w,
+         double xp, double yp, double zp,
+         signed long idx, signed long idy, signed long idz,
+         double step, octave_idx_type maxverts, Matrix& buffer,
+         octave_idx_type *nverts){
+
+  octave_idx_type i = 0;
+
+  buffer(i, 0) = xp;
+  buffer(i, 1) = yp;
+  buffer(i, 2) = zp;
+
+  while (1)
+    {
+      // exit if we leave the definition set
+      if ((idx < 0) || (idx >= (dX-1)) || (idy < 0) || (idy >= (dY-1)) ||
+          (idz < 0) || (idz >= (dZ-1)) )
+        {
+          break;
+        }
+      if ((xx(idx+1) == xx(idx)) || (yy(idy+1) == yy(idy)) ||
+          (zz(idz+1) == zz(idz)))
+        {
+          break;
+        }
+
+      // canonical coordinates of current point position
+      double nnx = 1.0 / (xx(idx+1) - xx(idx));
+      double nny = 1.0 / (yy(idy+1) - yy(idy));
+      double nnz = 1.0 / (zz(idz+1) - zz(idz));
+      double zeta = (xp - xx(idx)) * nnx;
+      double xi = (yp - yy(idy)) * nny;
+      double rho = (zp - zz(idz)) * nnz;
+
+      // linear interpolation onto the point (xp,yp,zp)
+      double vpx = IPO8(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),
+                        zeta, xi, rho);
+      double vpy = IPO8(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),
+                        zeta, xi, rho);
+      double vpz = IPO8(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),
+                        zeta, xi, rho);
+
+      // hit on singularity
+      if ((octave::math::isnan (vpx) || octave::math::isnan (vpy) ||
+           octave::math::isnan (vpz)) || ((vpx == 0) && (vpy == 0) && (vpz == 0)))
+        {
+          break;
+        }
+
+      // take inner derivative
+      double npx = vpx * nnx;
+      double npy = vpy * nny;
+      double npz = vpz * nnz;
+      // scale to step size
+      double norm = 1.0 / sqrt(npx * npx + npy * npy + npz * npz);
+      double dnpx = npx * step * norm;
+      double dnpy = npy * step * norm;
+      double dnpz = npz * step * norm;
+
+      // next point
+      double xpnxt = xp + dnpx * (xx(idx+1) - xx(idx));
+      double ypnxt = yp + dnpy * (yy(idy+1) - yy(idy));
+      double zpnxt = zp + dnpz * (zz(idz+1) - zz(idz));
+      idx = idx + floor ((xpnxt - xx(idx)) * nnx);
+      idy = idy + floor ((ypnxt - yy(idy)) * nny);
+      idz = idz + floor ((zpnxt - zz(idz)) * nnz);
+      xp = xpnxt;
+      yp = ypnxt;
+      zp = zpnxt;
+
+      // collect trajectory segments
+      i++;
+      buffer(i, 0) = xp;
+      buffer(i, 1) = yp;
+      buffer(i, 2) = zp;
+
+      // terminate at upper level of field line segments
+      if (i + 1 >= maxverts)
+        {
+          break;
+        }
+    }
+
+  *nverts = i + 1;
+}
+
+static octave_value
+streameuler2d (const octave_value_list& args)
+{
+
+  int nargin = args.length ();
+  if (nargin != 10)
+    print_usage ();
+
+  RowVector X = args(0).row_vector_value ();
+  RowVector Y = args(1).row_vector_value ();
+
+  Matrix U = args(2).matrix_value ();
+  Matrix V = args(3).matrix_value ();
+
+  double xp = args(4).double_value ();
+  double yp = args(5).double_value ();
+
+  double idx = args(6).double_value ();
+  double idy = args(7).double_value ();
+
+  double step = args(8).double_value ();
+  double maxverts = args(9).double_value ();
+
+  octave_idx_type rows = U.rows ();
+  octave_idx_type cols = U.columns ();
+
+  octave_idx_type nverts;
+  Matrix buffer (maxverts, 2);
+
+  euler2d (X, Y, cols, rows, U, V, xp, yp, idx, idy, step, maxverts,
+           buffer, &nverts);
+
+  Matrix verts = cropM (buffer, nverts, 2);
+
+  return  octave_value (verts);
+}
+
+static octave_value
+streameuler3d (const octave_value_list& args, const char *fcn)
+{
+
+  int nargin = args.length ();
+  if (nargin != 14)
+    print_usage ();
+
+  RowVector X = args(0).row_vector_value ();
+  RowVector Y = args(1).row_vector_value ();
+  RowVector Z = args(2).row_vector_value ();
+
+  NDArray U = args(3).array_value ();
+  NDArray V = args(4).array_value ();
+  NDArray W = args(5).array_value ();
+
+  double xp = args(6).double_value ();
+  double yp = args(7).double_value ();
+  double zp = args(8).double_value ();
+
+  double idx = args(9).double_value ();
+  double idy = args(10).double_value ();
+  double idz = args(11).double_value ();
+
+  double step = args(12).double_value ();
+  double maxverts = args(13).double_value ();
+
+  dim_vector dims = args(3).dims ();
+  int ndims = dims.ndims ();
+  if (ndims != 3)
+    {
+      error ("%s: dimension must be 3", fcn);
+    }
+
+  octave_idx_type nverts;
+  Matrix buffer (maxverts, 3);
+
+  euler3d (X, Y, Z, dims(1), dims(0), dims(2), U, V, W, xp, yp, zp,
+           idx, idy, idz, step, maxverts, buffer, &nverts);
+
+  Matrix verts = cropM (buffer, nverts, 3);
+
+  return  octave_value (verts);
+}
+
+DEFUN (streameuler2d, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} streameuler2d (@var{X}, @var{Y}, @var{U}, @var{V}, @var{xp}, @var{yp}, @var{idx}, @var{idy},  @var{step}, @var{maxverts})
+Compute the streamline starting from a point [@var{xp}, @var{yp}]
+by solving the trajectory equation using an explicit forward euler
+method.
+
+@seealso{streamline, stream2, stream3, streameuler3d}
+@end deftypefn */)
+{
+  return streameuler2d (args);
+}
+
+DEFUN (streameuler3d, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {} streameuler3d (@var{X}, @var{Y}, @var{Z}, @var{U}, @var{V}, @var{W}, @var{xp}, @var{yp}, @var{zp}, @var{idx}, @var{idy}, @var{idz}, @var{step}, @var{maxverts})
+Compute the streamline starting from a point [@var{xp}, @var{yp}, @var{zp}]
+by solving the trajectory equation using an explicit forward euler
+method.
+
+@seealso{streamline, stream2, stream3, streameuler2d}
+@end deftypefn */)
+{
+  return streameuler3d (args, "streameuler3d");
+}
+
diff -r 3fe26656e73c -r f78d71f81e38 scripts/plot/draw/module.mk
--- a/scripts/plot/draw/module.mk	Fri Oct 11 17:00:28 2019 -0400
+++ b/scripts/plot/draw/module.mk	Mon Oct 14 09:51:50 2019 +0200
@@ -94,6 +94,9 @@
   %reldir%/stem.m \
   %reldir%/stem3.m \
   %reldir%/stemleaf.m \
+  %reldir%/stream2.m \
+  %reldir%/stream3.m \
+  %reldir%/streamline.m \
   %reldir%/surf.m \
   %reldir%/surface.m \
   %reldir%/surfc.m \
diff -r 3fe26656e73c -r f78d71f81e38 scripts/plot/draw/stream2.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/stream2.m	Mon Oct 14 09:51:50 2019 +0200
@@ -0,0 +1,156 @@
+## Copyright (C) 2019 Chloros2
+##
+## 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{}, "options")
+## Compute 2D streamline data
+##
+## Calculates the stream-/field lines of a vector field given by
+## [@var{u}, @var{v}] which is defined over a rectangular grid given
+## by [@var{x}, @var{y}]. The streamlines start from the starting 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}]. It contains an increment
+## (default 0.1) for the iteration of the trajectory based on canonical
+## coordinates. The second parameter specifies an upper limit of
+## streamline segments (default 10000) gathered per 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 return value @var{xy} is a @nospell{nverts x 2} matrix containing the
+## coordinates of the field line segments.
+##
+## Trajectory integration is done by an explicite forward euler method.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y] = meshgrid (0:3, 0:3);
+## u = 2 * x;
+## v = y;
+## verts = stream2 (x, y, u, v, 1.0, 0.5);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream3}
+##
+## @end deftypefn
+
+## Author: Chloros2 <chloros2@gmx.de>
+
+function xy = stream2 (varargin)
+
+  options = [];
+  switch (nargin)
+    case (0)
+      print_usage ();
+    case {4,5}
+      if (nargin == 4)
+        [u,v,sx,sy] = varargin{:};
+      else
+        [u,v,sx,sy,options] = varargin{:};
+      end
+      [m,n] = size (u);
+      [x,y] = meshgrid (1:n, 1:m);
+    case (6)
+      [x,y,u,v,sx,sy] = varargin{:};
+    case (7)
+      [x,y,u,v,sx,sy,options] = varargin{:};
+    otherwise
+      error ('stream2: unknown input parameter count');
+  end
+
+  if ((~ isempty (options)) && (numel (options) == 2))
+    step = options(1);
+    maxverts = options(2);
+  elseif ((~ isempty (options)) && (numel (options) == 1))
+    step = options(1);
+    maxverts = 10000;
+  else
+    step = 0.1;
+    maxverts = 10000;
+  end
+
+  if (step == 0)
+    error ('stream2: step == 0 not allowed');
+  end
+
+  if (~ (isequal (size (u), size (v), size (x), size (y)) && ...
+         isequal (size (sx), size (sy))) )
+    error ('stream2: matrix dimensions must match');
+  end
+
+  # grid coordinates
+  xx = x(1,:);
+  yy = y(:,1)';
+  nvert = 0;
+  xy = [];
+  [s1, s2] = size (sx);
+  for i = 1:s1
+    for j = 1:s2
+      nvert = nvert+1;
+      xp = sx(i,j);
+      yp = sy(i,j);
+      # get index of starting point and ensure its inside the definition set
+      idx = find (xx<xp, 1, 'last');
+      idy = find (yy<yp, 1, 'last');
+      if (yy(1) == yp) idy = 1; end
+      if (xx(1) == xp) idx = 1; end
+      if (isempty (idx) || isempty (idy) || (idx == numel (xx)) || (idy == numel (yy)))
+        xy{nvert} = [];
+      else
+        verts = streameuler2d (xx, yy, u, v, xp, yp, idx-1, idy-1, step, maxverts);
+        xy{nvert} = verts;
+      end
+    end
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y] = meshgrid (-2:0.2:2, -2:0.2:2);
+%! u = -y - x/2;
+%! v = x - y/2;
+%! [sx, sy] = meshgrid (-2:2:2, -2:2:2);
+%! xy = stream2 (x, y, u, v, sx, sy);
+%! h = [];
+%! for k = 1:length (xy);
+%!   n = xy{k};
+%!   if (~ isempty (n))
+%!     h = [h ; line('xdata', n(:,1), 'ydata', n(:,2))];
+%!   end
+%! end
+%! set (h, 'Color', 'r');
+%! hold on;
+%! plot (sx, sy, '.r', 'MarkerSize', 15);
+%! quiver (x, y, u, v);
+%! title ('Spiral sink');
+%! grid on;
+%! 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);
\ No newline at end of file
diff -r 3fe26656e73c -r f78d71f81e38 scripts/plot/draw/stream3.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/stream3.m	Mon Oct 14 09:51:50 2019 +0200
@@ -0,0 +1,165 @@
+## Copyright (C) 2019 Chloros2
+##
+## 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{}, "options")
+## Compute 3D streamline data.
+##
+## Calculates the stream-/field lines of a vector field given by
+## [@var{u}, @var{v}, @var{w}] which is defined over a rectangular grid
+## given by [@var{x}, @var{y}, @var{z}]. The streamlines start from the
+## starting 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}]. It contains an increment
+## (default 0.1) for the iteration of the trajectory based on canonical
+## coordinates. The second parameter specifies an upper limit of
+## streamline segments (default 10000) gathered per 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 return value @var{xyz} is a @nospell{nverts x 3} matrix containing the
+## coordinates of the field line segments.
+##
+## Trajectory integration is done by an explicite forward euler method.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y, z] = meshgrid (0:3, 0:3, 0:3);
+## u = 2 * x;
+## v = y;
+## w = 3 * z;
+## verts = stream3 (x, y, z, u, v, w, 1.0, 0.5, 0.0);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream2}
+##
+## @end deftypefn
+
+## Author: Chloros2 <chloros2@gmx.de>
+
+function xyz = stream3 (varargin)
+
+  options = [];
+  switch (nargin)
+    case (0)
+      print_usage ();
+    case {6,7}
+      if (nargin == 6)
+        [u,v,w,sx,sy,sz] = varargin{:};
+      else
+        [u,v,w,sx,sy,sz,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,sx,sy,sz] = varargin{:};
+    case (10)
+      [x,y,z,u,v,w,sx,sy,sz,options] = varargin{:};
+    otherwise
+      error ('stream3: unknown input parameter count');
+  end
+
+  if ((~ isempty (options)) && (numel (options) == 2))
+    step = options(1);
+    maxverts = options(2);
+  elseif ((~ isempty (options)) && (numel (options) == 1))
+    step = options(1);
+    maxverts = 10000;
+  else
+    step = 0.1;
+    maxverts = 10000;
+  end
+
+  if (step == 0)
+    error ('stream3: step == 0 not allowed');
+  end
+
+  if (~ (isequal (size (u), size (v), size (w), size (x), size (y), size (z)) && ...
+         isequal (size (sx), size (sy))) )
+    error ('stream3: matrix dimensions must match');
+  end
+
+  # grid coordinates
+  xx = x(1,:,1);
+  yy = y(:,1,1)';
+  tmp = z(1,1,:);
+  zz = tmp(:)';
+  ssx = sx(:);
+  ssy = sy(:);
+  ssz = sz(:);
+  xyz = [];
+  for nvert = 1:length (ssx)
+    xp = ssx(nvert);
+    yp = ssy(nvert);
+    zp = ssz(nvert);
+    # get index of starting point and ensure its inside the definition set
+    idx = find (xx<xp, 1, 'last');
+    idy = find (yy<yp, 1, 'last');
+    idz = find (zz<zp, 1, 'last');
+    if (yy(1) == yp) idy = 1; end
+    if (xx(1) == xp) idx = 1; end
+    if (zz(1) == zp) idz = 1; end
+    if (isempty (idx) || isempty (idy) || isempty (idz) || ...
+       (idx == numel (xx)) || (idy == numel (yy)) || (idz == numel (zz)))
+      xyz{nvert} = [];
+    else
+      verts = streameuler3d (xx, yy, zz, u, v, w, xp, yp, zp, idx-1, idy-1, ...
+                             idz-1, step, maxverts);
+      xyz{nvert} = verts;
+    end
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y, z] = meshgrid (-3:3, -3:3, -3:3);
+%! u = -y - x/2;
+%! v = x - y/2;
+%! w = -z;
+%! [sx, sy, sz] = meshgrid (2, 0:1.5:1.5, 0:1.5:3);
+%! xyz = stream3 (x, y, z, u, v, w, sx, sy, sz);
+%! h = [];
+%! for k = 1:length(xyz);
+%!   n = xyz{k};
+%!   if (~ isempty (n))
+%!     h = [h ; line('xdata', n(:,1), 'ydata', n(:,2), 'zdata', n(:,3))];
+%!   end
+%! end
+%! 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 ('Critical spiral point in 3D');
+%! 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);
\ No newline at end of file
diff -r 3fe26656e73c -r f78d71f81e38 scripts/plot/draw/streamline.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/streamline.m	Mon Oct 14 09:51:50 2019 +0200
@@ -0,0 +1,156 @@
+## Copyright (C) 2019 Chloros2
+##
+## 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{u}, @var{v}, @var{sx}, @var{sy})
+## @deftypefnx {} {} streamline (@var{u}, @var{v}, @var{sx}, @var{sy})
+## @deftypefnx {} {} 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 (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
+## @deftypefnx {} {} streamline (@dots{}, "options")
+## @deftypefnx {} {} streamline (@var{hax}, @dots{})
+## @deftypefnx {} {[@var{h}] =} streamline (@dots{})
+## Plot streamlines for 2D or 3D vector fields
+##
+## Plot stream-/field lines 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 from the starting
+## 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}]. It contains an increment
+## (default 0.1) for the iteration of the trajectory based on canonical
+## coordinates. The second parameter specifies an upper limit of
+## streamline segments (default 10000) gathered per 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 hggroup
+## comprising the field lines.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y] = meshgrid (-2:0.2:2, -2:0.2:2);
+## u = -y - x/2;
+## v = x - y/2;
+## [sx, sy] = meshgrid (-2:2:2, -2:2:2);
+## h = streamline (x, y, u, v, sx, sy);
+## @end group
+## @end example
+##
+## @seealso{stream2, stream3}
+##
+## @end deftypefn
+
+## Author: Chloros2 <chloros2@gmx.de>
+
+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);
+  # check dimension and go
+  switch (numel (size (argval{:})))
+    case (2)
+      verts = stream2 (varargin{:});
+      for i = 1:length (verts)
+        n = verts{i};
+        if (~ isempty (n))
+          h = [h; line('xdata', n(:,1), 'ydata', n(:,2), 'color', [0 0 1], ...
+                       'parent', hax)];
+        end
+      end
+    case (3)
+      verts = stream3 (varargin{:});
+      for i = 1:length (verts)
+        n = verts{i};
+        if (~ isempty (n))
+          h = [h; line('xdata', n(:,1), 'ydata', n(:,2), 'zdata', n(:,3), ...
+                       'color', [0 0 1], 'parent', hax)];
+        end
+      end
+    otherwise
+      error ('streamline: either 2D or 3D input data');
+  end
+
+end
+
+%!demo
+%! clf;
+%! [x, y] = meshgrid (-2:0.2:2, -2:0.2:2);
+%! u = -y - x/2;
+%! v = x - y/2;
+%! [sx, sy] = meshgrid (-2:2:2, -2:2:2);
+%! h = streamline (x, y, u, v, sx, sy);
+%! set (h, 'Color', 'r');
+%! hold on;
+%! plot (sx, sy, '.r', 'MarkerSize', 15);
+%! quiver (x, y, u, v);
+%! title ('Spiral sink');
+%! grid on;
+%! axis equal;
+
+%!demo
+%! [x, y, z] = meshgrid (-3:3, -3:3, -3:3);
+%! u = -y - x/2;
+%! v = x - y/2;
+%! w = -z;
+%! [sx, sy, sz] = meshgrid (2, 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 ('Critical spiral point in 3D');
+%! grid on;
+%! axis equal;
+
+%!demo
+%! [x, y, z] = meshgrid (-3:3, -3:3, -3:3);
+%! u = -x;
+%! v = -y;
+%! w = -z;
+%! [sx, sy, sz] = meshgrid (3, -3:3:3, -3:3:3);
+%! h = streamline (x, y, z, u, v, w, sx, sy, sz, [0.1, 1000]);
+%! 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 ('Vanishing velocity at center point');
+%! grid on;
+%! axis equal;
+
+%!test
+%!error streamline ()
\ No newline at end of file
