function tubeplot3(x,y,z,varargin)
%tubeplot3(x,y,z)
%
%   plot curve as tube
%   x,y,z contain coordinates of points on curve as in plot3(x,y,z)
%
%optional arguments width, color: (can be given in any order)
%
% tubeplot(x,y,z,width,color)
%
%   width: thickness of tube
%   color:   'r', 'g', 'b', 'c', 'm', 'y', 'k'
%        or 'red','green','blue','cyan','magenta','yellow','black'
%        or [r g b]: e.g. [1 0 0] for red
%        or vector of the same size as x,y,z: color according to vector
%
%   default for color is rainbow colored according to z coordinate
%   ( same as tubeplot3(x,y,z,z) )
%
%Examples:
%  t = 0:.1:20;  
%  tubeplot3(cos(t), sin(t), t);  axis equal
%  tubeplot3(cos(t), sin(t), t, 2, 'r'); axis equal

width = [];
color = [];
colorv = [];

for i=1:length(varargin)
  a = varargin{i};
  if isa(a,'double') & length(a)==1
    width = a;
  else
    if isa(a,'double') & length(a)>3
      colorv = a;
    else
      color = a;
    end
  end
end

v{1} = [x(:),y(:),z(:)];
daspect([1 1 1])
if isempty(width)
%   streamtube0(v)
  streamtube0(colorv,v)
else
%   streamtube0(v,width)
  streamtube0(colorv,v,width)
end

view(3);
camlight
lighting gouraud

h = get(gca,'children');
if ~isempty(color)
  set(h(2),'FaceColor',color,'EdgeColor','none');
else
  set(h(2),'FaceColor','interp','EdgeColor','none');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hout=streamtube0(colorv,varargin)

% modified from streamtube
%STREAMTUBE  3D stream tube.

%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 1.6 $  $Date: 2001/04/15 12:04:13 $

[x, y, z, u, v, w, sx, sy, sz, verts, div, width, scale, n] = ...
    parseargs(nargin-1,varargin);

if any(isnan(u(:))) |  any(isnan(v(:))) |  any(isnan(w(:)))
  error('Vector data U,V,W must not contain any NaNs')
end

if isempty(n),     n=20;    end
if isempty(scale), scale=1; 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


dar = [];
f = get(0, 'currentfigure');
if ~isempty(f)
  ax = get(f, 'currentaxes');
  if ~isempty(ax) &  strcmp(get(ax, 'dataaspectratiomode'), 'manual')
    dar = get(ax, 'dataaspectratio');
  end
end

if isempty(dar)
  vv = cat(1,verts{:});
  if isempty(vv)
    dar = [1 1 1];
  else
    dar = max(vv,[],1) - min(vv,[],1);
    dar(dar==0) = 1;
  end
end

dar = dar/max(dar);

% Calculate and scale the widths from the divergence
if isempty(width)
  if ~isempty(u) & isempty(div)
    div = divergence(x,y,z,u,v,w);
  end
  for k = 1:length(verts)
    vv = verts{k};
    if ~isempty(vv)
      if ~isempty(div)
	if ~isempty(x)
	  divi{k} = interp3(x,y,z,div,vv(:,1), vv(:,2), vv(:,3));
	else
	  divi{k} = interp3(div,vv(:,1), vv(:,2), vv(:,3));
	end
      end
      vertminmax(k,:) = [min(vv,[],1) max(vv,[],1)];
    else
      divi{k} = [];
      vertminmax(k,:) = [nan nan nan nan nan nan];
    end
  end
  
  vertmin = min(vertminmax(:,1:3),[],1);
  vertmax = max(vertminmax(:,4:6),[],1);
  maxwidth = .05*max((vertmax-vertmin)./dar);
  if ~isempty(div)
    for k = 1:length(divi)
      if scale==0
	width{k} = divi{k};
      else
	divrange = max(divi{k})-min(divi{k});
	if divrange>0
	  width{k} = (divi{k}-min(divi{k}))/divrange *maxwidth*scale;
	else
	  width{k} = maxwidth*scale;
	end
      end
    end
  else
    width = maxwidth*scale; 
  end
  
else
  if iscell(width)
    for k = 1:length(width)
      width{k} = width{k}*scale;
    end
  else
    width = width*scale;
  end
end



h = [];
for k = 1:length(verts)
  vv = verts{k};
  if ~isempty(vv) & size(vv,1)>1
    if iscell(width)
      [xv yv zv] = tubecoords(vv,width{k},n,dar);
    else
      [xv yv zv] = tubecoords(vv,width,n,dar);
    end
        
    if isempty(colorv)    
      surf(xv,yv,zv);
    else
      cv = repmat(colorv(:),1,n+1);      
      surf(xv,yv,zv,cv);
    end
%    h = [h; surface(xv,yv,zv,cv)];
  end
end

if nargout>0 
  hout = h;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y,z]=tubecoords(verts,width,n,dar)

verts(:,1) = verts(:,1)/dar(1);
verts(:,2) = verts(:,2)/dar(2);
verts(:,3) = verts(:,3)/dar(3);

d1 = diff(verts);
zindex = find(~any(d1,2));
verts(zindex,:) = [];

if size(verts,1)<2
  x = []; y = []; , z = [];
  return;
end
		
d1 = diff(verts);

numverts = size(verts,1);
unitnormals = zeros(numverts,3);

% Radius of the tube.
if length(width)==1
  width = repmat(width, [numverts,1]);
else
  width(zindex) = [];
end
R = width/2;

d1(end+1,:) = d1(end,:);

x10 = verts(:,1)';
x20 = verts(:,2)';
x30 = verts(:,3)';
x11 = d1(:,1)';
x21 = d1(:,2)';
x31 = d1(:,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);
normb = norm(b); if normb~=0, b = b/norm(b); end
%b = b*R(1);

unitnormals(1,:) = b;

for j = 1:numverts-1;

  a = verts(j+1,:)-verts(j,:);
  c = crossSimple(a,b);
  b = crossSimple(c,a);
  normb = norm(b); if normb~=0, b = b/norm(b); end
%  b = b*R(j);
  unitnormals(j+1,:) = b;

end

unitnormal1 = unitnormals(:,1)';
unitnormal2 = unitnormals(:,2)';
unitnormal3 = unitnormals(:,3)';

speed = sqrt(x11.^2 + x21.^2 + x31.^2);

% And the binormal vector ( B = T x N )
binormal1 = (x21.*unitnormal3 - x31.*unitnormal2) ./ speed;
binormal2 = (x31.*unitnormal1 - x11.*unitnormal3) ./ speed;
binormal3 = (x11.*unitnormal2 - x21.*unitnormal1) ./ speed;

% s is the coordinate along the circular cross-sections of the tube:
s = (0:n)';
s = (2*pi/n)*s;

% Finally, the parametric surface.  
% Each of x1, x2, x3 is an (m+1)x(n+1) matrix.
% Rows represent coordinates along the tube.  Columns represent coordinates
% sgcfin each (circular) cross-section of the tube.

xa1 = ones(n+1,1)*x10;
xb1 = (cos(s)*unitnormal1 + sin(s)*binormal1);
xa2 = ones(n+1,1)*x20;
xb2 = (cos(s)*unitnormal2 + sin(s)*binormal2);
xa3 = ones(n+1,1)*x30;
xb3 = (cos(s)*unitnormal3 + sin(s)*binormal3);

R = repmat(R(:)',[n+1,1]);
x1 = xa1 + R.*xb1;
x2 = xa2 + R.*xb2;
x3 = xa3 + R.*xb3;
%x1 = xa1 + xb1;
%x2 = xa2 + xb2;
%x3 = xa3 + xb3;

x = x1' *dar(1);
y = x2' *dar(2);
z = x3' *dar(3);

%nx = unitnormal1;
%ny = unitnormal2;
%nz = unitnormal3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y, z, u, v, w, sx, sy, sz, verts, div, width, scale, n] = parseargs(nin, vargin)

x = [];
y = [];
z = [];
u = [];
v = [];
w = [];
sx = [];
sy = [];
sz = [];
verts = [];
div = [];
width = [];
scale = [];
n = [];

if (nin==6 | nin==7) & ~iscell(vargin{1})
                           % streamtube(u,v,w,sx,sy,sz) or streamtube(u,v,w,sx,sy,sz,scale)
  u = vargin{1};
  v = vargin{2};
  w = vargin{3};
  sx = vargin{4};
  sy = vargin{5};
  sz = vargin{6};
  if nin==7, scale = vargin{7}; end
elseif nin==9 | nin==10    % streamtube(x,y,z,u,v,w,sx,sy,sz) or streamtube(x,y,z,u,v,w,sx,sy,sz,scale)
  x = vargin{1};
  y = vargin{2};
  z = vargin{3};
  u = vargin{4};
  v = vargin{5};
  w = vargin{6};
  sx = vargin{7};
  sy = vargin{8};
  sz = vargin{9};
  if nin==10, scale = vargin{10}; end
elseif (nin==2 | nin==3 ) & ...
	    (length(size(vargin{2}))==3 | ...
	    iscell(vargin{2}) | ...
	    prod(size(vargin{2}))==1)
                         % streamtube(verts,divergence) or streamtube(verts,divergence,scale)
                         % streamtube(verts,width) or streamtube(verts,width,scale)
  verts = vargin{1};
  if length(size(vargin{2}))==3
    div = vargin{2};
  else
    width = vargin{2};
  end    
  if nin==3, scale = vargin{3}; end
elseif nin==5 | nin==6    % streamtube(verts,x,y,z,divergence) or streamtube(verts,x,y,z,divergence,scale)
  verts = vargin{1};
  x = vargin{2};
  y = vargin{3};
  z = vargin{4};
  div = vargin{5};
  if nin==6, scale = vargin{6}; end
elseif nin==1 | nin==2    % streamtube(verts)  or streamtube(verts,scale) 
  verts = vargin{1};
  if nin==2, scale = vargin{2}; end
else
  error('Wrong number of input arguments.'); 
end

if ~isempty(sx)
  sx = sx(:); 
  sy = sy(:); 
  sz = sz(:); 
end

if ~isempty(scale)
  if length(scale) == 2
    n = scale(2);
    scale = scale(1);
  else
    error('[SCALE N] must be a 2 element vector.'); 
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [msg,nx,ny,nz] = xyzuvwcheck(x,y,z,u,v,w)
%XYZUVWCHECK  Check arguments to 3D vector data routines.
%   [MSG,X,Y,Z] = XYZCHK(X,Y,Z,U,V,W) checks the input aguments
%   and returns either an error message in MSG or valid X,Y,Z.

%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 1.5 $  $Date: 2001/04/15 12:04:02 $

msg = '';
nx = x;
ny = y;
nz = z;

sz = size(u);
if ~isequal(size(v), sz) | ~isequal(size(w), sz)
  msg='U,V,W must all be the same size.';
  return
end

if ndims(u)~=3
  msg='U,V,W must all be a 3D arrays.';
  return
end
if min(sz)<2
  msg = 'U,V,W must all be size 2x2x2 or greater.'; 
  return
end

nonempty = ~[isempty(x) isempty(y) isempty(z)];
if any(nonempty) & ~all(nonempty)
  msg = 'X,Y,Z must all be empty or all non-empty.';
  return;
end

if ~isempty(nx) & ~isequal(size(nx), sz)
  nx = nx(:);
  if length(nx)~=sz(2)
    msg='The size of X must match the size of U or the number of columns of U.';
    return
  else
    nx = repmat(nx',[sz(1) 1 sz(3)]);
  end
end

if ~isempty(ny) & ~isequal(size(ny), sz)
  ny = ny(:);
  if length(ny)~=sz(1)
    msg='The size of Y must match the size of U or the number of rows of U.';
    return
  else
    ny = repmat(ny,[1 sz(2) sz(3)]);
  end
end

if ~isempty(nz) & ~isequal(size(nz), sz)
  nz = nz(:);
  if length(nz)~=sz(3)
    msg='The size of Z must match the size of U or the number of pages of U.';
    return
  else
    nz = repmat(reshape(nz,[1 1 length(nz)]),[sz(1) sz(2) 1]);
  end
end