clc;clear all;close all
addpath('/home/jhluna/tools/ncl_colormap/')
addpath('/data3/yiwen/matlab_toolbox/tight_subplot/')
addpath('/data3/yiwen/matlab_coawst/particle_tracking/')
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/'
% patht='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/Tide/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/'
% patht='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/Tide/'
%------%
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/Tide/_3hrs/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/Tide/+3hrs/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/Tide/_3hrs/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/Tide/+3hrs/'

pathi2='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/'
pathi3='/data3/yiwen/matlab_coawst/particle_tracking/Final_report/'
%====== Spring & Neap ======%
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/Tide/Spring/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/Tide/Neap/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/Tide/Spring/'
pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Aug_5days/Tide/Neap/'

%===== Reading files[for every day] =====%
% load([pathi,'POS_0712_0731_5D.mat']);load([pathi,'Tide_0712_0731.mat'])
load([pathi,'POS_0804_0823_5D.mat']);load([pathi,'Tide_0804_0823.mat'])


% Jul=9:28; %old case
% Aug=4:23;
%---- Check UV ----%
Mon=1:20; %[7/12 00:00 ~ 7/31 21:00 | 8/04 00:00 ~ 8/23 21:00]

allLONLAT_LT=all_lon_lat(Mon);
allLONLAT_HT=all_lon_lat2(Mon);
allLONLAT={allLONLAT_LT; allLONLAT_HT};

% batchSize = 1; %for running 6-days
batchSize =20;   %Total days
numBatches = floor(length(allLONLAT) / batchSize); 
mergedallLonLat = {};

for i=1:size(allLONLAT,1)
    mergedallLonLat{i} = cell2mat(allLONLAT{i}');
end
%======= Count the repeat(i,j) numbers =======%
all_lon_lat = {};
unique_lon_lat={};
idx={};Count_6D={};
sumValue={0 0};

for IDX = 1:size(mergedallLonLat,2)
%--------------------------------------%
all_lon_lat{IDX} = round(mergedallLonLat{IDX}(:,1:2),9);
% Determine unique latitude and longitude points and their number of occurrences
[unique_lon_lat{IDX}, ~, idx{IDX}] = unique(all_lon_lat{IDX}, 'rows', 'stable');
counts{IDX} = accumarray(idx{IDX}, 1);

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

%====== Count the amount of garbage ======% 
    for k=1:size(Count_6D{IDX},1)
     if (Count_6D{IDX}(k,2)<25.124)...
     || (Count_6D{IDX}(k,2)<25.124 && Count_6D{IDX}(k,1)>121.9155 )

        sumValue{IDX} = sumValue{IDX} + Count_6D{IDX}(k,3); 
     end     
    end

end

%======== Velocity =========%
% load([pathi2,'202407.mat']);
load([pathi2,'202408_0906.mat']);

load ([pathi3,'blue.mat']);
%======== Time index ==========%
startTime = datetime(2024, 8 ,04, 0, 0, 0); %7/12 | 8/04  
endTime = datetime(2024, 9, 05, 23, 54, 0);
timeStep = hours(3);                        %particles raleased timestep
timeSeries = startTime:timeStep:endTime;
%------ check val ------%
% T=sort(LT_idx);disp(timeSeries(T(7)))
% T=sort(HT_idx);disp(timeSeries(T(7)))
disp(nansum(Count_6D{1}(:,3)))
disp(nansum(Count_6D{2}(:,3)))

% Compute V2 matrix [LT HT]
% start = [4142 4262];  %Jul
% start = [1622+60 1422+60];  %Aug
% start = [7382 7202];        %Jul[Spring]   
% start = [4382 4262];        %Jul[Neap]   
% start = [6364 6662];        %Aug[Spring]
start = [9842 9182];        %Aug[Neap]:: 2nd case[7]

V2={};u_avg1={};v_avg1={};
for i = 1:size(Count_6D,2)
%----- N-points Avg. ------%
N1(i) = start(i)+20*1.5;    %3
N2(i) = start(i)+20*4.5;  %4.5

u_avg1{i} = nanmean(u6(:,:,N1(i):N2(i)),3);
v_avg1{i} = nanmean(v6(:,:,N1(i):N2(i)),3);
V2{i} = sqrt(u_avg1{i}.^2 + v_avg1{i}.^2);
%======= change UV =======%
vecx1{i}=u_avg1{i};
vecy1{i}=v_avg1{i};

%----(1) Jul Manual remove max vec.  ----%
% vecx1{i}(abs(vecx1{i})>.0085)=nan;  %.0085
% vecy1{i}(abs(vecy1{i})>.0085)=nan;

%----(2) Aug(x) ----%
end
clear u1 v1

%---- choose the section of the harbor ----%
range_xx = 121.9150; range_xx2=121.9153;
range_yy = 25.1242;  range_yy2=25.1245; 
%% [For Two fig.]
%------ [fix the colorbar] ------%
% bin_edges = 0:5:25;  %old case
bin_edges = 0:5:20;    %Check UV
num_bins = length(bin_edges) - 1;
colors = colormap(nclcmap('precip3_16lev', num_bins+1));
colors = colors(2:end, :);

% time_labels = {'Aug :Ebb Tide(LT) ','Aug :Flood Tide(HT) '};
time_labels = {'Aug :LT [Neap] ','Aug :HT [Neap] '};
high_idx={};idx={};

% Plotting 2 subplots
[ha, hfigure]= tight_subplot(1, 2, [.08 .08], [.15 .05], [.08 .08]);
for ii=1:2
    axes(ha(ii));
if lat>25.124
 Z_pcolorcen(lon,lat,V2{ii});shading flat;hold on;     
else    
 pcolor(lon,lat,V2{ii});shading flat;hold on; 
end
colormap(mycmap);caxis([0 1])
patch([range_xx, range_xx, range_xx2, range_xx2],[range_yy, range_yy2, range_yy2, range_yy], ... 
      'k','FaceAlpha',0,'EdgeColor', 'k','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3); 

[~, bin_idx{ii}] = histc(Count_6D{ii}(:, 3), bin_edges);
    for i = 1:num_bins
     idx{ii} = (bin_idx{ii} == i);
     scatter(Count_6D{ii}(idx{ii}, 1), Count_6D{ii}(idx{ii}, 2), 10, colors(i, :), 'filled');hold on
    %------ Value > 1000 ------%
     high_idx{ii} = (Count_6D{ii}(:, 3) >= 20); %maximum value
     plot(Count_6D{ii}(high_idx{ii}, 1), Count_6D{ii}(high_idx{ii}, 2), 'ro', 'MarkerSize', 1.5); hold on;
    end
%============== Surf. Vectors ===============%    
% scale_factor=0.075; %Jul
% scale_factor=.03;   %Aug
scale_factor=.04; %[Spring/Neap]

for i=1:6:200
for j=1:6:260
    if  ( (lat(i,j)<25.1242 && lon(i,j)>121.9151) ||...
          (lat(i,j)>25.124 && lat(i,j)<25.125 && lon(i,j)>121.9154) ) %remove max vec.
          q1=quiver(lon(i,j),lat(i,j),...
                vecx1{ii}(i,j)*scale_factor,vecy1{ii}(i,j)*scale_factor,...
                'AutoScale','off','color',[1 0.451 0],'Linewidth',1.5);hold on
          q1.ShowArrowHead = 'on';
          q1.MaxHeadSize = 10;      
    end
end
end
q2=quiver(121.9142,25.123,0.005*scale_factor,0,'AutoScale','off','color',[1 0.451 0],'Linewidth',1.5);
q2.ShowArrowHead = 'off';
q2.Marker = '<'
text(121.9142,25.1229,['0.5 cm/s'],'Color','k','Fontsize',16); 
% axis equal
axis image
axis([121.9135 121.9170 25.122 25.1252]) %focus on harbor    
%======== Amount of garbage ========%
text(121.9140,25.1225,['Amount of garbage:', num2str(sumValue{ii})],'Color','k','Fontsize',18)
set(gca, 'Fontsize', 16, 'Linewidth', 1, 'XMinorTick', 'on',...
    'YMinorTick', 'on', 'TickDir', 'out');  
      
% Format title to display 
title({time_labels{ii}},'Interpreter','none');
end
c=colorbar('location','southoutside','Position', [0.15 0.05 0.7 0.02],'colormap', colormap(mycmap))
c.Title.String = "m/s";
caxis([0 1])


return;
% ===== save to .jpg ======%
patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Figure/'
% fn=[patho 'Jul_Tide_Neap'];
fn=[patho 'Aug_Tide_Neap'];
% fn=[patho 'Trash_colorbar']
set(gcf, 'PaperPositionMode' , 'auto' );
print(fn, '-djpeg' , '-r300' )

%% [For One fig.]
%====== change maximum value:: max(Count(:,3)) ========%
% maximum:: 7/12[906] & 7/18[1731] & 7/26[723]
% ------[fix the colorbar]------ 
bin_edges = 0:5:25

num_bins = length(bin_edges) - 1;


colors = colormap(nclcmap('precip3_16lev',num_bins+1));  % GMT_seis /precip3_16lev/hotres
colors=colors(2:end,:);
num_bins = length(bin_edges) -1 ;  
IDX=1;

[~, bin_idx{IDX}] = histc(Count_6D{IDX}(:, 3), bin_edges);

set(0,'defaultfigurecolor','w')
%----- Shading ------%
% for P_color=1
ax1 = axes;
pcolor(ax1, lon, lat, V2{IDX}); shading flat; hold on; 
colormap(ax1, mycmap); axis square;
for i = 1:num_bins
idx{IDX} = (bin_idx{IDX} == i);
scatter(Count_6D{IDX}(idx{IDX}, 1), Count_6D{IDX}(idx{IDX}, 2),...
        20, colors(i, :), 'filled'); hold on;
    
end
    for k=1:size(Count_6D{1},1)
     %====== Count the amount of garbage ======% 
     if (Count_6D{IDX}(k,2)<25.124)...
     || (Count_6D{IDX}(k,2)<25.124 && Count_6D{IDX}(k,1)>121.9155 )
       
        sumValue = sumValue + Count_6D{IDX}(k,3); 
     end     
    end
%------ Value > 25 ------%
high_idx = (Count_6D{IDX}(:, 3) >=25);
plot(Count_6D{IDX}(high_idx, 1), Count_6D{IDX}(high_idx, 2),...
     'ro', 'MarkerSize', 4);hold on; axis square;

scale_factor=0.045;
for i=1:6:200
for j=1:6:260
%     if lat(i,j)<25.1268
%         if lon(i,j)>121.9142
    if  (lat(i,j)<25.124) || (lat(i,j)>25.124 && lon(i,j)>121.9155 )
          q1=quiver(lon(i,j),lat(i,j),...
                vecx1{1}(i,j)*scale_factor,vecy1{1}(i,j)*scale_factor,...
                'AutoScale','off','color',[1 0.451 0],'Linewidth',1.5);hold on
          q1.ShowArrowHead = 'on';
          q1.MaxHeadSize = 10;      
    end
%         end
%     end
end
end
q2=quiver(121.9142,25.123,0.005*scale_factor,0,'AutoScale','off','color',[1 0.451 0]);
q2.ShowArrowHead = 'off';
q2.Marker = '<'
text(121.9142,25.1229,['0.5 cm/s'],'Color','k');
axis equal
axis([121.9135 121.9170 25.122 25.1252]) %focus on harbor
%========= title ===========%
% title({'Jul [7/12 00:00 ~ 7/16 21:00]'})
% title({'High Tide[7/26 00:00 ~ 7/31 00:00]'})
title({'Jul :Low Tide '},'Fontsize',30)
% title({'Aug :Low Tide -3hrs'})
% title({'Aug :High Tide[Spring]-3hrs'})
% title({'Aug :[min. point]'})
% title({'Aug :High Tide[filter tide]'})

text(121.9140,25.1225,['Amount of garbage:', num2str(sumValue)],'Color','k','Fontsize',18)
set(gca, 'Fontsize', 16, 'Linewidth', 1, 'XMinorTick', 'on',...
    'YMinorTick', 'on', 'TickDir', 'out');
caxis([0 1])

return; %do not adding colorbar
c1 = colorbar('location','southoutside','Position', [.32 .03 .4 .02]); 
c1.Title.String = "m/s";
caxis([0 1]);


%====== second colorbar =======%
ax2 = axes('Position', ax1.Position, 'Color', 'none'); 
hold(ax2, 'on');
set(ax2, 'XTick', [], 'YTick', [], 'Visible', 'off'); 
colormap(ax2, colors); % do not affect pcolor
cbar2 = colorbar(ax2, 'location','eastoutside', 'Position', [.75 .1 .02 .8]); 
cbar2.Ticks = linspace(0, 1, num_bins+1); 
cbar2.TickLabels = arrayfun(@(x) sprintf('%d', bin_edges(x)), 1:num_bins+1, 'UniformOutput', false);
cbar2.Label.String = 'Amount of garbage';


% ===== save to .jpg ======%
% patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Figure/'
patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Figure/Check_UV/'
fn=[patho 'Jul_HT'];
% fn=[patho 'Aug_LT_-3hrs'];
% fn=[patho 'Aug_HT_Spring_-3hrs'];
% fn=[patho 'Jul_area'];
set(gcf, 'PaperPositionMode' , 'auto' );
print(fn, '-djpeg' , '-r300' )




