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/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/'
pathi2='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/'
pathi3='/data3/yiwen/matlab_coawst/particle_tracking/Final_report/'

%===== Reading files =====%
load([pathi,'POS_0804_0907_5D.mat'])   %for every day
% load([pathi,'POS_0712_0805_5D.mat'])
% load([pathi,'POS_0804_5D.mat']);allCount=Count;

load([pathi2,'202408_0914.mat']);
% load([pathi2,'202407.mat']);
load ([pathi3,'blue.mat']);

%===== refer to the date of POS.mat =====%
start_time = datetime(2024, 8, 04, 0, 0, 0); %change time [1-1] #If change .mat file#
end_time = datetime(2024, 9, 13, 23, 57, 0);
interval = hours(3);
T_series = start_time:interval:end_time;
%----- find the time index ------%
find_index = @(input_date) find(T_series == datetime(input_date, 'InputFormat', 'yyyy/MM/dd HH:mm'));
input_date = '2024/08/19 00:00';             %change time [1-2]
index = find_index(input_date)

% Count_6D=allCount(index:2:index+8*3); %1Days(8)
Count_6D=allCount(index:1:size(allCount,2)); %1Days(8)

% Compute V2 matrix
start = 8642;      %Start time
increment = 20*3;  %1D(480)[20*3]

V2={};u_avg1={};v_avg1={};
start_time = datetime(2024, 8, 19, 0, 0, 0); %change time [2-1]
end_time   = datetime(2024, 8, 21, 12, 0, 0);  %change time [2-2]
interval = hours(3); %to check every fig.
% interval = hours(6); %for report fig.
time_series = start_time:interval:end_time;


for i = 1:length(time_series)
%     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 + 20*3;  % 20*6 hrs
%----- N-points Avg. ------%
N1(i) = start+20*3;
N2(i) = start+20*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);

start = start + increment;
end

vecx1=u_avg1;
vecy1=v_avg1;

clear u1 v1 
%%
set(0,'defaultfigurecolor','w')
bin_edges = 0:5:25
num_bins = length(bin_edges) - 1;
colors = colormap(nclcmap('precip3_16lev', num_bins+1));
colors = colors(2:end, :);
high_idx={};idx={};
sumValue=zeros(1,21)

% Plotting N subplots
% [ha, hfigure]= tight_subplot(2, 5, [.05 .035], [.13 .05], [.04 .04]);
[ha, hfigure]= tight_subplot(3, 7, [.05 .035], [.13 .05], [.04 .04]);
for ii=1:size(u_avg1,2)
    axes(ha(ii));
    pcolor(lon,lat,V2{ii});shading flat;hold on; 
    colormap(mycmap);axis square;
    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), 7, colors(i, :), 'filled');hold on
        %------ Value > N ------%
         high_idx{ii} = (Count_6D{ii}(:, 3) >= 15); %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;
         hold on;axis square          
        end
    for k=1:size(Count_6D{ii},1)
     %====== Count the amount of garbage ======% 
     if (Count_6D{ii}(k,2)<25.124)...
     || (Count_6D{ii}(k,2)<25.125 && Count_6D{ii}(k,1)>121.9155)
       
        sumValue(ii) = sumValue(ii) + Count_6D{ii}(k,3); 
     end     
    end 
%======= Vectors =======% 
% scale_factor=0.009; %for 3hrs fig.
scale_factor=0.003;   %for 6hrs fig.
for i=1:8:200
for j=1:8:260
    if lat(i,j)<25.1268
        if lon(i,j)>121.9142
          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.1);hold on
          q1.ShowArrowHead = 'on';
          q1.MaxHeadSize = 6;       
        end
    end
end
end
q2=quiver(121.9145,25.123,0.075*scale_factor,0,'AutoScale','off','color',[1 0.451 0]);
q2.ShowArrowHead = 'off';
q2.Marker = '<'
text(121.9142,25.1229,['7.5 cm/s'],'Color','k');
axis equal
axis([121.9135 121.9170 25.122 25.1252]) %focus on harbor
%============%
    title(sprintf('[8/%d %02d:00]',...
          day(time_series(ii)), hour(time_series(ii))),'fontsize',5);
    text(121.9140,25.1225,{['Amount of'],['garbage:', num2str(sumValue(ii))]},'Color','k','Fontsize',10,'linewidth',1.2)  
axis equal
axis([121.9135 121.9170 25.122 25.1252]) %focus on harbor
set(gca,'fontsize',8,'linewidth',1.2);
end
c=colorbar('location','southoutside','Position', [0.15 .07 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/0113_tracking/3hrs/Figure/'
patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Figure/'
fn=[patho 'Typhoon_0723_sub'];
set(gcf, 'PaperPositionMode' , 'auto' );
print(fn, '-djpeg' , '-r300' )
%%
figure(2)
for i = 1:num_bins
 idx{1} = (bin_idx{1} == i);
 scatter(Count_6D{1}(idx{1}, 1), Count_6D{1}(idx{1}, 2), 50, colors(i, :), 'filled');
 hold on;axis square 
end

colormap(colors);
cbar2 = colorbar;
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';
set(gca,'Fontsize',12,'Linewidth',1,'XMinorTick','on','YMinorTick','on','TickDir','out');


















