% This script Finds the restitution coefficient in the Splash Function for
% the 8 different studies

% Created by Kristen Fauria: May 15, 2016
% We use a Splash Function with crater depth as the characteristic length
% scale

close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data must be loaded from SplashCompilation.xlsx (tab 3 - Matlab) as
% "data"

%Preliminaries %
g=9.8; % gravity (m/s2)
mew=0.5; %assumed tangent of the friction angle

% Remove some data with nan values
[row,col]=find(isnan(data(:,6)));%Remove rows with nan values for the number of suspended particles
data(row,:)=[];
[row,col]=find(isnan(data(:,3)));%Remove rows with nan values for the ratio of particle masses
data(row,:)=[];
%%
% Label each column from the spreadsheet

Study=data(:,1); %The study number where 1 = Mitha et al., 2 = Beladjine et al.,; 3 = Oger et al.,; 4 = Wu et al., (2013); 5 = This study; 6 = Werner and Haff; 7 = Anderson and Haff (1988); 8 = Willetts and Rice.
Method=data(:,2); % 1 = physical experiments; 2 = computer simulations
MassRatio=data(:,3); % Ratio of incident particle mass to ejected particle mass
VI=data(:,4); % Vertical incident speed in meters per second
VR=data(:,5); % Vertical rebound speed in meters per second
NumberEjected=data(:,6); % The total number of ejected particles
DE=data(:,7); % The diameter of ejected particles in cm
ImpactAngle=data(:,8); % The degrees from horizontal at which the impactor hit the substrate
IncidentSpeed=data(:,9); % The total incident speed of the impactor in meters per second
EjectionAngle=data(:,10); % The angle at which the impactor rebounded from the substrate
ReboundSpeed=data(:,11); % The speed at which the impactor rebounded from the substrate
AngleAssumption=data(:,12); % 0 = impacting particle ejection angle measured; 1 = imacting particle ejection angle assumed the same as the incident angle 
DI=data(:,13); % The incoming particle diamter in cm
DensityRatio=data(:,14); % Ratio of incident particle density to substrate or ejected particle density 


%%

%%%% Figure 3 in the manuscript %%%%
%%% First just plot our experimental results m_e/m_i vs. 1/2v2/(g*de) %%%


figure;
set(gca,'FontSize',20)
load('colorblind_colormap.mat');
MarkerType=['*','v','.','x','o','<','s','d','^','+'];
order=[1,8,6,7,2,3,4,5,9,10];
set(gca,'FontSize',20)
for i=8:10%number of studies
    k=order(i);
    hold on
    for j=1:length(data(:,1))
        if Study(j)==k% give each study a different symbol and plot in the order given by "order"
            hold on
            loglog(IncidentSpeed(j)^2/(g*DE(j)./100),NumberEjected(j)./MassRatio(j),MarkerType(k),'Color','k','MarkerSize',10,'LineWidth',2);       
        end   
    end     
end
xlabel('$\frac{V_{i}^{2}}{2gd_{e}}$','interpreter','latex','FontSize',30);
ylabel('$\frac{mass~ejected}{incident~particle~mass}$','interpreter','latex','FontSize',30);
%%
%%% Figure 4 in the manuscript %%%
%%% Plot to find the energy loss coefficient %%%
H=IncidentSpeed.^2/(2*g);% H in meters (the corresponding height from which an impactor is dropped)
DC=0.14*mew^-1*(DensityRatio).^(0.5).*(DI/100).^(2/3).*(H).^(1/3); % The depth of a crater made from the impact in meters
EL=(2*g.*NumberEjected.*DC)./(MassRatio.*IncidentSpeed.^2);

% Determine the average Energy loss coefficient for each study and the
% standard deviation
for i=1:max(data(:,1))%number of studies
    Temp=EL(Study==i);%Grab data from each study
    EL2(i,1)=mean(Temp);
    EL2(i,2)=std(Temp);
    clear Temp
end

% This is the order in which the studies will be plotted
order=[1,8,6,7,2,3,4,5,9,10];

figure;
load('colorblind_colormap.mat');
MarkerType=['*','v','.','x','o','<','s','d','^','+'];
Label={'Mitha et al. (1986)', 'Willetts and Rice (1986)','Werner and Haff (1986)','Anderson and Haff (1988)', 'Beladjine et al. (2007)', 'Oger et al. (2008)','Wu (2013)','This study: pumice and sand','This study: nylon and sand','This study:wood and sand'};

set(gca,'FontSize',20)
for i=1:max(data(:,1))%number of studies
    k=order(i);
    hold on
    for j=1:length(data(:,1))
        if Method(j)==1 %define color based on experimental or simulation study
            Color=colorblind(5,:);
        else
            Color=colorblind(4,:);
        end  
        if Study(j)==k % give each study a different symbol and plot in the order given by "order"
            hold on
            semilogy(i,EL(j),MarkerType(k),'Color',Color,'MarkerSize',15,'LineWidth',2);
            hold on     
            errorbar(i,EL2(k,1),EL2(k,2),'Color',colorblind(6,:),'Linewidth',2);hold on
        end
        
    end
       
end
xlabel('Study number');
%set(gca,'XTick',1:10,'XTickLabel',Label);
xlim([0 11]);
ylim([0 0.05]);
ylabel('$e_{n}^2=\frac{2N_{e}m_{e}gd_{c}}{m_{i}V_{i}^{2}}$','interpreter','latex','FontSize',30);

% Now label the figur with the mean restitution coefficient and standard
% deviation

for i=1:10
text(i,0.04,num2str(EL2(order(i),1),1));
text(i,0.03,num2str(EL2(order(i),2),1));
end

%%
%%% Figure 5 in the manuscript %%%

%Plot incoming kinetic energy vs. potential energy with derived En

figure;
MarkerType=['*','v','.','x','o','<','s','d','^','+'];
% Plot a 1:1 line
set(gca,'FontSize',20)
x=[10^-4 2 1000000];
y=[10^-4 2 1000000];
plot(x,y,'k-','LineWidth',1.5);hold on

for i=1:max(data(:,1))%number of studies
    hold on
    for j=1:length(data(:,1))
        if Method(j)==1 %define color based on experimental or simulation study
            Color=colorblind(5,:);
        else
            Color=colorblind(4,:);
        end  
        if Study(j)==i % give each study a different symbol
            hold on
            plot(EL2(i,1)*MassRatio(j)*IncidentSpeed(j)^2,2*NumberEjected(j)*g*DC(j),MarkerType(i),'Color',Color,'MarkerSize',8,'LineWidth',1.5);
            hold on     
            
        end
        
    end
       
end

xlabel('$\frac{e_{n}^2m_{i}V_{i}^{2}}{m_{e}}$','interpreter','latex','FontSize',30);
ylabel('$2N_{e}gd_{c}$','interpreter','latex','FontSize',25);
xlim([0.0001,10^6]);
ylim([0.0001,10^6]);