clc;clear all;close all;
addpath('/home/jhluna/tools/ncl_colormap/')
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/'
pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Check_UV/'
% load([pathi,'Jul_0727_5D.mat'])
% load([pathi,'Aug_0829_5D.mat'])
load([pathi,'Sep_0903_5D.mat'])


% xx(end) = [];
% yy(end) = [];
%====== [pos(step) =step+1 --> NAN] ====%
for IDX=1:size(xx,2)      %realsed time index
for ii=1:size(xx{1},1)    %particles  
for k=1:size(xx{IDX},2)-1  
 if xx{IDX}(ii,k)==xx{IDX}(ii,k+1)
    xx{IDX}(ii,k)=nan;
 elseif yy{IDX}(ii,k)==yy{IDX}(ii,k+1)  
    yy{IDX}(ii,k)=nan; 
 end
end
end
end

%======= <STEP1>find out the final position =======%
final_x = cell(size(xx,2), 1);
final_y = cell(size(xx,2), 1);
for IDX=1:size(xx,2)
    final_x{IDX} = NaN(size(xx{1},1), 1);
    final_y{IDX} = NaN(size(xx{1},1), 1);    
 for i = 1:size(xx{1},1)
    valid_indicesX = find(~isnan(xx{IDX}(i, :))); 
    valid_indicesY = find(~isnan(yy{IDX}(i, :))); 
    if ~isempty(valid_indicesX)
        final_x{IDX}(i) = valid_indicesX(end); 
    end    
    if ~isempty(valid_indicesY)
        final_y{IDX}(i) = valid_indicesY(end);
    end
 end
end

%======= <STEP2>find out the final position =======%
particle_data = cell(size(xx,2), 1);
for IDX=1:size(xx,2)
particle_data{IDX} = []; % [k, i(xx), j(yy)]

    for k=1:length(final_x{IDX})
         if ~isnan(final_x{IDX}(k)) 
            i_xx = xx{IDX}(k, final_x{IDX}(k));
         end
         if ~isnan(final_y{IDX}(k))
            j_yy = yy{IDX}(k, final_y{IDX}(k));    
         end
          particle_data{IDX} = [particle_data{IDX}; k, i_xx, j_yy];    
    end
end


%======= Count the repeat(i,j) numbers =======%
all_lon_lat = [];
lon_lat={};unique_lon_lat={};
idx={};Count={};
for IDX = 1:size(particle_data,1)
    current_data = particle_data{IDX};
    
    for k = 1:size(current_data, 1)      
%---- fix the 0. number[7] ------%
        lon_lat{IDX}(k,:) = round(current_data(k, 2:3), 5); 
%--------------------------------------%        
%       all_lon_lat = [all_lon_lat; lon_lat]; 

% Determine unique latitude and longitude points and their number of occurrences
[unique_lon_lat{IDX}, ~, idx{IDX}] = unique(lon_lat{IDX}, 'rows', 'stable');
counts{IDX} = accumarray(idx{IDX}, 1);

% Result: Count(i, j, number)
Count{IDX} = [unique_lon_lat{IDX}, counts{IDX}];
    end
end
% [unique_lon_lat, ~, idx] = unique(all_lon_lat, 'rows', 'stable');
% counts = accumarray(idx, 1);

% % Result: Count(i, j, number)
% Count = [unique_lon_lat, counts];

%===== save to mat =======%
patho=pathi;
fn = [patho 'POS_0903_5D.mat'];
save(fn,'Count','lon_lat')

%%
%====== figure ======%
% plot(Count(:,1), Count(:,2),'.','Color', 'r','Markersize',10);hold on

% bin_edges = 0:500:8500; 
% num_bins = length(bin_edges) - 1; 
% 
% [~, bin_idx] = histc(Count(:, 3), bin_edges);
% % colors = lines(num_bins); 
% colors = colormap(nclcmap('GMT_seis',num_bins)); 
% 
% figure;
% hold on;
% for i = 1:num_bins
%     idx = (bin_idx == i);
%     scatter(Count(idx, 1), Count(idx, 2), 50, colors(i, :), 'filled'); 
% end
% axis square;
% colormap(colors);
% cbar = colorbar;
% cbar.Ticks = linspace(0, 1, num_bins); 
% cbar.TickLabels = arrayfun(@(x) sprintf('%d-%d', bin_edges(x),...
%          bin_edges(x + 1) - 1), 1:num_bins, 'UniformOutput', false);
% 
















