%%PLOTFIGURES reproduces the figures in:
%
%Timothy T.K. Chan, Luis Blay Esteban, Sander G. Huisman, John Shrimpton and Bharathram Ganapathisubramani
%Settling behaviour of thin curved particles in quiescent fluid and turbulence
%Journal of Fluid Mechanics
%
%Date: 11 June, 2021

clearvars;

load('Data.mat')

%% Figure format
plot_color = [6 75 145; 76 173 212; 225 151 45; 190 0 0]./255;
plot_green = [58 130 0]./255;
markersymbol = {'s','d','^','o'};
markersz = [5 4.75 4.75 4.75];  %marker size
MarkerLineWidth = [1 1 1 1];
LineWidth = 0.75;

%%
%Figure 1
figure
phase1 = loglog(Quiescent.P_SR(:,1),Quiescent.P_SR(:,2),'k--','LineWidth',1);
hold on
phase2 = loglog(Quiescent.D_SP(:,1),Quiescent.D_SP(:,2),'k:','LineWidth',1);
phase3 = loglog(Quiescent.P_RA(:,1),Quiescent.P_RA(:,2),'k--','LineWidth',1);
phase4 = loglog(Quiescent.D_PC(:,1),Quiescent.D_PC(:,2),'k:','LineWidth',1);
loglog(Quiescent.D_CT(:,1),Quiescent.D_CT(:,2),'k:','LineWidth',1);
for i = 1:length(Quiescent.Re)
    phase_expt = loglog(Quiescent.Re(i),Quiescent.I_star(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
                   'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
end
legend([phase1 phase2],{'Plate','Disk'});
xlabel('$Re$','Interpreter','latex')
ylabel('$I^*$','Interpreter','latex')

%Figure 3
line_color = [191 68 142; 174 172 81; 218 78 86; 187 174 51; 84 53 134; 159 73 47;...
              214 145 56; 205 102 123; 210 73 121; 109 113 216; 210 73 121; 228 124 87;...
              215 123 183; 184 129 210; 164 124 57; 129 40 92; 94 139 213; 111 185 92;...
              189 88 185; 157 54 64; 101 176 111; 66 200 127; 133 180 63; 67 200 172; 72 112 34]./255;
figure
hold on
for i = 1:size(Quiescent.x_nondim_fig3,2)
    plot3(Quiescent.x_nondim_fig3(:,i),Quiescent.y_nondim_fig3(:,i),Quiescent.z_nondim_fig3(:,i),'-','color',line_color(i,:),'LineWidth',0.75)
end
grid on
zlim([0 inf])
xlabel('$x/D$','interpreter','latex')
ylabel('$y/D$','interpreter','latex')
zlabel('$z/D$','interpreter','latex')
ax=gca;
ax.ZDir='reverse';
ax.YDir='reverse';
view(10,10)
xlim([-2 2])
ylim([-2 2])
zlim([0 9])

%Figure 4
%(a),(b)
figure
plot3(Quiescent.traj_fig4(:,1)./Quiescent.D(2),Quiescent.traj_fig4(:,2)./Quiescent.D(2),(Quiescent.traj_fig4(:,3))./Quiescent.D(2),'k-','Linewidth',LineWidth)
xlabel('x/D')
ylabel('y/D')
zlabel('z/D')
view([-160 40]);
axis equal
grid on
ax=gca;
ax.ZDir='reverse';
ax.YDir='reverse';
zlim_plot = zlim;
%(c)
figure
subplot(3,1,1)
    plot(Quiescent.t_fig4,Quiescent.rh_fig4./Quiescent.D(2),'k-','Linewidth',LineWidth)
    ylabel('r_h/D')
    xlabel('t (s)')
subplot(3,1,2)
    plot(Quiescent.t_fig4,(Quiescent.traj_fig4(1:length(Quiescent.t_fig4),3))./Quiescent.D(2),'k-','Linewidth',LineWidth)
    ylabel('z/D')
    xlabel('t (s)')
subplot(3,1,3)
    plot(Quiescent.t_fig4,Quiescent.alpha_fig4,'k-','Linewidth',LineWidth)  %new. More accurate: 05-03-2021
    ylabel(['\alpha (' char(176) ')'])
    xlabel('t (s)')
%(d),(e)
figure
subplot(1,2,1)
    plot(Quiescent.V_fig4(:,1),Quiescent.V_fig4(:,2),'k-','Linewidth',LineWidth)
    xlabel('V_x/\langle V_z \rangle')
    ylabel('V_y/\langle V_z \rangle')
    axis equal
subplot(1,2,2)
    plot(Quiescent.Vh_fig4,Quiescent.V_fig4(:,3),'k-','Linewidth',LineWidth)
    xlabel('V_h/\langle V_z \rangle')
    ylabel('V_z/\langle V_z \rangle')
    axis equal
    
%Figure 5
figure
errorbar(Quiescent.D,Quiescent.Vel_Z,Quiescent.Std_VelZ,'k','Linewidth',LineWidth,'LineStyle','none')
hold on
    for i = 1:length(Quiescent.Vel_Z)
        plot(Quiescent.D(i),Quiescent.Vel_Z(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
set(gca,'fontsize',12)
xlabel('$D$ (mm)','FontSize',12,'interpreter','latex')
ylabel('$\langle V_z \rangle$ (mm/s)','FontSize',12,'interpreter','latex')

%Figure 6
figure
errorbar(Quiescent.Ar,Quiescent.Re,Quiescent.Std_Re,Quiescent.Std_Re,zeros(size(Quiescent.Ar)),zeros(size(Quiescent.Ar)),'k','LineWidth',LineWidth,'LineStyle','none')
hold on
for i = 1:length(Quiescent.Vel_Z)
    plot(Quiescent.Ar(i),Quiescent.Re(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
end
xlim([0.5E4 3E4])
xlabel('$Ar$','interpreter','latex')
ylabel('$Re$','interpreter','latex')

%Figure 7
%(a)-(d)
figure
subplot(2,2,1)
    errorbar(Quiescent.Ar,Quiescent.Lpend_nondim,Quiescent.Std_Lpendnondim,'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    for i = 1:length(Quiescent.Lpend_nondim)
        plot(Quiescent.Ar(i),Quiescent.Lpend_nondim(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
            'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    ylim([0.5 1.2])
    xlim([0.5E4 3E4])
    xlabel('Ar')
    ylabel('L_{pend}/D')
subplot(2,2,2)
    errorbar(Quiescent.Ar,Quiescent.alphamax,Quiescent.Std_alphamax,Quiescent.Std_alphamax,zeros(size(Quiescent.Ar)),zeros(size(Quiescent.Ar)),'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    for i = 1:length(Quiescent.alphamax)
        plot(Quiescent.Ar(i),Quiescent.alphamax(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
            'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    xlim([0.5E4 3E4])
    xlabel('Ar')
    ylabel(['\alpha_{max} (' char(176) ')'])
subplot(2,2,3)
    errorbar(Quiescent.Ar,Quiescent.A_nondim,Quiescent.Std_Anondim,'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    for i = 1:length(Quiescent.A_nondim)
        plot(Quiescent.Ar(i),Quiescent.A_nondim(i),markersymbol{i},'Color',plot_color(i,:),...
            'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    set(gca,'fontsize',12)
    xlabel('Ar')
    ylabel('A/D','FontSize',12)
subplot(2,2,4)
    hold on
    errorbar(Quiescent.Ar,Quiescent.St_r,Quiescent.Std_Str,'k','LineWidth',LineWidth,'LineStyle','none')
    for i = 1:length(Quiescent.Ar)
        plot(Quiescent.Ar(i),Quiescent.St_r(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    ylabel('St')
    xlabel('Ar')
%Inset of (b)
figure
errorbar(Quiescent.Ar,Quiescent.dAngle_mean,Quiescent.Std_dAngle,'k','Linewidth',LineWidth,'LineStyle','none')
hold on
    for i = 1:length(Quiescent.dAngle_mean)
        plot(Quiescent.Ar(i),Quiescent.dAngle_mean(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
          'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
xlim([0.5E4 3E4])
xlabel('Ar')
ylabel(['\Delta\alpha (' char(176) ')'])

%Figure 8
%(a),(b)
figure
subplot(1,2,1)
    errorbar(Quiescent.Ar,Quiescent.Vhmax,Quiescent.Std_Vhmax,'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    for i = 1:length(Quiescent.Ar)
        plot(Quiescent.Ar(i),Quiescent.Vhmax(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    xlim([0.5E4 3E4])
    ylim([1.5 2.7])
    yticks(1.5:0.3:2.7)
    set(gca,'fontsize',12)
    xlabel('$Ar$','FontSize',12,'interpreter','latex')
    ylabel('$V_{h,max}/\langle V_z \rangle$','interpreter','latex')
subplot(1,2,2)
    errorbar(Quiescent.Ar,Quiescent.amax,Quiescent.Std_amax,Quiescent.Std_amax,zeros(size(Quiescent.Ar)),zeros(size(Quiescent.Ar)),'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    errorbar(Quiescent.Ar,Quiescent.amax_pend,Quiescent.Stdamax_pend,Quiescent.Stdamax_pend,zeros(size(Quiescent.Ar)),zeros(size(Quiescent.Ar)),'k','LineWidth',LineWidth,'LineStyle','none')
    for i = 1:length(Quiescent.amax)
        plot(Quiescent.Ar(i),Quiescent.amax(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
        plot(Quiescent.Ar(i),Quiescent.amax_pend(i),markersymbol{i},'Color',[100 100 100]./255,'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',[100 100 100]./255);
    end
    xlim([0.5E4 3E4])
    xlabel('$Ar$','interpreter','latex')
    ylabel('$a_{z,max} (mm/s^2)$','interpreter','latex')
%(b) inset
figure
    errorbar(Quiescent.amax,Quiescent.alphamax,Quiescent.Std_alphamax,Quiescent.Std_alphamax,Quiescent.Std_amax,Quiescent.Std_amax,'k','LineWidth',LineWidth,'LineStyle','none')
    hold on
    for i = 1:length(Quiescent.alphamax)
        plot(Quiescent.amax(i),Quiescent.alphamax(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
    end
    plot([400 1000],0.0274.*[400 1000]+20.77223,'k-')
    ylim([30 50])
    xlabel('a_{z,max} (mm/s^2)')
    ylabel(['\alpha_{max} (' char(176) ')'])
    
%Figure 9
figure
hold on
plot(Quiescent.rh_fig9,Quiescent.z_nondim_fig9,'.','Color',plot_color(2,:),'MarkerSize',0.1)
plot(Quiescent.rh_mod_fig9,Quiescent.z_nondim_mod_fig9,'-','Color',[100 100 100]./255,'LineWidth',LineWidth)
xlabel('r_h/D')
ylabel('z/D')
xlim([-0.8 0.8])
ylim([0 9])
set(gca,'Ydir','reverse')

%Figure 10
figure
errorbar(Quiescent.Ar,Quiescent.C,Quiescent.Std_C,'k','LineWidth',LineWidth,'LineStyle','none','LineStyle','none')
hold on
for i = 1:length(Quiescent.Ar)
    plot(Quiescent.Ar(i),Quiescent.C(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
        'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
end
xlim([0.5E4 3E4])
ylim([0.12 0.28])
xlabel('Ar')
ylabel('C')
%inset
figure
errorbar(Quiescent.Ar,Quiescent.ma_star,Quiescent.Std_ma_star,'k','LineWidth',LineWidth,'LineStyle','none')
hold on
for i = 1:length(Quiescent.Ar)
    plot(Quiescent.Ar(i),Quiescent.ma_star(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
        'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
end
xlim([0.5E4 3E4])
ylim([2 6])
xlabel('Ar')
ylabel('m_a^*')

%Figure 11 (Replace SD with RR/VD/LG for other types of events)
figure
plot3(Turb.traj_SD(:,1),Turb.traj_SD(:,2),Turb.traj_SD(:,3),'k-','Linewidth',LineWidth)
hold on
plot3(Turb.traj_SD(Turb.event_SD,1),Turb.traj_SD(Turb.event_SD,2),Turb.traj_SD(Turb.event_SD,3),'ks','MarkerSize',2)
xlabel('x/D')
ylabel('y/D')
zlabel('z/D')
view([-120 10]);
axis equal
grid on
ax=gca;
ax.ZDir='reverse';
ax.YDir='reverse';

%Figure 12
figure
%(a)
subplot(1,2,1)
markersampleQ = 2:2:min(sum(Quiescent.edges_rnondim<1,1));   %downsample marker positions
markersampleT = 1:1:length(Turb.edges_rnondim);   %downsample marker positions
subplot(1,2,1)
    hold on
    plot(Quiescent.edges_rnondim(:,1),Quiescent.N_rnondim(:,1),'-','Color',plot_color(3,:),'Linewidth',LineWidth);
        plot(Quiescent.edges_rnondim(markersampleQ,1),Quiescent.N_rnondim(markersampleQ,1),markersymbol{3},...
            'Color',plot_color(3,:),'MarkerSize',3,'MarkerFaceColor',plot_color(3,:))
    for i = [1 3]
        plot(Turb.edges_rnondim(:,i),Turb.N_rnondim(:,i),'--','Color',plot_color(i,:),'Linewidth',LineWidth);
        plot(Turb.edges_rnondim(markersampleT,i),Turb.N_rnondim(markersampleT,i),markersymbol{i},...
            'Color',plot_color(i,:),'MarkerSize',3,'MarkerFaceColor','w')
    end
    xlim([0 6])
    ylim([0 0.9])
    xlabel('r_h/D')
    ylabel('PDF')
%(b)
subplot(1,2,2)
    errorbar(Quiescent.Ar,Quiescent.KE_fluc,Quiescent.Std_KEfluc,'k-','Linewidth',LineWidth,'LineStyle','none')
    hold on
        for i = 1:length(Quiescent.Ar)
            plot(Quiescent.Ar(i),Quiescent.KE_fluc(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
                'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
        end
    errorbar(Turb.Ar,Turb.KE_fluc,Turb.Std_KEfluc,Turb.Std_KEfluc,Turb.Std_Ar,Turb.Std_Ar,'k-','Linewidth',LineWidth,'LineStyle','none')
        for i = 1:length(Turb.Ar)
            plot(Turb.Ar(i),Turb.KE_fluc(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),...
                'Linewidth',MarkerLineWidth(i),'MarkerFaceColor','none');
        end
    set(gca,'fontsize',12)
    ax = gca;
    ax.YAxis.Exponent = 2;
    xlabel('Ar','FontSize',14)
    ylabel('E_{fluc} (mm^2/s^2)','FontSize',14)
    ylim([500 900])

%(a),inset
figure
hold on
for i = 1:size(Turb.edges_rnondim,2)
    plot(Turb.edges_rnondim(:,i),Turb.N_rnondim(:,i),'--','Color',plot_color(i,:),'Linewidth',LineWidth);
    plot(Turb.edges_rnondim(markersampleT,i),Turb.N_rnondim(markersampleT,i),markersymbol{i},...
        'Color',plot_color(i,:),'MarkerSize',3,'MarkerFaceColor','w')
end
xlim([0 4.5])
ylim([0 0.4])
xlabel('r_h/D')
ylabel('PDF')

%Figure 13
figure
%(a)
subplot(1,2,1)
    errorbar(Quiescent.D,Quiescent.Vel_Z,Quiescent.Std_VelZ,'k-','Linewidth',LineWidth,'LineStyle','none')
    hold on
        for i = 1:length(Quiescent.D)
            plot(Quiescent.D(i),Quiescent.Vel_Z(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
        end
    errorbar(Quiescent.D,Turb.Vel_Z,Turb.Std_VelZ,'k-','LineStyle','none')
        for i = 1:length(Quiescent.D)
            plot(Quiescent.D(i),Turb.Vel_Z(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor','none');
        end
    set(gca,'fontsize',12)
    xlabel('D (mm)','FontSize',12)
    ylabel('\langle V_z \rangle (mm/s)','FontSize',12,'interpreter','tex')
%(b)
subplot(1,2,2)
    errorbar(Quiescent.Ar,Quiescent.Re,Quiescent.Std_Re,Quiescent.Std_Re,zeros(size(Quiescent.Ar)),zeros(size(Quiescent.Ar)),'k-','Linewidth',LineWidth,'LineStyle','none')
    hold on
        for i = 1:length(Quiescent.Ar)
            plot(Quiescent.Ar(i),Quiescent.Re(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor',plot_color(i,:));
        end
    errorbar(Turb.Ar,Turb.Re,Turb.Std_Re,Turb.Std_Re,Turb.Std_Ar,Turb.Std_Ar,'k--','LineStyle','none')
        for i = 1:length(Turb.Ar)
            plot(Turb.Ar(i),Turb.Re(i),markersymbol{i},'Color',plot_color(i,:),'MarkerSize',markersz(i),'Linewidth',MarkerLineWidth(i),'MarkerFaceColor','none');
        end
    xlim([0.5E4 3E4])
    ylim([1000 3500])
    xlabel('Ar')
    ylabel('Re')
    
%Figure 14
%(a)
figure
for i = 1:4
    subplot(2,2,i)
    plot(Turb.fevent_norm(:,i),Turb.Vzevent_norm(:,i),'.','color',plot_color(i,:))
    xlabel('$f_{event}/f_q$','interpreter','latex')
    ylabel('$V_{event}/\langle V_q\rangle$','interpreter','latex')
    title(['Particle No. ' num2str(i)])
    xlim([0 4])
end
%(b)
figure
plot(Turb.feventbin,Turb.Vzbin,'ks','markersize',2)
xlim([0 3.5])
xlabel('$f_{event}/f_q$','interpreter','latex')
ylabel('$\langle V_{event}\rangle/\langle V_q\rangle$','interpreter','latex')

%Figure 15 (Replace SD with RR/VD/LG for other types of events)
figure
plot(Turb.Vz_SD,Turb.traj_SD(1:length(Turb.Vz_SD),3),'k-','Linewidth',LineWidth)
hold on
for i = 1:(length(Turb.event_SD)-1)
    plot(ones(1,2).*Turb.Vzevent_SD(i),[min(Turb.traj_SD(Turb.event_SD(i):Turb.event_SD(i+1),3)) max(Turb.traj_SD(Turb.event_SD(i):Turb.event_SD(i+1),3))],...
        'k--','Linewidth',LineWidth+0.25)
end
axis manual
plot([1 1],[0 1000],'Color',[0.5 0.5 0.5],'Linewidth',LineWidth-0.25)
xlabel('$V_z/\langle V_q \rangle$','interpreter','latex')
ylabel('z/D')
ax=gca;
ax.YDir='reverse';

%Figure 16
figure
plot(Turb.Veventbin,Turb.N_Vevent,'ks-','LineWidth',LineWidth,'MarkerSize',3,'MarkerFaceColor','none')
xlabel('$V_{event}/\langle V_q \rangle$','interpreter','latex')    % xlabel('U_{z,event}/U_{z,q}')
ylabel('PDF')

%Figure 17
nsep = 30;
%(a)
figure
imagesc(PIV.urms.x(2:end-1)-(PIV.urms.x(1)+PIV.urms.x(2))./2,PIV.urms.z(2:end-1)-(PIV.urms.z(1)+PIV.urms.z(2))./2,PIV.urms.data(2:end-1,2:end-1))
hold on
colormap(flipud(hot))
h = colorbar;
xlabel('x (mm)')
ylabel('z (mm)')
axis equal
xlim([0  150.9429])
ylim([0  226.4143])
caxis([14 18])
xticks([0 50 100 150])
yticks([0 50 100 150 200])
set(h,'YTick',[14 15 16 17 18])
set(get(h,'label'),'string','(|u_x| + 2|u_z|)/3 (mm/s)');
%(b)
UMEAN = (abs(PIV.umean.vx) + 2.*abs(PIV.umean.vz))'./3.*1000;
umean_xt = PIV.umean.vx';
umean_yt = PIV.umean.vz';
figure
imagesc(PIV.umean.x(2:end-1)-(PIV.umean.x(1)+PIV.umean.x(2))./2,PIV.umean.z(2:end-1)-(PIV.umean.z(1)+PIV.umean.z(2))./2,UMEAN(2:end-1,2:end-1))
hold on
quiver(PIV.umean.x(1+nsep:nsep:end-1)-(PIV.umean.x(1)+PIV.umean.x(2))./2,PIV.umean.z(2:nsep:end-1)-(PIV.umean.z(1)+PIV.umean.z(2))./2,...
    umean_xt(2:nsep:end-1,1+nsep:nsep:end-1),umean_yt(2:nsep:end-1,1+nsep:nsep:end-1))
colormap(flipud(hot))
h = colorbar;
xlabel('x (mm)')
ylabel('z (mm)')
axis equal
xlim([0  150.9429])
ylim([0  226.4143])
caxis([0 inf])
xticks([0 50 100 150])
yticks([0 50 100 150 200])
set(get(h,'label'),'string','(|U_x| + 2|U_z|)/3 (mm/s)');
set(h,'YTick',[0 2 4 6])

%Figure 18 (for the other directions, simply replace xx with xz/zx/zz)
figure
plot(PIV.Autocorr_xx.r,PIV.Autocorr_xx.data)
xlabel('r (mm)')
ylabel('\rho_{ij}')