% clc; clear all;

We = 540; lev = 11;
R0 = 0.05;
V0 = 0.5*sqrt(10*We);
L0 = 2;

t_i = 0; 
% t_i = 103.05; 
t_f = 0.0025;
% t_f = max(drop*2_data(length(drop_data),2), t_i+3); 
t_c = t_f;

del_11 = L0/2^14;

N_bins = 74; N_ens = 3; del = L0/2^lev;
dt = 1e-4; N_t = (t_f - t_i)/dt + 1;
edges_size = (2*del_11) * (0.05/(2*del_11))^(1/N_bins) .^ (0:N_bins);
edges_vel = linspace(0.1*V0, 5*V0, N_bins+1);

% edges_size = linspace(2*del_11, 0.05, N_bins+1);
for i=1:N_bins
    r_binned(i) = sqrt(edges_size(i)*edges_size(i+1));
%     r_binned(i) = 0.5*(edges_size(i)+edges_size(i+1));
%     v_binned(i) = sqrt(edges_vel(i)*edges_vel(i+1));
    v_binned(i) = 0.5*(edges_vel(i)+edges_vel(i+1));
end

for i=1:N_bins
  for j=1:N_ens
    N_c_ind(i,j) = 0; u_ave_ind(i,j) = 0; N_c_ind_vel(i,j) = 0; 
  end
end

for i=1:N_bins
    pdf_c(i) = 0; u_ave(i) = 0; uy_ave(i) = 0; ur_ave(i) = 0; 
    ur_ave_cf(i) = 0; uy_ave_cf(i) = 0; pdf_cf(i) = 0; N_cf(i) = 0;
    ur_var_cf(i) = 0; uy_var_cf(i) = 0; pdf_var_cf(i) = 0;

    pdf_c_vel(i) = 0;
end

numcount_global = 0; d_ave_global = 0; std_err_global = 0;

for ens=1:N_ens
% ens = 1;
filepath=strcat('/mnt/ResearchData_a/Cylinder Collision/Uniform_Disturb/Unsymmetric Disturb/Gravity_test/We_540_eta_0.06_Bo_125_unsym_L11_Nmax_80_enlarged_new/ens_', num2str(ens), '/');
% filepath=strcat('/mnt/ResearchData_a/Cylinder Collision/Uniform_Disturb/Unsymmetric Disturb/Gravity_test/We_280_eta_0.06_Bo_125_unsym_L11_enlarged_new/ens_2/');

drop_data = importdata(filepath + "droplets.dat");
drop_data = [zeros(size(drop_data,1),1) drop_data];
drop_data(:,11) = (drop_data(:,4)./(4*pi/3)).^(1/3);
% drop_data(:,12) = sqrt(drop_data(:,8).^2 + drop_data(:,9).^2 + drop_data(:,10).^2);
drop_data(:,12) = sqrt(drop_data(:,9).^2);



dnum_pt = 0; pt_f = 0; tmp_10 = 0; tmp_20 = 0; tmp_30 = 0;  numcount = 0;

for i=1:length(drop_data)
    tmp_10 = tmp_10 + drop_data(i,11);
    tmp_20 = tmp_20 + (drop_data(i,11))^2;
    tmp_30 = tmp_30 + (drop_data(i,11))^3;

    if (drop_data(i,11) > del) numcount = numcount + 1; end
    if ((drop_data(i,11) > 8*del) && (drop_data(i,11) < 0.2) && (drop_data(i,2) == t_c))
        numcount_global = numcount_global + 1;
        d_ave_global = d_ave_global + drop_data(i,11);
        std_err_global = std_err_global + drop_data(i,11).^2;
    end

    if ((i == length(drop_data)) || (drop_data(i+1,2) ~= drop_data(i,2)))
        dnum_pt = dnum_pt + 1;
        dnum(dnum_pt, 1) = drop_data(i,2); dnum(dnum_pt, 2) = numcount; numcount = 0;

        d10(dnum_pt, 1) = drop_data(i,2); d10(dnum_pt, 2) = tmp_10 / (drop_data(i,3)+1);
        d32(dnum_pt, 1) = drop_data(i,2); d32(dnum_pt, 2) = tmp_30 / tmp_20;
        tmp_10 = 0; tmp_20 = 0; tmp_30 = 0;

        if (abs(drop_data(i,2) - t_f) < 1e-3)
            pt_f = i;
        end
    end
end


% set(groot,'defaultAxesTickLabelInterpreter','latex');
% figure(1); box on; hold on;
% plot(dnum(:,1)-t_i, dnum(:,2), 'linewidth', 2.5);
% xlim([0 t_f-t_i]);
% set(gca, 'FontSize', 14); set(gca,'TickLength', [0.02,0.02]);
% xlabel('$t$', 'interpreter', 'latex');
% ylabel('$N_{drop}$', 'interpreter', 'latex');
% 
% figure(2); box on; hold on;
% plot(d10(:,1)-t_i, d10(:,2), 'linewidth', 2.5);
% plot(d32(:,1)-t_i, d32(:,2), 'linewidth', 2.5);
% xlim([0 t_f-t_i]);
% set(gca, 'FontSize', 14); set(gca,'TickLength', [0.02,0.02]);
% legend('$d_{10}$', '$d_{32}$', 'interpreter', 'latex');
% xlabel('$t$', 'interpreter', 'latex');
% ylabel('$d_{ij}$', 'interpreter', 'latex');



set(groot,'defaultAxesTickLabelInterpreter','latex');



for i=1:length(drop_data)
    if (t_i <= drop_data(i,2) && drop_data(i,2) <= t_f)
        for j=1:N_bins
            if (edges_size(j) <= drop_data(i,11) && drop_data(i,11) < edges_size(j+1) && drop_data(i,6) > L0*0.01 && drop_data(i,6) < L0*0.99)
                N_c_ind(j,ens) = N_c_ind(j,ens) + 1;
                u_ave(j) = u_ave(j) + drop_data(i,12);
                u_ave_ind(j,ens) = u_ave_ind(j,ens) + drop_data(i,12);
                uy_ave(j) = uy_ave(j) + drop_data(i,9);
                ur_ave(j) = ur_ave(j) + sqrt(drop_data(i,8).^2 + drop_data(i,10).^2);
                break;
            end

            if (edges_vel(j) <= drop_data(i,12) && drop_data(i,12) < edges_vel(j+1) && drop_data(i,6) > L0*0.01 && drop_data(i,6) < L0*0.99)
                N_c_ind_vel(j,ens) = N_c_ind_vel(j,ens) + 1;
            end
        end
    end
end

for i=1:length(drop_data)
    if (drop_data(i,2) == t_c)
        for j=1:N_bins
            if (edges_size(j) <= drop_data(i,11) && drop_data(i,11) < edges_size(j+1) && drop_data(i,6) > L0*0.01 && drop_data(i,6) < L0*0.99)
                N_cf(j) = N_cf(j) + 1;
                ur_ave_cf(j) = ur_ave_cf(j) + sqrt(drop_data(i,8).^2 + drop_data(i,10).^2);
                uy_ave_cf(j) = uy_ave_cf(j) + abs(drop_data(i,9));
                ur_var_cf(j) = ur_var_cf(j) + drop_data(i,8).^2 + drop_data(i,10).^2;
                uy_var_cf(j) = uy_var_cf(j) + drop_data(i,9).^2;
                break;
            end
        end
    end
end

for i=1:N_bins
  if (N_c_ind(i,ens) ~= 0) 
      u_ave_ind(i,ens) = u_ave_ind(i,ens) / N_c_ind(i,ens);
  else u_ave_ind(i,ens) = 0;
  end
end

end

for i=1:N_bins
  N_c_ens(i) = sum(N_c_ind(i,:));
  N_c_ens_vel(i) = sum(N_c_ind_vel(i,:));
end



for i=1:N_bins
  for j=1:N_ens
    pdf_c_ind(i,j) = N_c_ind(i,j)./sum(N_c_ind(:,j)) / (edges_size(i+1)-edges_size(i));
  end
end

pdf_c_ind_nonzero = sort(pdf_c_ind, 2, 'descend');
u_ave_ind_nonzero = sort(u_ave_ind, 2, 'descend');


for i=1:N_bins
  bin_top(i) = N_ens;
  for j=1:N_ens
    if (pdf_c_ind_nonzero(i,j) == 0) 
      bin_top(i) = j-1; break;
    end
  end
end



% for i=1:N_bins
%   if (bin_top(i) > 1) 
%     pdf_test_bds(i,:) = bootci(1000, @mean, pdf_c_ind_nonzero(i,1:bin_top(i)));
%     vel_test_bds(i,:) = bootci(1000, @mean, u_ave_ind_nonzero(i,1:bin_top(i)));
%   else
%     pdf_test_bds(i,:) = [NaN NaN];
%     vel_test_bds(i,:) = [NaN NaN];
%   end
% end

for j=1:N_bins
    pdf_c(j) = N_c_ens(j)/sum(N_c_ens)/(edges_size(j+1)-edges_size(j));
    dens_c(j) = N_c_ens(j)/(edges_size(j+1)-edges_size(j))/N_ens;

    pdf_max(j) = max(N_c_ind(j,:)./sum(N_c_ind(:,:)))/(edges_size(j+1)-edges_size(j));
    pdf_min(j) = min(N_c_ind(j,:)./sum(N_c_ind(:,:)))/(edges_size(j+1)-edges_size(j)); if (pdf_min(j) == 0) pdf_min(j) = pdf_c(j); end

    u_max(j) = max(u_ave_ind(j,:));
    u_min(j) = min(u_ave_ind(j,:));

    var_c(j) = var(N_c_ind(j,:)./sum(N_c_ind(j,:))./(edges_size(j+1)-edges_size(j)));
    
    u_ave(j) = u_ave(j) / N_c_ens(j);
    uy_ave(j) = uy_ave(j) / N_c_ens(j);
    ur_ave(j) = ur_ave(j) / N_c_ens(j);

    dens_c_vel(j) = N_c_ens_vel(j)/(edges_vel(j+1)-edges_vel(j))/N_ens;
end

err_c = sqrt(var_c)./sqrt(N_ens);

figure(7); box on; hold on;
pdf = dens_c.*R0./((t_f-t_i)/dt+1)./L0;
scatter(r_binned./R0, pdf, '+', 'LineWidth', 2.5);
% xlim([0.95*del 0.1]);
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');
set(gca, 'FontSize', 16); set(gca,'TickLength', [0.02,0.02]);
xlabel('$r/R_0$', 'interpreter', 'latex');
ylabel('$n(r/R_0)$', 'interpreter', 'latex');

figure(8); box on; hold on;
scatter(r_binned./R0, uy_ave./V0, '+', 'LineWidth', 2.5);
scatter(r_binned./R0, ur_ave./V0, 'x', 'LineWidth', 2.5);

% xlim([0.95*del 0.1]);
set(gca, 'XScale', 'log');
set(gca, 'FontSize', 14); set(gca,'TickLength', [0.02,0.02]);
xlabel('$d/R_0$', 'interpreter', 'latex');
ylabel('$\bar{v}(d/R_0)/U_0$', 'interpreter', 'latex');



for j=1:N_bins
    pdf_cf(j) = N_cf(j)/sum(N_cf)/(edges_size(j+1)-edges_size(j));
    uy_ave_cf(j) = uy_ave_cf(j)/N_cf(j); ur_ave_cf(j) = ur_ave_cf(j)/N_cf(j);

    uy_var_cf(j) = sqrt(uy_var_cf(j)/N_cf(j) - uy_ave_cf(j).^2)/sqrt(N_cf(j)); 
    ur_var_cf(j) = sqrt(ur_var_cf(j)/N_cf(j) - ur_ave_cf(j).^2)/sqrt(N_cf(j)); 
end

% figure(9); box on; hold on;
% scatter(r_binned./0.025, pdf_cf.*0.05, '+', 'LineWidth', 2.5);
% % xlim([0.9*del 2.1]);
% set(gca, 'XScale', 'log');
% set(gca, 'YScale', 'log');
% set(gca, 'FontSize', 16); set(gca,'TickLength', [0.02,0.02]);
% xlabel('$d/R_0$', 'interpreter', 'latex');
% ylabel('$f(d/R_0)$', 'interpreter', 'latex');
% 
% figure(10); box on; hold on;
% errorbar(r_binned, ux_ave_cf, ux_var_cf, '+', 'LineWidth', 2.5);
% % xlim([0.9*del 2.1]);
% set(gca, 'XScale', 'log');
% set(gca, 'FontSize', 14); set(gca,'TickLength', [0.02,0.02]);
% xlabel('$d/R_0$', 'interpreter', 'latex');
% ylabel('$\bar{u}_x(d/R_0)$', 'interpreter', 'latex');
% 
% figure(11); box on; hold on;
% errorbar(r_binned, ur_ave_cf, ur_var_cf, '+', 'LineWidth', 2.5);
% % xlim([0.9*del 2.1]);
% set(gca, 'XScale', 'log');
% set(gca, 'FontSize', 14); set(gca,'TickLength', [0.02,0.02]);
% xlabel('$d/R_0$', 'interpreter', 'latex');
% ylabel('$\bar{u}_r(d/R_0)$', 'interpreter', 'latex');

d_ave_global = d_ave_global / numcount_global
std_err_global = sqrt(std_err_global / numcount_global - d_ave_global^2)/sqrt(numcount_global)

figure(12); set(gca, 'YScale', 'log'); box on;
set(gca, 'FontSize', 16); set(gca,'TickLength', [0.02,0.02]);
hold on; scatter(v_binned./V0, dens_c_vel.*V0./((t_f-t_i)/dt+1), '+', 'LineWidth', 2.5);
xlabel('$u_{\rm frag}/U_0$', 'interpreter', 'latex'); ylabel('$n(u_{\rm frag}/U_0)$', 'interpreter', 'latex');