clear all, clc, clf, close all
% ------------- Material constants -------------- %
rho = 998; g = 9.81;
patm = 101325; dT = 2.5e-4;

pScale = 10^-3;
% ------------- Constant settings --------------- %
ln = 1; lnn = 1.5; fs = 14; fsL = 16;

% ------------- Read surface data --------------- %
Q = abs(readmatrix('postProcessing/flowRateInlet/0/surfaceFieldValue.dat', 'CommentStyle','#')*1e3);
time = Q(:,1)/1e3; dT = time(2)-time(1);
Q = Q;
% ------------- Read probe data ----------------- %
pProbes = readmatrix('postProcessing/probes/0/p', 'CommentStyle','#')*rho*pScale;

% --------- Read Numerical Wang data ------------ %
folder= [pwd '/validationData/'];
wangData = readmatrix([folder 'wangEtAl_numericalData.csv']);
wangData(:,1) = wangData(:,1) + 3;
wangData(:,3) = wangData(:,3)*rho*10*pScale;

expData = readmatrix([folder 'petitEtAl_experimentalData.txt'],'CommentStyle','#');
expData(:,2) = expData(:,2)*pScale;

% ----------------- Plot Data ------------------- %
figure(1)
	p(1) = plot(wangData(:,1), wangData(:,2)*10^3,'--r','linewidth',lnn); hold on
	p(2) = plot(time,Q(:,2),'k','linewidth',lnn);
	grid on
	set(gca, "linewidth", 1, "fontsize", fs)
	xlabel('Time (s)','fontsize', fsL,'Interpreter','latex')
	ylabel('Flow rate (l/s)','fontsize', fsL,'Interpreter','latex')
	xlim([0 30]); ylim([-20 60])
	xticks(0:5:30);
	legend(flip(p), 'headLossPressure $\quad$' ,'Wang et al. (2017)');
	set (legend, "fontsize", 12,"location", "northoutside", "orientation", "horizontal",'Interpreter','latex');	
	legend boxoff	
	ax = gca;
	set(ax,'ticklabelinterpreter','latex');

	% create a new pair of axes inside current figure
	axes('position',[.18 .65 .15 .20])
	box on % put box around new pair of axes
	plot(wangData(:,1), wangData(:,2)*10^3,'--r','linewidth',lnn); hold on % plot on new axes
	plot(time,Q(:,2),'k','linewidth',lnn); % plot on new axes
	xlim([9 19]); ylim([0 5]);
	axBox1 = gca;
	set(axBox1,'position',[.29 .17 .45 .1])
	set(axBox1,'ticklabelinterpreter','latex');
	set(gcf,'Position',[2700 100 600 450])
	
figure(2)
	p(1) = plot(wangData(:,1), wangData(:,3),'--r','linewidth',lnn); hold on	
	p(3) = plot(expData(:,1), expData(:,2),'-.c','linewidth',lnn);
	p(2) = plot(time,pProbes(:,2),'k','linewidth',lnn);
	grid on
	set(gca, "linewidth", 1, "fontsize", fs)
	xlabel('Time (s)','fontsize', fsL,'Interpreter','latex')
	ylabel('Pressure (kPa)','fontsize', fsL,'Interpreter','latex')
	legend([p(2) p(1) p(3)], 'headLossPressure' ,'Wang et al. (2017)', 'Petit et al. (2015)');
	set (legend, "fontsize", 12,"location", "north", "orientation", "vertical",'Interpreter','latex');	
	xlim([0 30]);
	ylim([0 50]);
	xticks(0:5:30);
	ax = gca;
	set(ax,'ticklabelinterpreter','latex');
	
	% create a new pair of axes inside current figure
	axes('position',[.18 .65 .15 .20])
	box on % put box around new pair of axes
	plot(wangData(:,1), wangData(:,3),'--r','linewidth',lnn); hold on % plot on new axes
	plot(time,pProbes(:,2),'k','linewidth',lnn); % plot on new axes
	plot(expData(:,1), expData(:,2),'-.c','linewidth',lnn)
	xlim([9.5 10.5]); ylim([29 31]);
	axBox1 = gca;
	set(axBox1,'position',[.17 .6 .15 .20])
	set(axBox1,'ticklabelinterpreter','latex');
	
	% create a new pair of axes inside current figure
	axes('position',[.18 .65 .15 .20])
	box on % put box around new pair of axes
	plot(wangData(:,1), wangData(:,3),'--r','linewidth',lnn); hold on % plot on new axes
	plot(time,pProbes(:,2),'k','linewidth',lnn); % plot on new axes
	plot(expData(:,1), expData(:,2),'-.c','linewidth',lnn)
	xlim([17.5 18.5]); ylim([32 34]);
	axBox2 = gca;
	set(axBox2,'position',[.72 .6 .15 .20])
	set(axBox2,'ticklabelinterpreter','latex');
	set(gcf,'Position',[2700 100 600 450])
	
if ~exist('Figures', 'dir')
   mkdir('Figures')
end

[~,folderName,~]=fileparts(pwd);
saveas (1, ['Figures/Q_' folderName '.eps'],'epsc')
saveas (2, ['Figures/p_' folderName '.eps'],'epsc')

