%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Project     : PON_HD
% Module      : ScatterColormap
% Description : constellation scatterplot with colormap
%
% References  :
%
%Version|Date    |author|Description
%-------|--------|------|-------------------------------------------------------
%01.01  |20200923|L,YM  |Initial version
%-------|--------|------|-------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ScatterColormap(ax,x,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% coeffients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%constraint max DataLength to MaxLength;
MaxLength  = 1e5; %1e4 for QuickCheck, 1e5 for accurate check.
DataLength = min(length(x), MaxLength);

x = x(end - DataLength + 1: end);
y = y(end - DataLength + 1: end);

%set radius
DivideNum = 30;
DataMax = max(max(abs(x)),max(abs(y)));
r = DataMax / DivideNum; %only use positive data.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Scatter Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% calculate density
%divide figure into small grids for higher plot speed.

Density = zeros(DataLength, 1);
for i1 = -DivideNum : DivideNum %x
  for i2 = -DivideNum : DivideNum %y
    LargeIndex = find(  x >  (i1 - 3/2) * r & x <= (i1 + 3/2) * r ...
                      & y >  (i2 - 3/2) * r & y <= (i2 + 3/2) * r);
    SmallIndex = find(  x >  (i1 - 1/2) * r & x <= (i1 + 1/2) * r ...
                      & y >  (i2 - 1/2) * r & y <= (i2 + 1/2) * r);

    if ~isempty(SmallIndex)
      for i3 = SmallIndex'
        Density(i3,1) = sum( (x(LargeIndex) - x(i3)).^2 ...
                            +(y(LargeIndex) - y(i3)).^2 < r.^2);
  %       Density(i3,1) = sum( abs(x(LargeIndex) - x(i3)) < r ...
  %                           &abs(y(LargeIndex) - y(i3)) < r);
      end
    end
    
  end
end

Density = Density ./ max(Density); %normalize
% Density = sqrt(Density); %re-mapping

%% plot

map = colormap;
MapLength = size(map,1);
Density = ceil(Density .* MapLength); %fit density to colormap

for i1 = 1 : MapLength
  if any(Density == i1) 
    plot(ax, x(Density == i1), y(Density == i1), ...
         'LineStyle','none', 'Color',map(i1,:), 'Marker','.'); 
    hold on;
  end
end

end%function