stream2

0001 ########################################################################
0002 ##
0003 ## Copyright (C) 2019-2022 The Octave Project Developers
0004 ##
0005 ## Added to EIDORS 2022 by Andy Adler
0006 ##   -- modifications to fix matlab compatibilities
0007 ##
0008 ## See the file COPYRIGHT.md in the top-level directory of this
0009 ## distribution or <https://octave.org/copyright/>.
0010 ##
0011 ## This file is part of Octave.
0012 ##
0013 ## Octave is free software: you can redistribute it and/or modify it
0014 ## under the terms of the GNU General Public License as published by
0015 ## the Free Software Foundation, either version 3 of the License, or
0016 ## (at your option) any later version.
0017 ##
0018 ## Octave is distributed in the hope that it will be useful, but
0019 ## WITHOUT ANY WARRANTY; without even the implied warranty of
0020 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0021 ## GNU General Public License for more details.
0022 ##
0023 ## You should have received a copy of the GNU General Public License
0024 ## along with Octave; see the file COPYING.  If not, see
0025 ## <https://www.gnu.org/licenses/>.
0026 ##
0027 ########################################################################
0028 
0029 ## -*- texinfo -*-
0030 ## @deftypefn  {} {@var{xy} =} stream2 (@var{x}, @var{y}, @var{u}, @var{v}, @var{sx}, @var{sy})
0031 ## @deftypefnx {} {@var{xy} =} stream2 (@var{u}, @var{v}, @var{sx}, @var{sy})
0032 ## @deftypefnx {} {@var{xy} =} stream2 (@dots{}, @var{options})
0033 ## Compute 2-D streamline data.
0034 ##
0035 ## Calculates streamlines of a vector field given by @code{[@var{u}, @var{v}]}.
0036 ## The vector field is defined over a rectangular grid given by
0037 ## @code{[@var{x}, @var{y}]}.  The streamlines start at the seed points
0038 ## @code{[@var{sx}, @var{sy}]}.  The returned value @var{xy} contains a cell
0039 ## array of vertex arrays.  If the starting point is outside the vector field,
0040 ## @code{[]} is returned.
0041 ##
0042 ## The input parameter @var{options} is a 2-D vector of the form
0043 ## @code{[@var{stepsize}, @var{max_vertices}]}.  The first parameter
0044 ## specifies the step size used for trajectory integration (default 0.1).  A
0045 ## negative value is allowed which will reverse the direction of integration.
0046 ## The second parameter specifies the maximum number of segments used to
0047 ## create a streamline (default 10,000).
0048 ##
0049 ## The return value @var{xy} is a @nospell{nverts x 2} matrix containing the
0050 ## coordinates of the field line segments.
0051 ##
0052 ## Example:
0053 ##
0054 ## @example
0055 ## @group
0056 ## [x, y] = meshgrid (0:3);
0057 ## u = 2 * x;
0058 ## v = y;
0059 ## xy = stream2 (x, y, u, v, 1.0, 0.5);
0060 ## @end group
0061 ## @end example
0062 ##
0063 ## @seealso{streamline, stream3}
0064 ## @end deftypefn
0065 
0066 ## References:
0067 ##
0068 ## @article{
0069 ##    title = {Particle Tracing Algorithms for 3D Curvilinear Grids},
0070 ##    year = {2000},
0071 ##    author = {Nielson, Gregory and Uller, H. and Sadarjoen, I. and Walsum, Theo and Hin, Andrea and Post, Frits}
0072 ## }
0073 ##
0074 ## @article{
0075 ##    title = {Sources of error in the graphical analysis of CFD results},
0076 ##    publisher = {Journal of Scientific Computing},
0077 ##    year = {1988},
0078 ##    volume = {3},
0079 ##    number = {2},
0080 ##    pages = {149--164},
0081 ##    author = {Buning, Pieter G.},
0082 ## }
0083 
0084 function xy = stream2 (varargin)
0085 
0086   options = [];
0087   switch (numel (varargin))
0088     case {4,5}
0089       if (numel (varargin) == 4)
0090         [u, v, spx, spy] = varargin{:};
0091       else
0092         [u, v, spx, spy, options] = varargin{:};
0093       endif
0094       [m, n] = size (u);
0095       [x, y] = meshgrid (1:n, 1:m);
0096     case 6
0097       [x, y, u, v, spx, spy] = varargin{:};
0098       [x, u] = fix_to_matrix(x, u, 1);
0099       [y, v] = fix_to_matrix(y, v, 2);
0100      
0101     case 7
0102       [x, y, u, v, spx, spy, options] = varargin{:};
0103       [x, u] = fix_to_matrix(x, u, 1);
0104       [y, v] = fix_to_matrix(y, v, 2);
0105     otherwise
0106       print_usage ();
0107   endswitch
0108 
0109   stepsize = 0.1;
0110   max_vertices = 10_000;
0111   if (! isempty (options))
0112     switch (numel (options))
0113       case 1
0114         stepsize = options(1);
0115       case 2
0116         stepsize = options(1);
0117         max_vertices = options(2);
0118       otherwise
0119         error ("stream2: OPTIONS must be a 1- or 2-element vector");
0120     endswitch
0121 
0122     if (! isreal (stepsize) || stepsize == 0)
0123       error ("stream2: STEPSIZE must be a real scalar != 0");
0124     endif
0125     if (! isreal (max_vertices) || max_vertices < 1)
0126       error ("stream2: MAX_VERTICES must be an integer > 0");
0127     endif
0128     max_vertices = fix (max_vertices);
0129   endif
0130 
0131   if (! (size_equal (u, v, x, y) && size_equal (spx, spy)))
0132     error ("stream2: matrix dimensions must match");
0133   endif
0134   if (iscomplex (u) || iscomplex (v) || iscomplex (x) || iscomplex (y)
0135       || iscomplex (spx) || iscomplex (spy))
0136     error ("stream2: all inputs must be real-valued");
0137   endif
0138 
0139   gx = x(1,:);
0140   gy = y(:,1).';
0141 
0142   ## Jacobian Matrix
0143   dx = diff (gx);
0144   dy = diff (gy);
0145   ## "<" used to check if the mesh is ascending
0146   if (any (dx <= 0) || any (dy <= 0)
0147       || any (isnan (dx)) || any (isnan (dy)))
0148     error ("stream2: non-monotonically increasing or NaN values found in mesh");
0149   endif
0150   tx = 1 ./ dx;
0151   ty = 1 ./ dy;
0152   ## "Don't cares" used for handling points located on the border
0153   tx(end + 1) = 0;
0154   ty(end + 1) = 0;
0155   dx(end + 1) = 0;
0156   dy(end + 1) = 0;
0157 
0158   px = spx(:);
0159   py = spy(:);
0160 
0161   for nseed = 1 : numel (px)
0162 
0163     xp = px(nseed);
0164     yp = py(nseed);
0165     idx = find (diff (gx <= xp), 1);
0166     if (gx(end) == xp)
0167       idx = numel (gx);
0168     endif
0169     idy = find (diff (gy <= yp), 1);
0170     if (gy(end) == yp)
0171       idy = numel (gy);
0172     endif
0173 
0174     if (isempty (idx) || isempty (idy))
0175       xy{nseed} = [];
0176     else
0177       ## Transform seed from P coordinates to C coordinates
0178       zeta = (idx - 1) + (xp - gx(idx)) * tx(idx);
0179       xi = (idy - 1) + (yp - gy(idy)) * ty(idy);
0180 
0181       C = __streameuler2d__ (u, v, tx, ty, zeta, xi, stepsize, max_vertices);
0182 
0183       ## Transform from C coordinates to P coordinates
0184       idu = floor (C(:,1));
0185       idv = floor (C(:,2));
0186       xy{nseed} = [gx(idu + 1).' + (C(:,1) - idu).*(dx(idu + 1).'), ...
0187                    gy(idv + 1).' + (C(:,2) - idv).*(dy(idv + 1).')];
0188     endif
0189 
0190   endfor
0191 
0192 endfunction
0193 
0194 # check if x is vector. Make full size
0195 # axix (=1 for x and =2 for y)
0196 function [x, u] = fix_to_matrix(x, u, axis);
0197   if size_equal(x,u)
0198     return
0199   endif
0200   if all(size(x)>1)
0201     return
0202   endif
0203    
0204   dirn = diff( x );
0205   if all(dirn>0)
0206     doflip = false;
0207   elseif all(dirn<0)
0208     doflip = true;
0209   else
0210     error ("stream2: axis not monotonic");
0211   end
0212   
0213   # extend x in the rows
0214   x = x(:).';
0215   x = x(ones(size(u,axis),1),:);
0216   if axis==2
0217     x = x.';
0218   endif
0219 
0220   if doflip % flipud if axis=2, rl if axis=1
0221     x = flip( x, 3-axis );
0222     u = flip( u, 3-axis );
0223   endif
0224 endfunction
0225 
0226 %!demo
0227 %! clf;
0228 %! [x, y] = meshgrid (-5:5, -4:4);
0229 %! u = x - 2 * y;
0230 %! v = 2 * x - 3 * y;
0231 %! sx = [3, 0, -1, -2, -3, 0, 1, 2];
0232 %! sy = [3, 3, 3, 3, -3, -3, -3, -3];
0233 %! h = streamline (x, y, u, v, sx, sy, 0.05);
0234 %! set (h, "color", "r");
0235 %! hold on;
0236 %! quiver (x, y, u, v);
0237 %! scatter (sx(:), sy(:), 20, "filled", "o", "markerfacecolor", "r");
0238 %! grid on;
0239 %! title ("Asymptotically Stable Equilibrium");
0240 %! axis equal;
0241 
0242 %!test
0243 %! 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]);
0244 %! assert (numel (xy{:}), 10);
0245 
0246 ## Test input validation
0247 %!error <Invalid call> stream2 ()
0248 %!error <Invalid call> stream2 (1)
0249 %!error <Invalid call> stream2 (1,2)
0250 %!error <Invalid call> stream2 (1,2,3)
0251 %!error <OPTIONS must be a 1- or 2-element> stream2 (1,2,3,4, [1,2,3])
0252 %!error <STEPSIZE must be a real scalar != 0> stream2 (1,2,3,4, [1i])
0253 %!error <STEPSIZE must be a real scalar != 0> stream2 (1,2,3,4, [0])
0254 %!error <MAX_VERTICES must be an integer> stream2 (1,2,3,4, [1, 1i])
0255 %!error <MAX_VERTICES must be an integer> stream2 (1,2,3,4, [1, 0])
0256 %!error <matrix dimensions must match> stream2 ([1 1],2,3,4)
0257 %!error <matrix dimensions must match> stream2 (1,[2 2],3,4)
0258 %!error <matrix dimensions must match> stream2 (1,2,[3 3],4)
0259 %!error <matrix dimensions must match> stream2 (1,2,3,[4 4])
0260 %!error <all inputs must be real-valued> stream2 (1i,2,3,4)
0261 %!error <all inputs must be real-valued> stream2 (1,2i,3,4)
0262 %!error <all inputs must be real-valued> stream2 (1,2,3i,4)
0263 %!error <all inputs must be real-valued> stream2 (1,2,3,4i)
0264 %!error <non-monotonically increasing or NaN values found in mesh>
0265 %! stream2 ([2 1], [1 2], [1 1], [2 2], [3 3], [4 4]);
0266 %!error <non-monotonically increasing or NaN values found in mesh>
0267 %! stream2 ([1 NaN], [1 2], [1 1], [2 2], [3 3], [4 4]);
0268 ## FIXME: vectors representing x, y mesh are not accepted.
0269 %#!error <non-monotonically increasing or NaN values found in mesh>
0270 %! stream2 ([1 2], [2 1], [1 1], [2 2], [3 3], [4 4]);
0271 %#!error <non-monotonically increasing or NaN values found in mesh>
0272 %! stream2 ([1 2], [1 NaN], [1 1], [2 2], [3 3], [4 4]);
0273 
0274