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/0101_tracking/'
% pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2024/1209_tracking/'
pathi='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/'

% load([pathi,'202301_10D.mat'])  
% load([pathi,'202307_10D.mat'])  
load([pathi,'202407.mat'])  
mask=u6(:,:,1);
mask(~isnan(mask))=1; %ocn.
mask(isnan(mask))=0;  %land

u6(isnan(u6))=0;
v6(isnan(v6))=0;

%----- 0402 -----%
for rep_time=1
% dstart = 9602;     %start date
% m = dstart+480*5;  
% 
% multiplier = 15;
% new_size = size(u6(:,:,dstart:m),3) * multiplier;
% U_expanded = zeros(size(u6,1), size(u6,2), new_size);
% V_expanded = zeros(size(u6,1), size(u6,2), new_size);
% U=u6(:,:,dstart:m);
% V=v6(:,:,dstart:m);
% for i = 1:size(U,3)
%     start_idx = (i-1) * multiplier + 1;
%     end_idx = i * multiplier;
%     U_expanded(:,:,start_idx:end_idx) = repmat(U(:,:,i), [1, 1, multiplier]);
%     V_expanded(:,:,start_idx:end_idx) = repmat(V(:,:,i), [1, 1, multiplier]);
% end
end

lon=lon(:,1);lat=lat(1,:);
%======================== case 1 :: harbor ===============================%
% x1=linspace(121.9147,121.9153,20);y1=linspace(25.1242,25.1250,20);
% x2=linspace(121.9151,121.9154,30);y2=linspace(25.1238,25.1244,30);
x1=linspace(121.9150,121.9153,12);y1=linspace(25.1242,25.1245,12);   % small size(12*12)
[x11,y11]=meshgrid(x1,y1);
%=========================================================================%
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);
    p=p+1;
    end
end
    mf=size(xpos,1);
%     xx(:,1) = xpos(:,1,1);
%     yy(:,1) = ypos(:,1,1);
    
    %Default:60*10 (time step [mins])
    %test for smaller dt  [W&N Bdy] 
    dtuv=3*60;
%     dtuv= 0.2*60; 
    
    ditt=1;   
    xpos1(:)=xpos(:,1,1);
    ypos1(:)=ypos(:,1,1);
%------ 0hrs = spunup for 9-days(7/12) ------%          
    dstart = 8642;       %start date
    m = dstart+480*5;  
    nstep = 180/dtuv;   

%======== tracking start =========%
 for  d=dstart:1:m       %timestep       
     
   for tt=1:180/dtuv     %let Model Dt == particle Dt
       disp(['timestep ', num2str(d), ', substep ', num2str(tt)])
   for i=1:mf             % position     
      um2=v6(:,:,d);
      vm2=u6(:,:,d);
%-------------------------------------------------------------%
      um2(isnan(um2))=0;vm2(isnan(vm2))=0;
      A=find((xpos1(i)-lon)>=0); 
      B=find((ypos1(i)-lat)>=0);
%   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
      [~,ix]=min(abs(xpos1(i)-lon(A)));
      [~,iy]=min(abs(ypos1(i)-lat(B)));
%    end
  % --------update the drifter positions--------- %
     [xpos1(i),ypos1(i)] = partical_B_new(xpos1(i),ypos1(i),um2,vm2,lon,lat',ix,iy,dtuv,mask);
   
    if tt == nstep       
      %========== check value ==========% 
      if d - dstart + 1 > 1
        % [xpos1(i) == xpos(i) && ypos1(i) == ypos(i)]     
        if abs(xpos1(i) - xpos(i)) < 1e-10 && abs(ypos1(i) - ypos(i)) < 1e-10
            disp(['Particle ', num2str(i), ' did not move at timestep ', num2str(d)]);
        end
      end
      %=================================%
%       if d - dstart + 1 > 1+1   
      %---- final pos ----%
      xx_1(i,ditt)=xpos1(i);  %ditt = d - dstart + 1
      yy_1(i,ditt)=ypos1(i);
      xx(i,d - dstart + 1)=xpos1(i);    
      yy(i,d - dstart + 1)=ypos1(i);
%       end
     end  % tt= last rep. time
   end    % i 
    ditt=ditt+1;     
    end  %timestep loop : 12(S)
 end     %timestep loop


%return;
%====== save to mat =======%
% patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/0101_tracking/'
% fn = [patho 'Jul_0709_high.mat'];
% save(fn,'xx','yy','xx2','yy2')

% patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2024/1209_tracking/'
% fn = [patho 'Jul_0712_03.mat'];
% save(fn,'xx','yy')
%====== KAEMI Typhoon ======%
patho = '/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/KAEMI_TYPHOON/'
fn = [patho 'Jul_0721_5D_T.mat'];
save(fn,'xx','yy','xx_1','yy_1','-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);
% % 
% set(0,'defaultfigurecolor','w')
% pcolor(lon,lat,V6');shading flat;colorbar;hold on;
% for index = 1:3:900
% %    plot(xx(index,1),yy(index,1),'o','Color','k','LineWidth',0.01);hold on
% %    plot(xx2(index,1),yy2(index,1),'o','Color','k','LineWidth',0.01);hold on
% %    plot(xx(index,end),yy(index,end),'Color','r','LineWidth',0.8);hold on  
%    plot(xx(index,:),yy(index,:),'Color','r','LineWidth',0.01);hold on 
%    plot(xx2(index,:),yy2(index,:),'Color','y','LineWidth',0.01);hold on 
% end
% axis square
% 
% % for index = 1:1:300
% % plot(xx(index,1),yy(index,1),'o','Color','k','Markersize',1);hold on
% % end
% colormap(brewermap(20, 'BuGn'));
% title({'801-900 pos'})
% set(gca,'Fontsize',14);set(gca,'Linewidth',2)
% caxis([0 1.]);axis square
