clc;clear all;close all;
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/colorbrewer')
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/m_map')
addpath('/home/yiwen/matlab_test/tight_subplot')

pathi1='/data3/yiwen/matlab_coawst/particle_tracking/1002_tracking/boundary_sp2/'
pathi2='/data3/yiwen/matlab_coawst/particle_tracking/0926_tracking/'
load([pathi2,'spun_up_10days'])

u_avg1=nanmean(u6(:,:,1152),3);v_avg1=nanmean(v6(:,:,1152),3); 
u_avg2=nanmean(u6(:,:,1153:1212),3);v_avg2=nanmean(v6(:,:,1153:1212),3); 
u_avg3=nanmean(u6(:,:,1213:1272),3);v_avg3=nanmean(v6(:,:,1213:1272),3); 
u_avg4=nanmean(u6(:,:,1273:1332),3);v_avg4=nanmean(v6(:,:,1273:1332),3);
V1=sqrt(u_avg1(:,2:end).^2+v_avg1(2:end,:).^2);
V2=sqrt(u_avg2(:,2:end).^2+v_avg2(2:end,:).^2);
V3=sqrt(u_avg3(:,2:end).^2+v_avg3(2:end,:).^2);
V4=sqrt(u_avg4(:,2:end).^2+v_avg4(2:end,:).^2);
Msize=linspace(1,2,17300);
%
x1=linspace(121.9125,121.9130,30);y1=linspace(25.1242,25.1285,30);
[x11,y11]=meshgrid(x1,y1);
% 
x2=linspace(121.9130,121.9175,30);y2=linspace(25.1280,25.1285,30);
[x22,y22]=meshgrid(x2,y2);
x11=reshape(x11,900,1);y11=reshape(y11,900,1);
x22=reshape(x22,900,1);y22=reshape(y22,900,1);


%%
set(0,'defaultfigurecolor','w')
[ha, hfigure]= tight_subplot(1, 4, [0.005 0.005], [0.13 0.05], [0.05 0.01]);
load([pathi1,'case2_0hr_A1']);load([pathi1,'case2_0hr_A2']);
xx=[x11(:,1) xx(:,1:end-1)];yy=[y11(:,1) yy(:,1:end-1)];
xx2=[x22(:,1) xx2(:,1:end-1)];yy2=[y22(:,1) yy2(:,1:end-1)];
xx=xx(3:3:end,:);yy=yy(3:3:end,:);
xx2=xx2(3:3:end,:);yy2=yy2(3:3:end,:);
%  itt=0;
%  for k=1:5:30
%      A(itt*30+1:itt*30+30)=k*30+1:k*30+30;
%      itt=itt+1;
%  end
%  xx=xx(A,:);yy=yy(A,:);
%  xx2=xx2(A,:);yy2=yy2(A,:);
%  
nt=size(xx,2);nt1=linspace(1,nt,nt*10);
for i=1:size(xx,1)
xx_new(i,:) = interp1([1:1:nt],xx(i,:),nt1,'linear');
yy_new(i,:) = interp1([1:1:nt],yy(i,:),nt1,'linear');
xx2_new(i,:) = interp1([1:1:nt],xx2(i,:),nt1,'linear');
yy2_new(i,:) = interp1([1:1:nt],yy2(i,:),nt1,'linear');
end

axes(ha(1))
pcolor(lon(2:end,2:end),lat(2:end,2:end),V1);shading flat;hold on;
m_map=m_colmap('jet',size(xx_new,2));
for k=2:2:size(xx_new,2);
plot(xx_new(:,k), yy_new(:,k) ,'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
plot(xx2_new(:,k), yy2_new(:,k),'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
end 
load slice_cmaps;colormap(cmap)
set(gca,'Fontsize',14);set(gca,'Linewidth',2)
title({'Velocity field' ;'at 4/9 00:00:00' ; 'throw particles'},'Fontsize',16)
axis square;caxis([0 1.])

clear xx_new yy_new xx2_new yy2_new 
%======%
load([pathi1,'case2_6hr_A1']);load([pathi1,'case2_6hr_A2']);
xx=[x11(:,1) xx(:,1:end-1)];yy=[y11(:,1) yy(:,1:end-1)];
xx2=[x22(:,1) xx2(:,1:end-1)];yy2=[y22(:,1) yy2(:,1:end-1)];
xx=xx(3:3:end,:);yy=yy(3:3:end,:);
xx2=xx2(3:3:end,:);yy2=yy2(3:3:end,:);
nt=size(xx,2);nt1=linspace(1,nt,nt*10);

for i=1:size(xx,1)
xx_new(i,:) = interp1([1:1:nt],xx(i,:),nt1,'linear');
yy_new(i,:) = interp1([1:1:nt],yy(i,:),nt1,'linear');
xx2_new(i,:) = interp1([1:1:nt],xx2(i,:),nt1,'linear');
yy2_new(i,:) = interp1([1:1:nt],yy2(i,:),nt1,'linear');
end

axes(ha(2))
pcolor(lon(2:end,2:end),lat(2:end,2:end),V2);shading flat;hold on;
for k=2:2:size(xx_new,2);
plot(xx_new(:,k), yy_new(:,k) ,'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
plot(xx2_new(:,k), yy2_new(:,k),'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
end   
set(gca,'yticklabel',[])
load slice_cmaps;colormap(cmap)
set(gca,'Fontsize',14);set(gca,'Linewidth',2)
title({'Velocity field' ;'at 4/9 06:00:00' ; 'throw particles'},'Fontsize',16)
axis square;caxis([0 1.])


clear xx_new yy_new xx2_new yy2_new 
%=======%
load([pathi1,'case2_12hr_A1']);load([pathi1,'case2_12hr_A2']);
xx=[x11(:,1) xx(:,1:end-1)];yy=[y11(:,1) yy(:,1:end-1)];
xx2=[x22(:,1) xx2(:,1:end-1)];yy2=[y22(:,1) yy2(:,1:end-1)];
xx=xx(3:3:end,:);yy=yy(3:3:end,:);
xx2=xx2(3:3:end,:);yy2=yy2(3:3:end,:);
nt=size(xx,2);nt1=linspace(1,nt,nt*10);
for i=1:size(xx,1)
xx_new(i,:) = interp1([1:1:nt],xx(i,:),nt1,'linear');
yy_new(i,:) = interp1([1:1:nt],yy(i,:),nt1,'linear');
xx2_new(i,:) = interp1([1:1:nt],xx2(i,:),nt1,'linear');
yy2_new(i,:) = interp1([1:1:nt],yy2(i,:),nt1,'linear');
end
axes(ha(3))
pcolor(lon(2:end,2:end),lat(2:end,2:end),V2);shading flat;hold on;
for k=2:2:size(xx_new,2);
plot(xx_new(:,k), yy_new(:,k) ,'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
plot(xx2_new(:,k), yy2_new(:,k),'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
end   
load slice_cmaps;colormap(cmap);set(gca,'yticklabel',[])
set(gca,'Fontsize',14);set(gca,'Linewidth',2)
title({'Velocity field' ;'at 4/9 12:00:00' ; 'throw particles'},'Fontsize',16)
axis square;caxis([0 1.])


clear xx_new yy_new xx2_new yy2_new 
%============%
load([pathi1,'case2_18hr_A1']);load([pathi1,'case2_18hr_A2']);
xx=[x11(:,1) xx(:,1:end-1)];yy=[y11(:,1) yy(:,1:end-1)];
xx2=[x22(:,1) xx2(:,1:end-1)];yy2=[y22(:,1) yy2(:,1:end-1)];
xx=xx(3:3:end,:);yy=yy(3:3:end,:);
xx2=xx2(3:3:end,:);yy2=yy2(3:3:end,:);
nt=size(xx,2);nt1=linspace(1,nt,nt*10);
for i=1:size(xx,1)
xx_new(i,:) = interp1([1:1:nt],xx(i,:),nt1,'linear');
yy_new(i,:) = interp1([1:1:nt],yy(i,:),nt1,'linear');
xx2_new(i,:) = interp1([1:1:nt],xx2(i,:),nt1,'linear');
yy2_new(i,:) = interp1([1:1:nt],yy2(i,:),nt1,'linear');
end
axes(ha(4))
pcolor(lon(2:end,2:end),lat(2:end,2:end),V2);shading flat;hold on;
for k=2:2:size(xx_new,2);
plot(xx_new(:,k), yy_new(:,k) ,'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
plot(xx2_new(:,k), yy2_new(:,k),'o','markerfacecolor', m_map(k,:),...
   'markeredgecolor','none','Markersize',Msize(k));hold on
end   
load slice_cmaps;colormap(cmap);set(gca,'yticklabel',[])
set(gca,'Fontsize',14);set(gca,'Linewidth',2)
title({'Velocity field' ;'at 4/9 18:00:00' ; 'throw particles'},'Fontsize',16)
axis square;caxis([0 1.])

c=colorbar('location','southoutside','Position', [0.23 0.1 0.6 0.02],...
 'colormap', colormap(cmap));
c.Title.String = "m/s"; 

%% ======= save to .jpg ========%
set(gcf, 'PaperPositionMode' , 'auto' ) 
print( 'case2_1020' , '-djpeg' , '-r500' )
% print( 'case2_18hr' , '-djpeg' , '-r500' )







