function [ worked ] = PlotColorPoints( x, y, z, cmap, n_bands )

%Plots the data in different colors - similiar to scatter, but much
%faster.
%   Inputs:  x, y, z, colormap2, n
%      x, y, z should be obvious
%      cmap is colormap to use for point colors
%      n_bands is the number of color bands to use
%
%The one trick to this is that rather than a proportional spacing of the
%color bands on the extent of z, instead this uses a histogram of z to
%allocate bands of equal number of points per band...

% Make sure no NaNs...
ok = find(~isnan(z) & ~isnan(y) & ~isnan(x));
x = x(ok);
y = y(ok);
z = z(ok);

% Band the data by color...
[sorted_z, index] = sort(z);  % First order the indices by increasing z
n_data = length(index);       % Length of data array
ci(1:n_bands, 1) = (floor(1+n_data*(0:n_bands-1)/n_bands)); % Now find lower index in each of n color band
ci(1:n_bands, 2) = (floor(n_data*(1:n_bands)/n_bands));     % Find upper index in each of n color bands
nn = length(cmap);  % Find the length of the user-provided color map

% Set up for calibrating color bar...
ii1 = round(1);
ii2 = round(n_data);
z1 = sorted_z(ii1); % Lower z limit
z2 = sorted_z(ii2); % Upper z limit
zbar(1:(n_bands+1)) = sorted_z(1+floor((n_data-1)*(0:n_bands)/n_bands));

for i = n_bands:-1:1
    cmap_index = (1 + floor((nn-1)*(i-1)/(n_bands-1)));
    ok = index(ci(i,1):ci(i,2)); % Indices corresponding to the range of data we are plotting
    plot(x(ok), y(ok), '.', 'color', cmap(cmap_index, :));
    
    cbar_index = floor(1+99*(zbar(i)-z1)/(z2-z1)):floor(1+99*(zbar(i+1)-z1)/(z2-z1));
    cc(cbar_index,1) = cmap(cmap_index, 1);
    cc(cbar_index,2) = cmap(cmap_index, 2);
    cc(cbar_index,3) = cmap(cmap_index, 3);
    hold on
end

colormap(cc);
colorbar;
caxis([z1 z2]);
hold off


