close all;clear all;clc;
datestr(now,'HH:MM')
addpath(genpath('/home/yckuo/MATLAB/m_map/'))
addpath('/data3/yiwen/matlab_coawst/particle_tracking/')
%---------- set the initial drifter positions --------------%
%==================================%
% loading spun_up :: [U & V] 
pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/'
load([pathi,'202407.mat'])  

lon=lon(:,1);lat=lat(1,:);
%===== KAEMI Typhoon ======%
x1=linspace(121.9154,121.9157,10);  y1=linspace(25.1247,25.1249,10);
x2=linspace(121.9154,121.91585,10); y2=linspace(25.1235,25.1239,10);
x3=linspace(121.9157,121.9160,10);  y3=linspace(25.1243,25.1245,10);
x4=linspace(121.91505,121.91535,10);y4=linspace(25.12375,25.12395,10);
x5=linspace(121.91605,121.91615,10);y5=linspace(25.1222,25.1225,10);
x6=linspace(121.9162,121.9163,10);  y6=linspace(25.12235,25.12265,10);

[x11,y11]=meshgrid(x1,y1);
[x22,y22]=meshgrid(x2,y2);
[x33,y33]=meshgrid(x3,y3);
[x44,y44]=meshgrid(x4,y4);
[x55,y55]=meshgrid(x5,y5);
[x66,y66]=meshgrid(x6,y6);
%=========================================================================%
p=1; 
for i=1:size(x11,1)
    for j=1:size(x11,2)    
    xpos(p,1,1)=x11(i,j);ypos(p,1,1)=y11(i,j);
    xpos2(p,1,1)=x22(i,j);ypos2(p,1,1)=y22(i,j);  
    xpos3(p,1,1)=x33(i,j);ypos3(p,1,1)=y33(i,j); 
    xpos4(p,1,1)=x44(i,j);ypos4(p,1,1)=y44(i,j);
    xpos5(p,1,1)=x55(i,j);ypos5(p,1,1)=y55(i,j);
    xpos6(p,1,1)=x66(i,j);ypos6(p,1,1)=y66(i,j);
    p=p+1;
    end
end
    mf=size(xpos,1);
    %Default:60*10 (time step [mins])
    %test for smaller dt  [W&N Bdy] 
%     dtuv=1*60; 
    dtuv=3*60;
    testmon0=1;
for testmon=testmon0
    ditt=1;   
    xpos1(:)=xpos(:,1,1);ypos1(:)=ypos(:,1,1);
    xpos2(:)=xpos2(:,1,1);ypos2(:)=ypos2(:,1,1);
    xpos3(:)=xpos3(:,1,1);ypos3(:)=ypos3(:,1,1);
    xpos4(:)=xpos4(:,1,1);ypos4(:)=ypos4(:,1,1);
    xpos5(:)=xpos5(:,1,1);ypos5(:)=ypos5(:,1,1);
    xpos6(:)=xpos6(:,1,1);ypos6(:)=ypos6(:,1,1);
    TT=0;
%------ 0hrs = spunup for 9days(7/12) ------%       
    dstart = 9602; 
    m = dstart+480*2; 
    mon1=1;     

 for  d=dstart:1:m %timestep   
  for i=1:mf % position          
      um2=u6(:,:,d);vm2=v6(:,:,d);
%-------------------------------------------------------------%
      um2(isnan(um2))=0;vm2(isnan(vm2))=0;
      A=find((xpos1(i)-lon)>=0); B=find((ypos1(i)-lat)>=0);
      A2=find((xpos2(i)-lon)>=0); B2=find((ypos2(i)-lat)>=0);
      A3=find((xpos3(i)-lon)>=0); B3=find((ypos3(i)-lat)>=0);
      A4=find((xpos4(i)-lon)>=0); B4=find((ypos4(i)-lat)>=0);
      A5=find((xpos5(i)-lon)>=0); B5=find((ypos5(i)-lat)>=0);
      A6=find((xpos6(i)-lon)>=0); B6=find((ypos6(i)-lat)>=0);
for N_point=1      
  if (isempty(A==1)|isempty(B==1)|length(A)==length(lon)|length(B)==length(lat) ) 
         xpos(i,d,mon1)=xpos1(i);ypos(i,d,mon1)=ypos1(i);
  else
  [C,ix]=min(abs(xpos1(i)-lon(A)));[C,iy]=min(abs(ypos1(i)-lat(B)));       
  [xpos1(i),ypos1(i)] = partical_B_new(xpos1(i),ypos1(i),um2,vm2,lon,lat',ix,iy,dtuv);       
  end       
         
  if (isempty(A2==1)|isempty(B2==1)|length(A2)==length(lon)|length(B2)==length(lat)) 
         xpos2(i,d,mon1)=xpos2(i);ypos2(i,d,mon1)=ypos2(i);
  else       
  [C2,ix2]=min(abs(xpos2(i)-lon(A2)));[C2,iy2]=min(abs(ypos2(i)-lat(B2)));
  [xpos2(i),ypos2(i)] = partical_B_new(xpos2(i),ypos2(i),um2,vm2,lon,lat',ix2,iy2,dtuv);
  end
  
  if (isempty(A3==1)|isempty(B3==1)|length(A3)==length(lon)|length(B3)==length(lat))
         xpos3(i,d,mon1)=xpos3(i);ypos3(i,d,mon1)=ypos3(i);
  else       
  [C3,ix3]=min(abs(xpos3(i)-lon(A3)));[C3,iy3]=min(abs(ypos3(i)-lat(B3)));
  [xpos3(i),ypos3(i)] = partical_B_new(xpos3(i),ypos3(i),um2,vm2,lon,lat',ix3,iy3,dtuv);    
  end
  
  if (isempty(A4==1)|isempty(B4==1)|length(A4)==length(lon)|length(B4)==length(lat)) 
         xpos4(i,d,mon1)=xpos4(i);ypos4(i,d,mon1)=ypos4(i);
  else
  [C4,ix4]=min(abs(xpos4(i)-lon(A4)));[C4,iy4]=min(abs(ypos4(i)-lat(B4)));
  [xpos4(i),ypos2(i)] = partical_B_new(xpos4(i),ypos4(i),um2,vm2,lon,lat',ix4,iy4,dtuv);   
  end
  
  if (isempty(A5==1)|isempty(B5==1)|length(A5)==length(lon)|length(B5)==length(lat)) 
         xpos5(i,d,mon1)=xpos5(i);ypos5(i,d,mon1)=ypos5(i);
  else
  [C5,ix5]=min(abs(xpos5(i)-lon(A5)));[C5,iy5]=min(abs(ypos5(i)-lat(B5)));
  [xpos5(i),ypos5(i)] = partical_B_new(xpos5(i),ypos5(i),um2,vm2,lon,lat',ix5,iy5,dtuv);  
  end
  
  if (isempty(A6==1)|isempty(B6==1)|length(A6)==length(lon)|length(B6)==length(lat)) 
         xpos6(i,d,mon1)=xpos6(i);ypos6(i,d,mon1)=ypos6(i);         
  else
  [C6,ix6]=min(abs(xpos6(i)-lon(A6)));[C6,iy6]=min(abs(ypos6(i)-lat(B6)));
% c--------update the drifter positions-------------- %
 [xpos6(i),ypos6(i)] = partical_B_new(xpos6(i),ypos6(i),um2,vm2,lon,lat',ix6,iy6,dtuv); 
  end
end 

    xx(i,ditt)=xpos1(i);yy(i,ditt)=ypos1(i);
    xx2(i,ditt)=xpos2(i);yy2(i,ditt)=ypos2(i);
    xx3(i,ditt)=xpos3(i);yy3(i,ditt)=ypos3(i);
    xx4(i,ditt)=xpos4(i);yy4(i,ditt)=ypos4(i);
    xx5(i,ditt)=xpos5(i);yy5(i,ditt)=ypos5(i);
    xx6(i,ditt)=xpos6(i);yy6(i,ditt)=ypos6(i);
  end
  ditt=ditt+1;  
end
   
end

XX={xx,xx2,xx3,xx4,xx5,xx6};
YY={yy,yy2,yy3,yy4,yy5,yy6};
%return;

%====== KAEMI Typhoon ======%
patho = '/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/KAEMI_TYPHOON/'
fn = [patho 'Jul_0724_2D.mat'];
save(fn,'XX','YY','-v7.3')

%%
%---------------------plot figure------------------------------
%surface layer:N=6
u_avg=nanmean(u6,3);v_avg=nanmean(v6,3); 
V6=sqrt(u_avg.^2+v_avg.^2);
%---- choose the section of the harbor ----%
A1_x = 121.9154; A1_x2=121.9157;
A1_y = 25.1247; A1_y2=25.1249;   

A2_x = 121.9154; A2_x2=121.91585;
A2_y = 25.1235; A2_y2=25.1239; 

A3_x = 121.9157; A3_x2=121.9160;
A3_y = 25.1243; A3_y2=25.1245;

A4_x = 121.91505; A4_x2=121.91535;
A4_y = 25.12375; A4_y2=25.12395;

A5_x = 121.91605; A5_x2= 121.91615;
A5_y = 25.1222; A5_y2= 25.1225;

A6_x = 121.9162; A6_x2= 121.9163;
A6_y = 25.12235; A6_y2= 25.12265;
% 
set(0,'defaultfigurecolor','w')
pcolor(lon,lat,V6');shading flat;colorbar;hold on;
patch([A1_x, A1_x, A1_x2, A1_x2],[A1_y, A1_y2, A1_y2, A1_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3); 
patch([A2_x, A2_x, A2_x2, A2_x2],[A2_y, A2_y2, A2_y2, A2_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);   
patch([A3_x, A3_x, A3_x2, A3_x2],[A3_y, A3_y2, A3_y2, A3_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);  
patch([A4_x, A4_x, A4_x2, A4_x2],[A4_y, A4_y2, A4_y2, A4_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);   
patch([A5_x, A5_x, A5_x2, A5_x2],[A5_y, A5_y2, A5_y2, A5_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);   
patch([A6_x, A6_x, A6_x2, A6_x2],[A6_y, A6_y2, A6_y2, A6_y], ... 
      'b','FaceAlpha',0,'EdgeColor', 'r','LineWidth', 2.5,'LineStyle', '-','EdgeAlpha',0.3);    
axis square  



