--- a/NEWS	Thu Apr 02 12:22:56 2020 +0200
+++ b/NEWS	Sat Apr 04 07:29:49 2020 +0200
@@ -104,6 +104,7 @@
 * `memory`
 * `rng`
 * `startsWith`
+* `streamribbon`
 
 ### Deprecated functions and properties
 

--- a/doc/interpreter/plot.txi	Thu Apr 02 12:22:56 2020 +0200
+++ b/doc/interpreter/plot.txi	Sat Apr 04 07:29:49 2020 +0200
@@ -245,6 +245,8 @@
 
 @DOCSTRING(quiver3)
 
+@DOCSTRING(streamribbon)
+
 @DOCSTRING(streamtube)
 
 @DOCSTRING(ostreamtube)

--- a/scripts/help/__unimplemented__.m	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/help/__unimplemented__.m	Sat Apr 04 07:29:49 2020 +0200
@@ -1188,7 +1188,6 @@
   "stopasync",
   "str2mat",
   "streamparticles",
-  "streamribbon",
   "streamslice",
   "string",
   "strings",

--- a/scripts/plot/draw/module.mk	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/plot/draw/module.mk	Sat Apr 04 07:29:49 2020 +0200
@@ -99,6 +99,7 @@
   %reldir%/stream2.m \
   %reldir%/stream3.m \
   %reldir%/streamline.m \
+  %reldir%/streamribbon.m \
   %reldir%/streamtube.m \
   %reldir%/surf.m \
   %reldir%/surface.m \

--- a/scripts/plot/draw/ostreamtube.m	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/plot/draw/ostreamtube.m	Sat Apr 04 07:29:49 2020 +0200
@@ -71,7 +71,7 @@
 ## @end group
 ## @end example
 ##
-## @seealso{stream3, streamline, streamtube}
+## @seealso{stream3, streamline, streamribbon, streamtube}
 ## @end deftypefn
 
 ## References:

--- a/scripts/plot/draw/stream3.m	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/plot/draw/stream3.m	Sat Apr 04 07:29:49 2020 +0200
@@ -58,7 +58,7 @@
 ## @end group
 ## @end example
 ##
-## @seealso{stream2, streamline, streamtube, ostreamtube}
+## @seealso{stream2, streamline, streamribbon, streamtube, ostreamtube}
 ## @end deftypefn
 
 ## References:

--- a/scripts/plot/draw/streamline.m	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/plot/draw/streamline.m	Sat Apr 04 07:29:49 2020 +0200
@@ -61,7 +61,7 @@
 ## @end group
 ## @end example
 ##
-## @seealso{stream2, stream3, streamtube, ostreamtube}
+## @seealso{stream2, stream3, streamribbon, streamtube, ostreamtube}
 ## @end deftypefn
 
 function h = streamline (varargin)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/streamribbon.m	Sat Apr 04 07:29:49 2020 +0200
@@ -0,0 +1,444 @@
+########################################################################
+##
+## Copyright (C) 2020 The Octave Project Developers
+##
+## See the file COPYRIGHT.md in the top-level directory of this
+## distribution or <https://octave.org/copyright/>.
+##
+## 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  {} {} streamribbon (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
+## @deftypefnx {} {} streamribbon (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
+## @deftypefnx {} {} streamribbon (@var{xyz}, @var{x}, @var{y}, @var{z}, @var{anlr_spd}, @var{lin_spd})
+## @deftypefnx {} {} streamribbon (@var{xyz}, @var{anlr_spd}, @var{lin_spd})
+## @deftypefnx {} {} streamribbon (@var{xyz}, @var{anlr_rot})
+## @deftypefnx {} {} streamribbon (@dots{}, @var{width})
+## @deftypefnx {} {} streamribbon (@var{hax}, @dots{})
+## @deftypefnx {} {@var{h} =} streamribbon (@dots{})
+## Calculate and display streamribbons.
+##
+## The streamribbon is constructed by rotating a normal vector around a
+## streamline according to the angular rotation of the vector field.
+##
+## The vector field is given by @code{[@var{u}, @var{v}, @var{w}]} and is
+## defined over a rectangular grid given by @code{[@var{x}, @var{y}, @var{z}]}.
+## The streamribbons start at the seed points
+## @code{[@var{sx}, @var{sy}, @var{sz}]}.
+##
+## @code{streamribbon} can be called with a cell array that contains
+## pre-computed streamline data.  To do this, @var{xyz} must be created with
+## the @code{stream3} function.  @var{lin_spd} is the linear speed of the
+## vector field and can be calculated from @code{[@var{u}, @var{v}, @var{w}]}
+## by the square root of the sum of the squares.  The angular speed
+## @var{anlr_spd} is the projection of the angular velocity onto the velocity
+## of the normalized vector field and can be calculated with the @code{curl}
+## command.  This option is useful if you need to alter the integrator step
+## size or the maximum number of streamline vertices.
+##
+## Alternatively, ribbons can be created from an array of vertices @var{xyz} of
+## a path curve.  @var{anlr_rot} contains the angles of rotation around the
+## edges between adjacent vertices of the path curve.
+##
+## The input parameter @var{width} sets the width of the streamribbons.
+##
+## Streamribbons are colored according to the total angle of rotation along the
+## ribbon.
+##
+## 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 plot objects
+## created for each streamribbon.
+##
+## Example:
+##
+## @example
+## @group
+## [x, y, z] = meshgrid (0:0.2:4, -1:0.2:1, -1:0.2:1);
+## u = - x + 10;
+## v = 10 * z.*x;
+## w = - 10 * y.*x;
+## streamribbon (x, y, z, u, v, w, [0, 0], [0, 0.6], [0, 0]);
+## view (3);
+## @end group
+## @end example
+##
+## @seealso{streamline, stream3, streamtube, ostreamtube}
+##
+## @end deftypefn
+
+## References:
+##
+## @inproceedings{
+##    title = {Feature Detection from Vector Quantities in a Numerically Simulated Hypersonic Flow Field in Combination with Experimental Flow Visualization},
+##    author = {Pagendarm, Hans-Georg and Walter, Birgit},
+##    year = {1994},
+##    publisher = {IEEE Computer Society Press},
+##    booktitle = {Proceedings of the Conference on Visualization ’94},
+##    pages = {117–123},
+## }
+##
+## @article{
+##    title = {Efficient streamline, streamribbon, and streamtube constructions on unstructured grids},
+##    author = {Ueng, Shyh-Kuang and Sikorski, C. and Ma, Kwan-Liu},
+##    year = {1996},
+##    month = {June},
+##    publisher = {IEEE Transactions on Visualization and Computer Graphics},
+## }
+##
+## @inproceedings{
+##    title = {Visualization of 3-D vector fields - Variations on a stream},
+##    author = {Dave Darmofal and Robert Haimes},
+##    year = {1992}
+## }
+##
+## @techreport{
+##    title = {Parallel Transport Approach to Curve Framing},
+##    author = {Andrew J. Hanson and Hui Ma},
+##    year = {1995}
+## }
+##
+## @article{
+##    title = {There is More than One Way to Frame a Curve},
+##    author = {Bishop, Richard},
+##    year = {1975},
+##    month = {03},
+##    volume = {82},
+##    publisher = {The American Mathematical Monthly}
+## }
+
+function h = streamribbon (varargin)
+
+  [hax, varargin, nargin] = __plt_get_axis_arg__ ("streamribbon", varargin{:});
+
+  width = [];
+  xyz = [];
+  anlr_spd = [];
+  lin_spd = [];
+  anlr_rot = [];
+  switch (nargin)
+    case 0
+      print_usage ();
+    case 2
+      [xyz, anlr_rot] = varargin{:};
+    case 3
+      if (numel (varargin{3}) == 1)
+        [xyz, anlr_rot, width] = varargin{:};
+      else
+        [xyz, anlr_spd, lin_spd] = varargin{:};
+        [m, n, p] = size (anlr_spd);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 4
+      [xyz, anlr_spd, lin_spd, width] = varargin{:};
+      [m, n, p] = size (anlr_spd);
+      [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+    case 6
+      if (iscell (varargin{1}))
+        [xyz, x, y, z, anlr_spd, lin_spd] = varargin{:};
+      else
+        [u, v, w, spx, spy, spz] = varargin{:};
+        [m, n, p] = size (u);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 7
+      if (iscell (varargin{1}))
+        [xyz, x, y, z, anlr_spd, lin_spd, width] = varargin{:};
+      else
+        [u, v, w, spx, spy, spz, width] = varargin{:};
+        [m, n, p] = size (u);
+        [x, y, z] = meshgrid (1:n, 1:m, 1:p);
+      endif
+    case 9
+      [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
+    case 10
+      [x, y, z, u, v, w, spx, spy, spz, width] = varargin{:};
+    otherwise
+      error ("streamribbon: invalid number of inputs");
+  endswitch
+
+  if (isempty (xyz))
+    xyz = stream3 (x, y, z, u, v, w, spx, spy, spz);
+    anlr_spd = curl (x, y, z, u, v, w);
+    lin_spd = sqrt (u.*u + v.*v + w.*w);
+  endif
+
+  ## Derive scale factor from the bounding box diagonal
+  if (isempty (width))
+    mxx = mnx = mxy = mny = mxz = mnz = [];
+    j = 1;
+    for i = 1 : length (xyz)
+      sl = xyz{i};
+      if (! isempty (sl))
+        slx = sl(:,1); sly = sl(:,2); slz = sl(:,3);
+        mxx(j) = max (slx); mnx(j) = min (slx);
+        mxy(j) = max (sly); mny(j) = min (sly);
+        mxz(j) = max (slz); mnz(j) = min (slz);
+        j += 1;
+      end
+    endfor
+    dx = max (mxx) - min (mnx);
+    dy = max (mxy) - min (mny);
+    dz = max (mxz) - min (mnz);
+    width = sqrt (dx*dx + dy*dy + dz*dz) / 25;
+  elseif (! isreal (width) || width <= 0)
+    error ("streamribbon: WIDTH must be a real scalar > 0");
+  endif
+
+  if (! isempty (anlr_rot))
+    for i = 1 : length (xyz)
+      if (rows (anlr_rot{i}) != rows (xyz{i}))
+        error ("streamribbon: ANLR_ROT must have same length as XYZ");
+      endif
+    endfor
+  endif
+
+  if (isempty (hax))
+    hax = gca ();
+  else
+    hax = hax(1);
+  endif
+
+  ## Angular speed of a paddle wheel spinning around a streamline in a fluid
+  ## flow "V":
+  ## dtheta/dt = 0.5 * <curl(V), V/norm(V)>
+  ##
+  ## Integration along a streamline segment with the length "h" yields the
+  ## rotation angle:
+  ## theta = 0.25 * h * <curl(V), V(0)/norm(V(0))^2) + V(h)/norm(V(h))^2)>
+  ##
+  ## Alternative approach using the curl angular speed "c = curl()":
+  ## theta = 0.5 * h * (c(0)/norm(V(0)) + c(h)/norm(V(h)))
+  ##
+  ## Hints:
+  ## i. ) For integration use trapezoidal rule
+  ## ii.) "V" can be assumend to be piecewise linear and curl(V) to be
+  ##      piecewise constant because of the used linear interpolation
+
+  h = [];
+  for i = 1 : length (xyz)
+    sl = xyz{i};
+    num_vertices = rows (sl);
+    if (! isempty (sl) && num_vertices > 1)
+      if isempty (anlr_rot)
+        ## Plot from vector field
+        ## Interpolate onto streamline vertices
+        [lin_spd_sl, anlr_spd_sl, max_vertices] = ...
+                                  interp_sl (x, y, z, lin_spd, anlr_spd, sl);
+        if (max_vertices > 1)
+          ## Euclidean distance between two adjacent vertices
+          stepsize = vecnorm (diff (sl(1:max_vertices, :)), 2, 2);
+          ## Angular rotation around edges between two adjacent sl-vertices
+          ## Note: Potential "division by zero" is checked in interp_sl()
+          anlr_rot_sl = 0.5 * stepsize.*(anlr_spd_sl(1:max_vertices - 1)./ ...
+                                         lin_spd_sl(1:max_vertices - 1) + ...
+                                         anlr_spd_sl(2:max_vertices)./ ...
+                                         lin_spd_sl(2:max_vertices));
+
+          htmp = plotribbon (hax, sl, anlr_rot_sl, max_vertices, 0.5 * width);
+          h = [h; htmp];
+        endif
+      else
+          ## Plot from vertice array
+          anlr_rot_sl = anlr_rot{i};
+
+          htmp = plotribbon (hax, sl, anlr_rot_sl, num_vertices, 0.5 * width);
+          h = [h; htmp];
+      endif
+    endif
+  endfor
+
+endfunction
+
+function h = plotribbon (hax, sl, anlr_rot_sl, max_vertices, width2)
+
+  total_angle = cumsum (anlr_rot_sl);
+  total_angle = [0; total_angle];
+
+  ## 1st streamline segment
+  X0 = sl(1,:);
+  X1 = sl(2,:);
+  R = X1 - X0;
+  RE = R / norm (R);
+
+  ## Initial vector KE which is to be transported along the vertice array
+  KE = get_normal2 (RE);
+  XS10 = - width2 * KE + X0;
+  XS20 = width2 * KE + X0;
+
+  ## Apply angular rotation
+  cp = cos (anlr_rot_sl(1));
+  sp = sin (anlr_rot_sl(1));
+  KE = rotation (KE, RE, cp, sp).';
+
+  XS1 = - width2 * KE + X1;
+  XS2 = width2 * KE + X1;
+
+  px = zeros (2, max_vertices);
+  py = zeros (2, max_vertices);
+  pz = zeros (2, max_vertices);
+  pc = zeros (2, max_vertices);
+
+  px(:,1) = [XS10(1); XS20(1)];
+  py(:,1) = [XS10(2); XS20(2)];
+  pz(:,1) = [XS10(3); XS20(3)];
+  pc(:,1) = total_angle(1) * [1; 1];
+
+  px(:,2) = [XS1(1); XS2(1)];
+  py(:,2) = [XS1(2); XS2(2)];
+  pz(:,2) = [XS1(3); XS2(3)];
+  pc(:,2) = total_angle(2) * [1; 1];
+
+  for i = 3 : max_vertices
+
+    ## Next streamline segment
+    X0 = X1;
+    X1 = sl(i,:);
+    R = X1 - X0;
+    RE = R / norm (R);
+
+    ## Project KE onto RE and get the difference in order to transport
+    ## the normal vector KE along the vertex array
+    Kp = KE - RE * dot (KE, RE);
+    KE = Kp / norm (Kp);
+
+    ## Apply angular rotation to KE
+    cp = cos (anlr_rot_sl(i - 1));
+    sp = sin (anlr_rot_sl(i - 1));
+    KE = rotation (KE, RE, cp, sp).';
+
+    XS1 = - width2 * KE + X1;
+    XS2 = width2 * KE + X1;
+
+    px(:,i) = [XS1(1); XS2(1)];
+    py(:,i) = [XS1(2); XS2(2)];
+    pz(:,i) = [XS1(3); XS2(3)];
+    pc(:,i) = total_angle(i) * [1; 1];
+
+  endfor
+
+  h = surface (hax, px, py, pz, pc);
+
+endfunction
+
+## Interpolate speed and divergence onto the streamline vertices and
+## return the first chunck of valid samples until a singularity /
+## zero is hit or the streamline vertex array "sl" ends
+function [lin_spd_sl, anlr_spd_sl, max_vertices] = ...
+                               interp_sl (x, y, z, lin_spd, anlr_spd, sl)
+
+  anlr_spd_sl = interp3 (x, y, z, anlr_spd, sl(:,1), sl(:,2), sl(:,3));
+  lin_spd_sl = interp3 (x, y, z, lin_spd, sl(:,1), sl(:,2), sl(:,3));
+
+  is_singular_anlr_spd = find (isnan (anlr_spd_sl), 1, "first");
+  is_zero_lin_spd = find (lin_spd_sl == 0, 1, "first");
+
+  max_vertices = rows (sl);
+  if (! isempty (is_singular_anlr_spd))
+    max_vertices = min (max_vertices, is_singular_anlr_spd - 1);
+  endif
+  if (! isempty (is_zero_lin_spd))
+    max_vertices = min (max_vertices, is_zero_lin_spd - 1);
+  endif
+
+endfunction
+
+## N normal to X, so that N is in span ([0 0 1], X)
+## If not possible then span ([1 0 0], X)
+function N = get_normal2 (X)
+
+  if ((X(1) == 0) && (X(2) == 0))
+    A = [1, 0, 0];
+  else
+    A = [0, 0, 1];
+  endif
+
+  ## Project A onto X and get the difference
+  N = A - X * dot (A, X) / (norm (X)^2);
+  N /= norm (N);
+
+endfunction
+
+## Rotate X around U where |U| = 1
+## cp = cos (angle), sp = sin (angle)
+function Y = rotation (X, U, cp, sp)
+
+  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 (0:0.2:4, -1:0.2:1, -1:0.2:1);
+%! u = - x + 10;
+%! v = 10 * z.*x;
+%! w = - 10 * y.*x;
+%! sx = [0, 0];
+%! sy = [0, 0.6];
+%! sz = [0, 0];
+%! streamribbon (x, y, z, u, v, w, sx, sy, sz);
+%! hold on;
+%! quiver3 (x, y, z, u, v, w);
+%! colormap (jet);
+%! shading interp;
+%! camlight ("headlight");
+%! view (3);
+%! axis tight equal off;
+%! set (gca, "cameraviewanglemode", "manual");
+%! hcb = colorbar;
+%! title (hcb, "Angle");
+%! title ("Streamribbon");
+
+%!demo
+%! clf;
+%! t = (0:pi/50:2*pi).';
+%! xyz{1} = [cos(t), sin(t), 0*t];
+%! twist{1} = ones (numel (t), 1) * pi / (numel (t) - 1);
+%! streamribbon (xyz, twist, 0.5);
+%! colormap (jet);
+%! view (3);
+%! camlight ("headlight");
+%! axis tight equal off;
+%! title ("Moebius Strip");
+
+## Test input validation
+%!error streamribbon ()
+%!error <invalid number of inputs> streamribbon (1)
+%!error <invalid number of inputs> streamribbon (1,2,3,4,5)
+%!error <invalid number of inputs> streamribbon (1,2,3,4,5,6,7,8)
+%!error <invalid number of inputs> streamribbon (1,2,3,4,5,6,7,8,9,10,11)
+%!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,1i)
+%!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,0)
+%!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,-1)
+%!error <ANLR_ROT must have same length as XYZ> streamribbon ({[1,1,1;2,2,2]},{[1,1,1]})

--- a/scripts/plot/draw/streamtube.m	Thu Apr 02 12:22:56 2020 +0200
+++ b/scripts/plot/draw/streamtube.m	Sat Apr 04 07:29:49 2020 +0200
@@ -66,7 +66,7 @@
 ## The optional return value @var{h} is a graphics handle to the plot objects
 ## created for each tube.
 ##
-## @seealso{stream3, streamline, ostreamtube}
+## @seealso{stream3, streamline, streamribbon, ostreamtube}
 ## @end deftypefn
 
 function h = streamtube (varargin)