% Determine runout distance for a range of Temperatures and different
% conditions (no sedimentation, no particle entrainment, hot entrainment,
% cold entrainment)
% Created by Kristen Fauria 
% February 2016


clear all

%%% Preliminaries %%%

CA=1.005;%heat capacity, (kJ/(kg K))
CR=0.84;%heat capacity of basalt, (kJ/(kg K));
B=3.43*10^-3;%expansion coefficient, air (1/K)
R=8.314;%specific gas constant, (J/(mol K)); (get rid of 10^3?
P=10^5;%atmospheric pressure (Pa)
M=1/(28.97/1000); %molar mass of air (kg/mol) (need to invert?) - yes such that M = 1/(28.97/1000)
E = 0.1; %Entrainment rate
G = 9.8; % gravity m/s2
Fr=1.2; %Froude number
DT=5; %time step in seconds

%%% Initial conditions for all currents %%%

% Temperatures
TA=20+273; % amibent temperature = 20 degrees C.
RHOA = P/(M*R*TA); % density of the ambient air kg/m3

% Mass of particles in different size fractions %
MP1 = [100 10]; %Change this to change the initial mass of particles throughout the models (kg)
D = [0.0001 0.01]; % diameter of each particle size fraction (m)
RHOP = [2400 1000]; % density of each particle size fraction (kg/m^3)
PVEL = [0.1 10]; %prescribed particle sedimentatin velocity m/s
prescribe = 0; % prescribe = 0 allows settling velocities to be calculated via Dufek et al. 2009 code, prescribe = 1 (or anything else transfers the values on the previous line)


% Initial mass of air in the current %
MA1=30; %initial air mass for all models

%%% Temperature range %%%

TR1=380+273; %initial particle temperature (K) , 300
TR2=1100+273; %initial particle temperature (K) ,1200
Temperature=TR1:20:TR2;
Runout=zeros(length(Temperature),5);
Runout(:,1)=Temperature;
%%
for ind=1:length(Temperature); % initial T of rocks in K
    TR=Temperature(ind);
    %save initial temperature
    for type=1:4
        [RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL,i]=PDC_1D(type,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe);
        Runout(ind,type+1)=RO; % save runout distance
        TempC(ind,type)=TM(1);%initial temperature of the mixture
    end
end

%%
% Plot runout distance
%close all
figure;
hold on
set(gca,'FontSize',20);
load('colorblind_colormap.mat');
F=[colorblind(5,:);colorblind(4,:);colorblind(1,:);colorblind(6,:)];

for type = 1:4
plot(TempC(:,type)-273,Runout(:,type+1),'Color',F(type,:),'LineWidth',3);%no settling or particle entrainment
end
xlabel('Initial current temperature (^{\circ}C)');
ylabel('Runout distance (m)');
ylim([0 3000]);
xlim([250 900]);

legend('No sedimentation','Sedimentation only', 'Cold particle entrainment','Hot particle entrainment');