clc;clear all;close all
addpath('/data3/yiwen/matlab_toolbox/t_tide/')
file1 = '/data3/yiwen/matlab_coawst/particle_tracking/ROMS_log_files/Bitou_WIND/TS_const/Check_UV/20240704_0731/Bitou_ocean_his.nc'
file2 = '/data3/yiwen/matlab_coawst/particle_tracking/ROMS_log_files/Bitou_WIND/TS_const/old_test/20240801_0913/Bitou_ocean_his.nc'
pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/'
pathi2='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/'
pathi3='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Total_Tide/'


for loading=1
%======== Velocity [7/04 00:00 ~ 9/12 23:57]=========%
load([pathi,'202407.mat']);U07=u6;V07=v6;
load([pathi,'202408_0913.mat']);U08=u6;V08=v6;
clear u1 v1
u6=cat(3,U07,U08);v6=cat(3,V07,V08);
clear U07 V07 U08 V08 

lon=ncread(file1,'lon_rho');
lat=ncread(file1,'lat_rho');
%======== SSH [7/04 00:00 ~ 9/12 23:57]=========%
SSH1=ncread(file1,'zeta'); 
SSH2=ncread(file2,'zeta'); 
SSH=cat(3,SSH1,SSH2);
clear SSH1 SSH2
%======== Vwind [7/04 00:00 ~ 9/12 23:57]=========%
Vwind1=ncread(file1,'Vwind'); Vwind2=ncread(file2,'Vwind'); Vwind=cat(3,Vwind1,Vwind2);
clear Vwind1 Vwind2 
Uwind1=ncread(file1,'Uwind'); Uwind2=ncread(file2,'Uwind'); Uwind=cat(3,Uwind1,Uwind2);
clear Uwind1 Uwind2
end

% choose the point of the harbar
i= 108:115; j= 80;
% i=102:108;j=93;

for k=1:length(i)
VVel_M{k}=squeeze(v6(i(k),j,:)).*(-1);
UVel_M{k}=squeeze(u6(i(k),j,:)).*(-1);
SSH_M{k}=squeeze(SSH(i(k),j,:));
VWind_M{k}=squeeze(Vwind(i(k),j,:)).*(-1);
UWind_M{k}=squeeze(Uwind(i(k),j,:)).*(-1);
end
clear u6 v6 SSH Vwind Uwind
% Velo=sqrt(VVel_M{1}.^2+UVel_M{1}.^2); %total velocity
% Velo2=sqrt(VWind_M{1}.^2+UWind_M{1}.^2); 

%======= Bitou_his.nc total time =========%
startTime = datetime(2024, 7, 04, 0, 0, 0);
endTime = datetime(2024, 9, 12, 23, 57, 0);
timeStep = minutes(3); %Roms output timestep
timeSeries = startTime:timeStep:endTime;

%======= Tidal date =======%
load([pathi3, 'Tide_07_09.mat'])
start = 3842;             %7/12(3842) 
End = 34080 ;
indices = start:20*3:End; %every 3hrs(20*3)
SSH_3hrs=SSH_M{1}(indices);
V_3hrs=VVel_M{1}(indices);U_3hrs=UVel_M{1}(indices);
V2_3hrs=VWind_M{1}(indices);U2_3hrs=UWind_M{1}(indices);

startTime = datetime(2024, 7, 12, 0, 0, 0);
endTime = datetime(2024, 9, 12, 23, 55, 0);
timeStep = hours(3); 
timeSeries2 = startTime:timeStep:endTime;


%------- Tide index ------%
for Tide=1
y1 = -0.6; y2 =  0.6; 
x1 = datetime(2024, 7, 12, 0, 0, 0);
x2 = datetime(2024, 7, 16, 0, 0, 0);
x3 = datetime(2024, 7, 25, 0, 0, 0);
x4 = datetime(2024, 7, 29, 0, 0, 0);
x5 = datetime(2024, 8, 01, 0, 0, 0);

x6 = datetime(2024, 8, 04, 0, 0, 0);
x7 = datetime(2024, 8, 09, 0, 0, 0);
x8 = datetime(2024, 8, 18, 12, 0, 0);
x9 = datetime(2024, 8, 23, 0, 0, 0);
x10 = datetime(2024, 8, 24, 0, 0, 0);
Neap=[0.9290 0.6940 0.1250];
Spring=[0.6350 0.0780 0.1840];
end

for tide_testing=1
% %======== T_tide analysis[1] =========%
% % === (High Tide) ===
% [HighTide, HighTimeIdx] = findpeaks(SSH_M{1}); 
% HighTimes = timeSeries(HighTimeIdx); 
% 
% % === (Low Tide) ===
% [LowTide, LowTimeIdx] = findpeaks(-SSH_M{1}); 
% LowTide = -LowTide;       
% LowTimes = timeSeries(LowTimeIdx); 
% 
% %======== T_tide analysis[2] =========%
% mean_u=nanmean(UVel_M{1});
% mean_v=nanmean(VVel_M{1});
% Current=(UVel_M{1}-mean_u)+sqrt(-1).*(VVel_M{1}-mean_v);
% 
% % interval : hrs
% [Name_M2,Fq_M2,Tidecon_M2,Xoui_M2]=t_tide(VVel_M{1},'interval',0.05,'latitude',25.124,...
%                        'start time',[2024,07,04],'rayleigh', ['M2 '; 'S2 '; 'N2 '; 'K2 ']) 
% [Name_K1,Fq_K1,Tidecon_K1,Xoui_K1]=t_tide(VVel_M{1},'interval',0.05,'latitude',25.124,...
%                        'start time',[2024,07,04],'rayleigh',['K1';'O1';'P1';'Q1'])
%                    
% t = datenum(timeSeries);                    
% SSH_pred = t_predic(t, Name_K1, Fq_K1, Tidecon_K1); 
% 
% % === Find High tide ===
% [LowTide, LT_idx] = findpeaks(SSH_pred);     % low tide (極大值)
% [HighTide, HT_idx] = findpeaks(-SSH_pred);   % high tide(極小值)
% HighTide = -HighTide;                        % 轉回正值
% 
% % === transfer time format(hourly) === %
% HighTimes = timeSeries(HT_idx); 
% LowTimes = timeSeries(LT_idx);
end

%%
set(0,'defaultfigurecolor','w')
start= 3842;
% indices = start:60:start+480*5; %every 3 hrs
indices = start:60:size(timeSeries,2);

subplot(2,1,1)
for SSH=1
o1=plot(timeSeries,SSH_M{1},'Linewidth',2);hold on 
h1=plot(timeSeries2(HT_idx), SSH_3hrs(HT_idx), 'o',...
        'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
h2=plot(timeSeries2(LT_idx), SSH_3hrs(LT_idx), 'o',...
        'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'g'); 
%---------%    
h3=fill([x1 x2 x2 x1], [y1 y1 y2 y2], Neap, 'FaceAlpha', 0.3, 'EdgeColor', 'none');  
h4=fill([x2 x3 x3 x2], [y1 y1 y2 y2], Spring, 'FaceAlpha', 0.3, 'EdgeColor', 'none');
fill([x3 x4 x4 x3], [y1 y1 y2 y2], Neap, 'FaceAlpha', 0.3, 'EdgeColor', 'none'); 
fill([x4 x5 x5 x4], [y1 y1 y2 y2], Spring, 'FaceAlpha', 0.3, 'EdgeColor', 'none');
fill([x6 x7 x7 x6], [y1 y1 y2 y2], Neap, 'FaceAlpha', 0.3, 'EdgeColor', 'none');
fill([x7 x8 x8 x7], [y1 y1 y2 y2], Spring, 'FaceAlpha', 0.3, 'EdgeColor', 'none');
fill([x8 x9 x9 x8], [y1 y1 y2 y2], Neap, 'FaceAlpha', 0.3, 'EdgeColor', 'none');
fill([x9 x10 x10 x9], [y1 y1 y2 y2], Spring, 'FaceAlpha', 0.3, 'EdgeColor', 'none');

legend([h1(1) h2(1)],'High tide','Low tide') 
ylabel('m');title('Sea surface height(SSH)')
set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out');ytickformat('%.2f')
set(gca,'fontsize',14,'linewidth',2);grid on;
end
% xlim([datetime(2024, 7, 12, 0, 0, 0), datetime(2024, 8, 24, 0, 0, 0)])
xlim([datetime(2024, 7, 12, 0, 0, 0), datetime(2024, 8, 24, 0, 0, 0)]);xLimits = xlim();    
xticks(xLimits(1):days(5):xLimits(2));xtickformat('MMM dd')
ylim([-.6 .6])

%----------------------------%
subplot(2,1,2)
for r=1
% plot(timeSeries,Var2,'Linewidth',2);hold on
% h1=plot(timeSeries(indices), Var2(indices), 'o',...
%     'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
% ylabel('m/s');title('Surface V-Velocity')
% ylim([-.04 .04])
% set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out');ytickformat('%.3f')
% set(gca,'fontsize',14,'linewidth',2);grid on;
% xlim([datetime(2024, 8, 09, 0, 0, 0), datetime(2024, 8, 20, 0, 0, 0)])
end

C={[0 0.4470 0.7410] [0.8500 0.3250 0.0980] [0.3010 0.7450 0.9330] [0.9290 0.6940 0.1250]}

for Surface_V=1
scale_factor=3;
% X=datenum(timeSeries);
% X2=datenum(timeSeries2);
yyaxis left
P1=plot(timeSeries(indices),VVel_M{1}(indices),'Linewidth',2,'Linestyle','-','color',C{1});hold on
P2=plot(timeSeries(indices),UVel_M{1}(indices),'Linewidth',2,'Linestyle','-','color',C{2});hold on
h1=plot(timeSeries2(HT_idx), V_3hrs(HT_idx), 'o',...
        'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
h2=plot(timeSeries2(LT_idx), V_3hrs(LT_idx), 'o',...
        'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'g');    
ylabel('m/s');
ylim([-.1 .1])

%V2_3hrs=VWind_M{1}
yyaxis right
P3=plot(timeSeries(indices),VWind_M{1}(indices),'Linewidth',2,'Linestyle','-','color',C{3});hold on
P4=plot(timeSeries(indices),UWind_M{1}(indices),'Linewidth',2,'Linestyle','-','color',C{4});hold on
ylabel('m/s');
title('Surface current')
legend([P1(1) P2(1) P3(1) P4(1)],'Surf. V-Velocity','Surf. U-Velocity',...
       'Surf. V-Wind','Surf. U-Wind');legend boxoff
ylim([-60 20])

% datetick('x', 'mmm-dd HH:MM', 'keepticks'); % show ticks
set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out');ytickformat('%.2f')
set(gca,'fontsize',14,'linewidth',2);grid on;
end
% xlim(datenum([datetime(2024, 8, 04, 0, 0, 0), datetime(2024, 8, 24, 0, 0, 0)]));
xlim([datetime(2024, 7, 12, 0, 0, 0), datetime(2024, 8, 24, 0, 0, 0)]);xLimits = xlim();    
xticks(xLimits(1):days(5):xLimits(2));xtickformat('MMM dd')

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

%%
for testing=1
C=colormap(lines);

subplot(2,1,1)
for k=1:size(Vel_M,2)
H{k}=plot(timeSeries,Vel_M{k},'Linewidth',2.5,'Color',C(k,:));hold on
end
            
set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out');ytickformat('%.3f')
set(gca,'fontsize',14,'linewidth',2);grid on;
ylabel('m/s');title('Surface V-Velocity')
xlim([datetime(2024, 7, 04, 0, 0, 0), datetime(2024, 7, 31, 0, 0, 0)])             
ylim([-0.08 0.05])  

subplot(2,1,2)
for k=1:size(Vel_M,2)
H{k}=plot(timeSeries,SSH_M{k},'Linewidth',2.5,'Color',C(k,:));hold on
end
h=legend([H{1}(1),H{2}(1),H{3}(1),H{4}(1),H{5}(1),...
        H{6}(1),H{7}(1),H{8}(1)],'(i,j)=(108,80)','(i,j)=(109,80)','(i,j)=(110,80)',...
                 '(i,j)=(111,80)','(i,j)=(112,80)','(i,j)=(113,80)',...
                 '(i,j)=(114,80)','(i,j)=(115,80)');legend boxoff;
set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out');ytickformat('%.3f')
set(gca,'fontsize',14,'linewidth',2);grid on;
ylabel('m/s');title('Subsurface V-Velocity')
% ylabel('N/m^2');title('Surface V-momentum stress')
xlim([datetime(2024, 7, 04, 0, 0, 0), datetime(2024, 7, 31, 0, 0, 0)])      
ylim([-0.08 0.05])  
end


for region=1
%---- choose the section of the harbor ----%
% range_x = lon(108,80); range_x2= lon(115,80);
% range_y = lat(108,80); range_y2= lat(115,80);   %fix the lat
% 
% 
% %------- figure -------%
% pcolor(lon,lat,Vwind(:,:,1));shading flat;hold on 
% patch([range_x, range_x, range_x2, range_x2],[range_y, range_y2, range_y2, range_y], ... 
%       'b','FaceAlpha',0,'EdgeColor', 'k','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);  
% axis square
end

  
  
  
  
  
  
  