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

%Plots the data in different colors - similiar to scatter, but much
%faster.
%   Inputs:  x, y, z, colormap2, n_bands, cscale
%      x, y, z should be obvious
%      cmap is colormap to use for point colors
%      n_bands is the number of color bands to use
%      cscale is 'linear', 'proportional', or 'log' for color scaling
%
%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...
%
% Useful routine for collecting bad data segments.
% for i=1:11, 
%    bb(i) = [eval(['b' num2str(i) '1.Position(1)']):eval(['b' num2str(i) '1.Position(1)'])]; 
% end;

% Set up defaults;
cscale = 'linear';
cmap = jet;
nbands = 10;

% Deal with the variable arguements
for i = 4:nargin
    if iscellstr(varargin(i-3))
        if(strcmp(varargin(i-3),'linear') | strcmp(varargin(i-3),'proportional')| strcmp(varargin(i-3),'log'))
            cscale = varargin{i-3}
        else
            cmap = eval(varargin{i-3});
        end;
    else
        n_bands = varargin{i-3};
    end;
end;

% 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
nn = length(cmap);  % Find the length of the user-provided color map

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

% Select range of variable to plot in each color band depending on user
% defined scale.
switch cscale
    case 'linear'
        z_bands = (0:n_bands)*(z2-z1)/n_bands + z1;
        for i_bands = 1:n_bands
            ci(i_bands, 1) = find(sorted_z > z_bands(i_bands), 1,'first'); % Now find lower index in each of n color band
            ci(i_bands, 2) = find(sorted_z < z_bands(i_bands+1), 1,'last'); % Find upper index in each of n color bands
        end;
        ci(1,1) = 1;
        ci(n_bands, 2) = n_data;
    case 'proportional'
        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
        z_bands(1:n_bands) = sorted_z(ci(1:n_bands,1));
        z_bands(n_bands+1) = sorted_z(ci(n_bands, 2));
    case 'log'
        z_bands = (0:n_bands)*(log10(z2)-log10(z1))/n_bands + log10(z1);
        z_bands = 10.^z_bands;
        
        for i_bands = 1:n_bands
            ci(i_bands, 1) = find(sorted_z > z_bands(i_bands), 1,'first'); % Now find lower index in each of n color band
            ci(i_bands, 2) = find(sorted_z < z_bands(i_bands+1), 1,'last'); % Find upper index in each of n color bands
        end;
        ci(1,1) = 1;
        ci(n_bands, 2) = n_data;

end;

% Set the color for each band...
cmap_index = round(1+(nn-1)*((0:n_bands-1)/(n_bands-1)));

for i = n_bands:-1:1 %1:n_bands
    % Pick the color, select the subset of data, and the plot the data
    % points in that color...
    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(i), :), 'MarkerSize', 5);
    
    % Fill in the colormap that will be used in the colorbar
    cbar_index = floor(1+99*(z_bands(i)-z1)/(z2-z1)):floor(1+99*(z_bands(i+1)-z1)/(z2-z1));
    if (cbar_index > 0 & cbar_index <101)
        cc(cbar_index,1) = cmap(cmap_index(i), 1);
        cc(cbar_index,2) = cmap(cmap_index(i), 2);
        cc(cbar_index,3) = cmap(cmap_index(i), 3);
    end;
    hold on
end

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


