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

%====== Spring & Neap ======%
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Jul_6days/Tide/Spring/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Jul_6days/Tide/Neap/'
% patht='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/Spring/'
% patht='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/Neap/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/Spring/_3hrs/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/Neap/_3hrs/'
%====== V_index ======%
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Jul_6days/Tide/V_index/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/V_index/'

%==== rm M2 & K1 ====%
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Jul_6days/Tide/rm_tide/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Aug_6days/Tide/rm_tide/'

pathi2='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/'
pathi3='/data3/yiwen/matlab_coawst/particle_tracking/Final_report/'

%===== Reading files[for every day] =====%
% load([pathi,'POS_0712_0805_6D.mat'])
% load([pathi,'POS_0804_0907_6D.mat']) 
load([pathi,'POS_0712_0731_6D.mat']);load([patht,'Tide_0712_0731.mat'])
% load([pathi,'POS_0804_0823_6D.mat']);load([pathi,'Tide_0804_0823.mat'])
% load([pathi,'POS_0712_0731_6D.mat']);
% load([pathi,'POS_0804_0823_6D.mat']);

%---- choose 805~828(every 4days) ----%
% allCount=Count(15:end-5);
% allLONLAT=all_lon_lat(15:end-5);
% allCount=Count(1);
% allLONLAT=all_lon_lat(1);

Jul=9:28;
% Aug=4:23;
allLONLAT_LT=all_lon_lat(Jul);
allLONLAT_HT=all_lon_lat2(Jul);
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),6);
% 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']);
% U07=u6;V07=v6;
% load([pathi2,'202408_0913.mat']);
% U08=u6;V08=v6;
% u6=cat(3,U07,U08);v6=cat(3,V07,V08);
% clear U07 V07 U08 V08
load ([pathi3,'blue.mat']);

%======== Time index ==========%
startTime = datetime(2024, 7 ,04, 0, 0, 0); %8/01(spring) / 8/04(V_index) 
endTime = datetime(2024, 7, 31, 23, 54, 0);
timeStep = hours(3);                        %particles raleased timestep
timeSeries = startTime:timeStep:endTime;
disp(timeSeries(LT_idx(1)))
nansum(Count_6D{1}(:,3))

% Compute V2 matrix
start = [3842 3962];  %[LT HT]
End = start + 120;
increment = 480*6;

V2={};u_avg1={};v_avg1={};
for i = 1:size(Count_6D,2)
%     u_avg1 = nanmean(u6(:,:,start:End),3);
%     v_avg1 = nanmean(v6(:,:,start:End),3);
%     V2{i} = sqrt(u_avg1.^2 + v_avg1.^2);
%     
%     start = start + increment;
%     End = start + 120;
%----- N-points Avg. ------%
N1(i) = start(i)+20*1.5;
N2(i) = start(i)+20*3;

u_avg1{i} = (u6(:,:,N1(i))+u6(:,:,N2(i)))./2;
v_avg1{i} = (v6(:,:,N1(i))+v6(:,:,N2(i)))./2;
V2{i} = sqrt(u_avg1{i}.^2 + v_avg1{i}.^2);

start(i) = start(i) + increment;

%======= change UV =======%
vecx1{i}=u_avg1{i}.*(-1);
vecy1{i}=v_avg1{i}.*(-1);

vecx1{i}(abs(vecx1{i})>0.03)=nan;
vecy1{i}(abs(vecy1{i})>0.03)=nan;
end


clear u1 v1
%% [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' )

%%
% ------[fix the colorbar]------ %
bin_edges = 0:5:25;
num_bins = length(bin_edges) - 1;
colors = colormap(nclcmap('precip3_16lev', num_bins+1));
colors = colors(2:end, :);

time_labels = {'Jul :Low Tide','Jul :High Tide'};
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));
pcolor(lon,lat,V2{ii});shading flat;hold on; 
colormap(mycmap);caxis([0 1])
[~, 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), 8.5, colors(i, :), 'filled');hold on
    %------ Value > 1000 ------%
     high_idx{ii} = (Count_6D{ii}(:, 3) >= 25); %max:450 / 260 /145
     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.045;
for i=1:6:200
for j=1:6:260
    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{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([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])






