% Plot different current variables with distance for specific initial
% conditions
% Created by Kristen Fauria
% May 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));
P=10^5;%atmospheric pressure (Pa)
M=1/(28.97/1000); %molar mass of air (kg/mol)
E = 0.1; %Entrainment rate
G = 9.8; % gravity m/s2
Fr=1.2; %Froude number
DT=0.1; %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) (fine volcanic ash < 0.063 mm)
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

%%% Initial temperature of the rocks in the mixture %%%
TR = 600 + 273;


% Current types (1 = no settling, 2 = no entraining particles, 3 = cold
% entrainment, 4 = hot entrainment)
% plot variables vs distance


% plot just density and temperature

load('colorblind_colormap.mat');

figure; hold on
F=[colorblind(5,:);colorblind(4,:);colorblind(1,:);colorblind(6,:)];
set(gca,'FontSize',20);
for type = 1:4
[RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL,i,PSAVE]=PDC_1D(type,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe);
%print
disp(PSAVE);
subplot(2,1,1); set(gca,'FontSize',20); hold on
semilogx(X,RHOM./RHOA,'Color',F(type,:),'LineWidth',3); 
ylabel('density/ambient density');
xlabel('distance (m)');


subplot(2,1,2); set(gca,'FontSize',20); hold on
semilogx(X,TM-273,'Color',F(type,:),'LineWidth',3);
ylabel('temperature (^{\circ}C)');
xlabel('distance (m)');
end

legend('No sedimentation','Sedimentation only', 'cold particle entrainment','hot particle entrainment');
%%

figure; hold on
F=['k',[0.5 0.5 0.5],'b','r'];
set(gca,'FontSize',20);
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);
subplot(3,2,1); set(gca,'FontSize',20); hold on
semilogx(X,RHOM./RHOA,F(type),'LineWidth',2); 
ylabel('density/ambient density');
xlabel('distance (m)');
subplot(3,2,2); set(gca,'FontSize',20); hold on
semilogx(X,H,F(type),'LineWidth',2);
ylabel('height (m)');
xlabel('distance (m)');
subplot(3,2,3); set(gca,'FontSize',20); hold on
semilogx(X,MASS,F(type),'LineWidth',2);
ylabel('mass (kg)');
xlabel('distance (m)');
subplot(3,2,4); set(gca,'FontSize',20); hold on
semilogx(X,TM-273,F(type),'LineWidth',2);
ylabel('temperature (C)');
xlabel('distance (m)');
subplot(3,2,5); set(gca,'FontSize',20); hold on
semilogx(X,U,F(type),'LineWidth',2);
ylabel('velocity (m s^-1)');
xlabel('distance (m)');
subplot(3,2,6); set(gca,'FontSize',20); hold on
semilogx(X,VRVT,F(type),'LineWidth',2);
ylabel('Rock volume fraction');
xlabel('distance (m)');
end

%%
% plot just density and temperature

figure; hold on
F=['k',[0.5 0.5 0.5],'b','r'];
set(gca,'FontSize',20);
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);

subplot(2,1,1); set(gca,'FontSize',20); hold on
semilogx(X,RHOM./RHOA,F(type),'LineWidth',2); 
ylabel('density/ambient density');
xlabel('distance (m)');

subplot(2,1,2); set(gca,'FontSize',20); hold on
semilogx(X,TM-273,F(type),'LineWidth',2);
ylabel('temperature (C)');
xlabel('distance (m)');
end

% 
% hold on
% [RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL]=PDC_1D(2,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe);
% subplot(3,2,1);hold on
% plot(X,RHOM,'.g');
% ylabel('density');
% xlabel('distance (m)');
% subplot(3,2,2);hold on
% plot(X,H,'.g');
% ylabel('height');
% xlabel('distance (m)');
% subplot(3,2,3);hold on
% plot(X,MASS,'.g');
% ylabel('mass');
% xlabel('distance (m)');
% subplot(3,2,4);hold on
% plot(X,TM,'.g');
% ylabel('temperature (K');
% xlabel('distance (m)');
% subplot(3,2,5);hold on
% plot(X,U,'.g');
% ylabel('velocity (m/s)');
% xlabel('distance (m)');
% subplot(3,2,6);
% plot(X,VRVT,'.g');
% ylabel('Rock volume fraction');
% xlabel('distance (m)');
% 
% 
% hold on
% [RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL]=PDC_1D(3,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe);
% subplot(3,2,1);hold on
% plot(X,RHOM,'.b');
% ylabel('density');
% xlabel('distance (m)');
% subplot(3,2,2);hold on
% plot(X,H,'.b');
% ylabel('height');
% xlabel('distance (m)');
% subplot(3,2,3);hold on
% plot(X,MASS,'.b');
% ylabel('mass');
% xlabel('distance (m)');
% subplot(3,2,4);hold on
% plot(X,TM,'.b');
% ylabel('temperature (K');
% xlabel('distance (m)');
% subplot(3,2,5);hold on
% plot(X,U,'.b');
% ylabel('velocity (m/s)');
% xlabel('distance (m)');
% subplot(3,2,6);
% plot(X,VRVT,'.b');
% ylabel('Rock volume fraction');
% xlabel('distance (m)');
% 
% hold on
% [RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL]=PDC_1D(4,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe);
% subplot(3,2,1);hold on
% plot(X,RHOM,'.r');
% ylabel('density');
% xlabel('distance (m)');
% subplot(3,2,2);hold on
% plot(X,H,'.r');
% ylabel('height');
% xlabel('distance (m)');
% subplot(3,2,3);hold on
% plot(X,MASS,'.r');
% ylabel('mass');
% xlabel('distance (m)');
% subplot(3,2,4);hold on
% plot(X,TM,'.r');
% ylabel('temperature (K');
% xlabel('distance (m)');
% subplot(3,2,5);hold on
% plot(X,U,'.r');
% ylabel('velocity (m/s)');
% xlabel('distance (m)');
% subplot(3,2,6);
% plot(X,VRVT,'.r');
% ylabel('Rock volume fraction');
% xlabel('distance (m)');

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