function hout=streamribbon(varargin)
    %STREAMRIBBON  3D stream ribbon.
    %   STREAMRIBBON(X,Y,Z,U,V,W,STARTX,STARTY,STARTZ) draws stream
    %   ribbons from vector data U,V,W. The arrays X,Y,Z define the
    %   coordinates for U,V,W and must be monotonic and 3D plaid (as if
    %   produced by MESHGRID). STARTX, STARTY, and STARTZ define the
    %   starting positions of the streamlines at the center of the
    %   ribbons.  The twist of the ribbons is proportional to the
    %   curl of the vector field. The width of the ribbons is
    %   calculated  automatically.
    %
    %   STREAMRIBBON(U,V,W,STARTX,STARTY,STARTZ) assumes
    %      [X Y Z] = meshgrid(1:N, 1:M, 1:P) where [M,N,P]=SIZE(U).
    %
    %   STREAMRIBBON(VERTICES,X,Y,Z,CAV,SPEED) assumes precomputed
    %   streamline vertices, curl angular velocity, and flow
    %   speed. VERTICES is cell array of streamline vertices (as if
    %   produced by stream3). X,Y,Z, CAV and SPEED are 3D arrays.
    %
    %   STREAMRIBBON(VERTICES,CAV,SPEED)  assumes
    %      [X Y Z] = meshgrid(1:N, 1:M, 1:P) where [M,N,P] = SIZE(CAV).
    %
    %   STREAMRIBBON(VERTICES,TWISTANGLE) uses the cell array of
    %   vectors TWISTANGLE for the twist of the ribbons (in
    %   radians). The size of each corresponding element of VERTICES
    %   and TWISTANGLE must be equal.
    %
    %   STREAMRIBBON(...,WIDTH) sets the width of the ribbons to be
    %   WIDTH.
    %
    %   STREAMRIBBON(AX,...) plots into AX instead of GCA.
    %
    %   H = STREAMRIBBON(...)  returns a vector of handles (one per
    %   start point) to SURFACE objects.
    %
    %   Example 1:
    %      load wind
    %      [sx sy sz] = meshgrid(80, 20:10:50, 0:5:15);
    %      daspect([1 1 1])
    %      h=streamribbon(x,y,z,u,v,w,sx,sy,sz);
    %      axis tight
    %      shading interp;
    %      view(3);
    %      camlight; lighting gouraud
    %
    %   Example 2:
    %      load wind
    %      [sx sy sz] = meshgrid(80, 20:10:50, 0:5:15);
    %      daspect([1 1 1])
    %      verts = stream3(x,y,z,u,v,w,sx,sy,sz);
    %      cav = curl(x,y,z,u,v,w);
    %      spd = sqrt(u.^2 + v.^2 + w.^2);
    %      h=streamribbon(verts,x,y,z,cav,spd);
    %      axis tight
    %      shading interp;
    %      view(3);
    %      camlight; lighting gouraud
    %
    %   Example 3:
    %      t = 0:.15:15;
    %      verts = {[cos(t)' sin(t)' (t/3)']};
    %      twistangle = {cos(t)'};
    %      daspect([1 1 1])
    %      h=streamribbon(verts,twistangle);
    %      axis tight
    %      shading interp;
    %      view(3);
    %      camlight; lighting gouraud
    %
    %   Example 4:
    %     xmin = -7; xmax = 7;
    %     ymin = -7; ymax = 7;
    %     zmin = -7; zmax = 7;
    %     x = linspace(xmin,xmax,30);
    %     y = linspace(ymin,ymax,20);
    %     z = linspace(zmin,zmax,20);
    %     [x y z] = meshgrid(x,y,z);
    %     u = y; v = -x; w = 0*x+1;
    %
    %     daspect([1 1 1]);
    %     [cx cy cz] = meshgrid(linspace(xmin,xmax,30),linspace(ymin,ymax,30),[-3 4]);
    %     h2=coneplot(x,y,z,u,v,w,cx,cy,cz, 'q');
    %     h2.Color = 'k';
    %
    %     [sx sy sz] = meshgrid([-1 0 1],[-1 0 1],-6);
    %     p = streamribbon(x,y,z,u,v,w,sx,sy,sz);
    %     [sx sy sz] = meshgrid([1:6],[0],-6);
    %     p2 = streamribbon(x,y,z,u,v,w,sx,sy,sz);
    %     shading interp
    %
    %     view(-30,10) ; axis off tight
    %     camproj p; camva(66); camlookat; camdolly(0,0,.5,'f')
    %     camlight
    %
    %   See also CURL, STREAMTUBE, STREAMLINE, STREAM3.
    
    %   Copyright 1984-2018 The MathWorks, Inc.
    
    [cax,args,nargs] = axescheck(varargin{:});
    [x, y, z, u, v, w, sx, sy, sz, verts, cav, spd, twistangle, width] = ...
        parseargs(nargs,args);
    
    if any(isnan(u(:))) ||  any(isnan(v(:))) ||  any(isnan(w(:)))
        error(message('MATLAB:streamribbon:NonNanData'))
    end
    
    % Calculate the vertices
    if isempty(verts)
        [msg, x, y, z] = xyzuvwcheck(x,y,z,u,v,w);  error(msg) 
        verts = stream3(x,y,z,u,v,w,sx,sy,sz);
    end
    
    if isempty(cax)
        f = get(0, 'currentfigure');
        if ~isempty(f)
            cax = get(f, 'currentaxes');
        end
    end
    
    % Get the axes if none has been specified yet
    if isempty(cax)
        cax = newplot;
    end
    
    % Calculate the twistangles
    if isempty(twistangle)
        if isempty(cav)
            cav = curl(x,y,z,u,v,w);
            spd = sqrt(u.^2 + v.^2 + w.^2 );
        end
        
        cavi = cell(1,length(verts));
        spdi = cell(1,length(verts));
        for k = 1:length(verts)
            vv = verts{k};
            if ~isempty(vv)
                if ~isempty(x)
                    cavi{k} = interp3(x,y,z,cav,vv(:,1), vv(:,2), vv(:,3));
                    spdi{k} = interp3(x,y,z,spd,vv(:,1), vv(:,2), vv(:,3));
                else
                    cavi{k} = interp3(cav,vv(:,1), vv(:,2), vv(:,3));
                    spdi{k} = interp3(spd,vv(:,1), vv(:,2), vv(:,3));
                end
            end
        end
        
    end
    
    % Calculate a width if none was passed in.
    if isempty(width)
        for k = 1:length(verts)
            vv = verts{k};
            if ~isempty(vv)
                vertminmax(k,:) = [min(vv,[],1) max(vv,[],1)];
            else
                vertminmax(k,:) = [nan nan nan nan nan nan];
            end
        end
        
        vertmin = min(vertminmax(:,1:3),[],1);
        vertmax = max(vertminmax(:,4:6),[],1);
        width = .05*max(vertmax-vertmin);
    end
    
    
    h = [];
    for k = 1:length(verts)
        vv = verts{k};
        if ~isempty(vv) && size(vv,1)>1
            if isempty(twistangle)
                %fvc = ribboncoords(vv,cavi{k},spdi{k},width);
                %h = [h; patch(fvc)];
                [xyzc] = ribboncoords(vv,cavi{k},spdi{k},width);
                xyzc.Parent = cax;
                h = [h; surface(xyzc)];
            else
                %fvc = ribboncoords(vv,twistangle{k},width);
                %h = [h; patch(fvc)];
                [xyzc] = ribboncoords(vv,twistangle{k},width);
                xyzc.Parent = cax;
                h = [h; surface(xyzc)];
            end
        end
    end
    
    % Disable data tips
    for i = 1:numel(h)
        set(hggetbehavior(cax, 'DataCursor'), 'Enable', false);
        setInteractionHint(h(i), 'DataCursor', false);
    end

    
    if nargout>0
        hout = h;
    end
    
    % Register handles with MATLAB code generator
    if ~isempty(h)
        if ~isdeployed
            makemcode('RegisterHandle',h,'IgnoreHandle',h(1),'FunctionName','streamribbon');
        end
    end
end

function ret=ribboncoords(verts,arg1,arg2,arg3)
    
    if nargin==4
        cav = arg1;
        spd = arg2;
        twistangle = [];
        width = arg3;
    else
        cav = [];
        spd = [];
        twistangle = arg1;
        width = arg2;
    end
    
    R = width/2;
    
    numverts = size(verts,1);
    srv = zeros(numverts*2,3);
    
    
    a = verts(2,:) - verts(1,:);
    b = [0 0 1];
    c = crossSimple(a,b);
    if ~any(c)
        b = [1 0 0];
        c = crossSimple(a,b);
    end
    b = crossSimple(c,a);
    if ~isempty(twistangle)
        b = rotategeom(b, a,  twistangle(1));
    end
    normb = norm(b); if normb~=0, b = b/norm(b); end
    b = b*R;
    
    srv(1,:) = verts(1,:)+b;
    srv(2,:) = verts(1,:)-b;
    
    
    if isempty(twistangle)
        
        dis = sqrt(sum(diff(verts).^2,2)); %distance
        q = zeros(numverts,1); % the currenty angle
        spd = (spd(1:numverts-1)+spd(2:numverts))/2;
        cav = (cav(1:numverts-1)+cav(2:numverts))/2;
        for j = 1:numverts-1
            t = dis(j)/spd(j);
            
            q(j+1) = q(j)+cav(j)*t;
            
            a = verts(j+1,:)-verts(j,:);
            nb = rotategeom(b, a,  q(j+1)-q(j));
            c = crossSimple(a,nb);
            b = crossSimple(c,a);
            normb = norm(b); if normb~=0, b = b/norm(b); end
            b = b*R;
            
            srv(2*j+1,:) = verts(j+1,:)+b;
            srv(2*j+2,:) = verts(j+1,:)-b;
        end
        
    else
        
        q = cumsum(twistangle(:));
        for j = 1:numverts-1
            a = verts(j+1,:)-verts(j,:);
            nb = rotategeom(b, a, twistangle(j+1));
            c = crossSimple(a,nb);
            b = crossSimple(c,a);
            normb = norm(b); if normb~=0, b = b/norm(b); end
            b = b*R;
            
            srv(2*j+1,:) = verts(j+1,:)+b;
            srv(2*j+2,:) = verts(j+1,:)-b;
        end
        
    end
    
    %fvc = [q';q'];
    %f = (1:2:size(srv,1)-2)';
    %f = [f f+1 f+3 f+2];
    
    %ret.faces = f;
    %ret.vertices = srv;
    %ret.facevertexcdata = fvc(:);
    
    ret.XData = reshape(srv(:,1)',[2,numverts])';
    ret.ydata = reshape(srv(:,2)',[2,numverts])';
    ret.zdata = reshape(srv(:,3)',[2,numverts])';
    ret.cdata = [q q];
end

% simple cross product
function c=crossSimple(a,b)
    c(1) = b(3)*a(2) - b(2)*a(3);
    c(2) = b(1)*a(3) - b(3)*a(1);
    c(3) = b(2)*a(1) - b(1)*a(2);
end

function vout=rotategeom(v,azel,alpha)
    nrmazel = norm(azel);
    if nrmazel==0 || alpha==0
        vout = v;
    else
        u = azel(:)/nrmazel;
        cosa = cos(alpha);
        sina = sin(alpha);
        vera = 1 - cosa;
        x = u(1);
        y = u(2);
        z = u(3);
        rot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; ...
            x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; ...
            x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';
        
        x = v(:,1);
        y = v(:,2);
        z = v(:,3);
        
        [m,n] = size(x);
        newxyz = [x(:), y(:), z(:)];
        newxyz = newxyz*rot;
        newx = reshape(newxyz(:,1),m,n);
        newy = reshape(newxyz(:,2),m,n);
        newz = reshape(newxyz(:,3),m,n);
        
        vout = [newx newy newz];
    end
end

function [x, y, z, u, v, w, sx, sy, sz, verts, cav, spd, twistangle, width] = parseargs(nin, vargin)
    
    [x, y, z, u, v, w, sx, sy, sz, verts, cav, spd, twistangle, width] = deal([]);
    
    if (nin==6 || nin==7) && ~iscell(vargin{1})
        % streamribbon(u,v,w,sx,sy,sz) or streamribbon(u,v,w,sx,sy,sz,width)
        [u, v, w, sx, sy, sz] = deal(vargin{1:6});
        if nin==7, width = vargin{7}; end
    elseif nin==9 || nin==10    % streamribbon(x,y,z,u,v,w,sx,sy,sz) or streamribbon(x,y,z,u,v,w,sx,sy,sz,width)
        [x, y, z, u, v, w, sx, sy, sz] = deal(vargin{1:9});
        if nin==10, width = vargin{10}; end
    elseif (nin==2 || nin==3 || nin==4)
        % streamribbon(verts,cav,spd) or streamribbon(verts,cav,spd,width)
        % streamribbon(verts,twistangle) or streamribbon(verts,twistangle,width)
        verts = vargin{1};
        if length(size(vargin{2}))==3
            cav = vargin{2};
            spd = vargin{3};
            if nin==4, width = vargin{4}; end
        else
            twistangle = vargin{2};
            if nin==3, width = vargin{3}; end
        end
    elseif nin==6 || nin==7    % streamribbon(verts,x,y,z,cav,spd) or streamribbon(verts,x,y,z,cav,spd,width)
        [verts, x, y, z, cav, spd] = deal(vargin{1:6});
        if nin==7, width = vargin{7}; end
    else
        error(message('MATLAB:streamribbon:WrongNumberOfInputs'));
    end
    
    if ~isempty(sx)
        sx = sx(:);
        sy = sy(:);
        sz = sz(:);
    end
end

