function quiverC3D(x,y,z,u,v,w,varargin)
%quiverC3D creates a 3D quiver plot and adds a color coding. The color coding is
%given by the magnitudes of the component vectors. Large values result in colors 
%from the upper end of the used colormap. 
% 
%   INPUT:
%       x - array, x components of initial points
%       y - array, y components of initial points
%       z - array, z components of initial points
%       u - array, x components of arrows
%       v - array, y components of arrows
%       w - array, z components of arrows
%       scale - number, if > 0, automatic scaling is used to avoid overlap
%           of adjacent arrows. If scale = 0, the scaling is disabled.
%       MaxNumArrows - a positive integer (non-integer should work as well)
%           number limiting the maximum number of plotted arrows. Since vectors
%           of length zero aren't plotted and the matrices are sampled
%           uniformly, it's possible that fewer arrows will be visible in the
%           end. If MaxNumArrows is not given its value is set to inf.
% 
%   OUTPUT:
%       none
% 
% WARNING!: This function might not create arrows with the same length as
%   would be obtained using quiver3(x, y, z, u, v, w). I do not know
%   whether the scaling method is the same is in quiver3.
% 
% 
% --------------------------------------------------------------------------------------
% 
%   EXAMPLE:
%       load wind
%       quiverC3D(x,y,z,u,v,w,1,1000);
%       view(-45,20)
%       set(gca, 'Color', 'black');
%       axis tight
%       axis equal
%   
% --------------------------------------------------------------------------------------
% 
%   See also: QUIVER3, LINESPEC, COLORMAP.
% 

%% parse arguments

p = inputParser;
addRequired(p, 'x', @isnumeric);
addRequired(p, 'y', @isnumeric);
addRequired(p, 'z', @isnumeric);
addRequired(p, 'u', @isnumeric);
addRequired(p, 'v', @isnumeric);
addRequired(p, 'w', @isnumeric);
addOptional(p, 'scale', 1, @isscalar);
addOptional(p, 'CLim', [nan nan], @isnumeric);
addOptional(p, 'MaxNumArrows', inf, @isscalar);
addParameter(p, 'Normalize', false);
addParameter(p, 'LineWidth', 0.5, @isscalar);
addParameter(p, 'MaxHeadSize', 1, @isscalar);
addParameter(p, 'Colorbar', false, @islogical);
addParameter(p, 'Parent', gca(), @isgraphics);

parse(p, x, y, z, u, v, w, varargin{:});

scale = 0;%p.Results.scale;
maxnumarrows = p.Results.MaxNumArrows;
normalize = p.Results.Normalize;
lw = p.Results.LineWidth;
hs = p.Results.MaxHeadSize;
ax = p.Results.Parent;
c_lim = p.Results.CLim;

assert(scale >= 0);

if ~isequal(size(x),size(y),size(z),size(u),size(v),size(w))
    error('X,Y,Z,U,V,W have to be arrays of the same size.');
end

%% initialization
if numel(u) > maxnumarrows
    if isvector(u)
        N = ceil(numel(u)/maxnumarrows);
    else
        N = ceil(nthroot(numel(u)/maxnumarrows, ndims(u)));
    end
    
    if isvector(x)
        idx = {1:N:length(x)};
    else
        idx = arrayfun(@(x) 1:N:x, size(x), 'UniformOutput', false);
    end
    
    x = x(idx{:});
    y = y(idx{:});
    z = z(idx{:});
    u = u(idx{:});
    v = v(idx{:});
    w = w(idx{:});
end


%% colormap definition
C = colormap(ax);
ncolors = size(C, 1);
I = sqrt(u.^2 + v.^2 + w.^2);

if isnan(c_lim)
	if verLessThan('matlab','9.6')
		Imax = max(I(:), [], 'omitnan');
		Imin = min(I(:), [], 'omitnan');
	else
		Imax = max(I, [], 'all', 'omitnan');
		Imin = min(I, [], 'all', 'omitnan');
	end
else
	Imax=c_lim(2);
	Imin=c_lim(1);
end
% assume that the colormap matrix contains the colors in its rows
Ic = round((I-Imin)/(Imax-Imin)*ncolors);
Ic(Ic == 0) = 1;

%  scale to avoid overlaps
if scale > 0
    [~,dxi] = uniquetol(x);
    dx = diff(x(sort(dxi)));
    [~,dyi] = uniquetol(y);
    dy = diff(y(sort(dyi)));
    [~,dzi] = uniquetol(z);
    dz = diff(z(sort(dzi)));
    
    % avg difference
    dr_mean = mean(abs([dx; dy; dz]));
    
    lenscale = scale*dr_mean/Imax;
else
    lenscale = 1;
end

if normalize
    lenscale = lenscale./I*Imax;
end

u = u.*lenscale;
v = v.*lenscale;
w = w.*lenscale;


%% plotting
ishold_flag = ishold(ax);
hold(ax, 'on');

if numel(u) > ncolors
    % let's take an intelligent approach: group all values which belong to
    % the same color value
    for k = 1:ncolors
        mask = (Ic == k);    
        quiver3(ax, x(mask), y(mask), z(mask), u(mask), v(mask), w(mask), 0, ...
            'Color', C(k, :), 'LineWidth', lw, 'MaxHeadSize', hs);
    end
else
    % if there are so few values, it will be more efficient to plot each
    % arrow individually
    % linear indexing!
    for k = 1:numel(u)
        quiver3(ax, x(k), y(k), z(k), u(k), v(k), w(k), 0, 'Color', C(Ic(k),:), ...
            'LineWidth', lw, 'MaxHeadSize', hs);
    end
end

if ~ishold_flag
    hold(ax, 'off');
end

if p.Results.Colorbar
    cb = colorbar(ax);
    tlabels = cb.TickLabels;
    tlabels = sprintfc('%1.3g ',linspace(Imin, Imax, length(tlabels)));
    cb.TickLabels = tlabels;
end

end


