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/0519_tracking/'
% load([pathi,'202407.mat']) ;U07=u6;V07=v6;
% load([pathi,'202408_0906.mat']) ;U08=u6;V08=v6;
% clear u6 v6 u1 v1
% u6=cat(3,U07,U08);
% v6=cat(3,V07,V08);
load([pathi,'202407.mat'])
% load([pathi,'202408_0913.mat'])
 
mask=u6(:,:,1);
mask(~isnan(mask))=1; %ocn.
mask(isnan(mask))=0;  %land
u6(isnan(u6))=0;
v6(isnan(v6))=0;

lon=lon(:,1);lat=lat(1,:);
%======================== case 1 :: harbor ===============================%
% x1=linspace(121.9147,121.9153,20);y1=linspace(25.1242,25.1250,20); % Original(20*20)
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);
    %Default:60*[10] (time step [mins])
    dtuv=3*60; %[12]                    
%   dtuv = 0.2*60;

%---------------------------------%
    start = 6242;
    dstart_range = start:20*3:start+480*5; %[every three hour/ 1hr=20]
    nstep = 180/dtuv;
    mon1 = 1;
    
    xx = cell(1, length(dstart_range));
    yy = cell(1, length(dstart_range));
%---------------------------------------------------%    
for DS = 1:length(dstart_range)   
    dstart = dstart_range(DS);   
    m = dstart+480*5;          % Running for 6days  
    xpos1(:) = xpos(:, 1, 1);
    ypos1(:) = ypos(:, 1, 1);
%     ditt = 1; 

for d = dstart:1:m    % time step [4/3]
 for tt=1:180/dtuv    %Model & particle DT
     disp(['timestep ', num2str(d), ', substep ', num2str(tt)])
    for i = 1:mf      % loc. index  
      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);

        if isempty(A) || isempty(B) || length(A) == length(lon) || length(B) == length(lat)
            xpos_T(i, d, mon1) = xpos1(i);
            ypos_T(i, d, mon1) = ypos1(i);
        else
        [~, ix] = min(abs(xpos1(i) - lon(A)));
        [~, 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 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
%         xx{DS}(i, ditt) = xpos1(i);
%         yy{DS}(i, ditt) = ypos1(i);
        xx{DS}(i, d-dstart+1) = xpos1(i);
        yy{DS}(i, d-dstart+1) = ypos1(i);
    end 
    end
%     ditt = ditt + 1; 
end
    clear um2 vm2   
end
end

% return;
%====== save to mat =======%l
% patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0113_tracking/3hrs/Jul_5days/'
patho='/data3/yiwen/matlab_coawst/particle_tracking/Offline_output/2025/0519_tracking/3hrs/Jul_5days/'
fn = [patho 'Jul_0717_5D.mat'];
% fn = [patho 'Aug_0804_6D.mat'];
save(fn,'xx','yy')

%%
%  for testmon=testmon0
% %     ditt=1;   
%     xpos1(:)=xpos(:,1,1);ypos1(:)=ylpos(:,1,1);
%     TT=0;
% %------ 0hrs = spunup for 9days(7/12) ------%       
%     %(7/12 00:00=5282 / (7/31 00:00 = 13440) / ) %     
%     dstart=5282; 
%     dstart=dstart:20:m;
%     m=5400;mon1=1;
%     
% for DS=1:length(dstart);
% %    for i=1:mf    
% %      if i>=10 
% %          dstart=dstart+20;
%   if dstart>m ,dstrat=m;end       
% %          ditt=i-10+1+1;
% % %      if i>=10 ,dstart=dstart+20;ditt=1;
% % %        else ,d=1
% %      else, dstart=dstart; 
% %      end
% 
%  for  d=dstart: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);
%   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)));
% 
% % c--------update the drifter positions-------------- %
%  [xpos1(i),ypos1(i)] = partical_B_new(xpos1(i),ypos1(i),um2,vm2,lon,lat',ix,iy,dtuv);
% 
%   end
%     xx(i,ditt)=xpos1(i);yy(i,ditt)=ypos1(i);
% 
% % ditt=ditt+1;  
% 
%   end
%   ditt=ditt+1;  
% %    i
% end
%    
%  end
% 
% end %dstart
