% This script compares measurements from eight different studies to six
% different Splash Functions

% 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 
%%
% compare our results to existing Splah Functions

% Create figure 
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)
x=[10^-3 10 10^6];
y=[10^-3 10 10^6];
plot(x,y,'k-','LineWidth',1.5);hold on
en=sqrt(0.008);

% Calculate expected crater depth
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);

%EL=(2*g.*NumberEjected.*DC)./(MassRatio.*IncidentSpeed.^2);

XLabels={'$3.36\sin{\theta_i}(0.00572V_{i}-0.915)$';'$0.55(V_i-40\sqrt(gd_i))$';'$(0.4V_i-4)+(V_i-25)\sin{\theta_i}$';'$[\frac{d_{i}}{d_{e}}]^{3}\frac{V_{i}^2}{gd_{i}}$';'$0.3\frac{V_{i}^{2}}{a^{2}{gd}}$';'$\frac{\beta_2{V_{in}}^{2}}{{gd}}$';'$\frac{\frac{1}{2}e_{n}^{2}m_{i}V_{i}^{2}}{0.14m_{e}g\mu^{-1}\sqrt{\frac{\rho_{i}}{\rho_{s}}}d_{i}^{\frac{2}{3}}\big(\frac{V_{i}^{2}}{2g}\big)^{\frac{1}{3}}}$'};


% Calculate the Splash Function for each study
for i=1:10%number of studies
    k=order(i);
    hold on
    for j=1:length(data(:,1))
         % give each study a different symbol and plot in the order given by "order"
            hold on
            if Method(j)==1 %define color based on experimental or simulation study
            Color=colorblind(5,:);
            else
            Color=colorblind(4,:);
            end  
            
            Splash(j,1)= 3.36*sin(deg2rad(ImpactAngle(j)))*(0.00572*(IncidentSpeed(j)-0.915)); %Werner (1990)
            Splash(j,2)= 0.55*(IncidentSpeed(j)-40*sqrt(g*DI(j)/100)); %Oger (2008), first
            Splash(j,3)= (0.4*IncidentSpeed(j)-4)+(IncidentSpeed(j)-25)*sin(ImpactAngle(j)); % Oger (2008) second
            Splash(j,4)= (DI(j)/DE(j))^3*(IncidentSpeed(j)^2/(g*DE(j)/100)); % Anderson (1987)
            Splash(j,5)= 0.3*(IncidentSpeed(j)^2/(g*DI(j)/100)); %Androetti (2004)
            Splash(j,6)= (en)^2*(VI(j)^2/(2*g*DI(j)/100)); %Wu (2013)
            Splash(j,7)= (en^2)*MassRatio(j)*IncidentSpeed(j)^2/(2*g*DC(j)); %our splash function with constant en
            
            for ij =1:7
            set(gca,'FontSize',16);
            subplot(3,3,ij)
            xlim([10^-2 10^6]);
            ylim([10^2 10^6]);
            hold on
            if Study(j)==k;
            loglog(Splash(j,ij),NumberEjected(j),MarkerType(k),'MarkerEdgeColor',Color,'MarkerSize',5,'LineWidth',2);
            end
            plot(x,y,'k','LineWidth',1.5);
            ylabel('N_e');
            ax.XTick = [0.01 10^0 10^2 10^4 10^6];
            ax.YTick = [0.01 10^0 10^2 10^4 10^6];
            xlabel(XLabels{ij},'interpreter','latex','FontSize',20);%label x-axis
            end
    end
end