# 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
+.
+
+*/
+
+#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
+## .
+
+## -*- 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
+
+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.
+
+## -*- 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
+
+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.
+
+## -*- 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
+
+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