%% Polynomial fitting %%

clc
close all

%% Data

T = readtable('ExpData/24.3.21.xlsx');

start = 1500;   % Values before this measurement are discarded
Qo = T{start:2150,2};
betaU = T{start:2150,3};
betaIn = 300;
betaIn2 = 420;
Zo = T{start:2150,4};

%% Plotting experimental results

set(groot, 'defaultTextInterpreter','latex');
set(groot, 'defaultAxesTickLabelInterpreter','latex');
set(groot, 'defaultLegendInterpreter','latex');


%% Data 2

T2 = readtable('ExpData/26.3.21.xlsx');

start2 = 810;
Qo2 = T2{start2:1650,2};
betaU2 = T2{start2:1650,3};
Zo2 = T2{start2:1650,4};

%% DataCombined (For Plotting)

QoCom = [Qo; Qo2];
betaUCom = [betaU; betaU2];
ZoCom = [Zo; Zo2];

%% Time (for plotting)

time1 = T{start:2150,6};
time2 = T2{start2:1650,6};
time2 = T2{start2:1650,6}+max(time1)-min(time2);
time = [time1;time2];


%% Plotting experimental results COMBINED

figure('visible','on');

subplot(3,1,1)
plot(time, QoCom, 'b')
ax=gca;
ax.XAxis.Exponent=0;
xlim([min(time), max(time)])
ylabel('$Q_O$ [m$^3$/h]')
% title('Overflow Flowrate')

subplot(3,1,2)
plot(time, betaUCom, 'b')
ax=gca;
ax.XAxis.Exponent=0;
xlim([min(time), max(time)])
ylabel('$\beta_U$ [ppm]')
% title('Underflow Oil Concentration')

subplot(3,1,3)
plot(time, ZoCom, 'b')
ax=gca;
ax.XAxis.Exponent=0;
xlim([min(time), max(time)])
ylim([0 25])
xlabel('$t$ [s]')
ylabel('$Z_O$ [\%]')
% title('Overflow Valve Opening')

% sgtitle('Experimental Data - Polynomial Fitting')

saveas(gcf, [pwd '/Figures/ExperimentalDataPolynomial'],'epsc')

%% Calculating average values

QoAvg1 = mean([Qo(1:67)]);
QoAvg2 = mean([Qo(126:177)]);
QoAvg3 = mean([Qo(197:264)]);
QoAvg4 = mean([Qo(312:431)]);
QoAvg5 = mean([Qo(466:540)]);

QoAvgList = [QoAvg1 QoAvg2 QoAvg3 QoAvg4 QoAvg5];

% Discarded values:
% QoAvg6 = mean([Qo(677:762)]);
% QoAvg7 = mean([Qo(803:850)]);
% QoAvg8 = mean([Qo(893:906)]);
% QoAvg9 = mean([Qo(944:1018)]);
% QoAvg10 = mean([Qo(1044:1068)]);
% 
% QoAvgList2 = [QoAvg6 QoAvg7 QoAvg8 QoAvg9 QoAvg10];

betaUAvg1 = mean([betaU(1:67)]);
betaUAvg2 = mean([betaU(127:175)]);
betaUAvg3 = mean([betaU(214:249)]);
betaUAvg4 = mean([betaU(364:400)]);
betaUAvg5 = mean([betaU(433:535)]);

betaUAvgList = [betaUAvg1 betaUAvg2 betaUAvg3 betaUAvg4 betaUAvg5];

% Discarded values:
% 
% betaUAvg6 = mean([betaU(677:762)]);
% betaUAvg7 = mean([betaU(803:850)]);
% betaUAvg8 = mean([betaU(893:906)]);
% betaUAvg9 = mean([betaU(944:1018)]);
% betaUAvg10 = mean([betaU(1044:1068)]);
% 
% betaUAvgList2 = [betaUAvg6 betaUAvg7 betaUAvg8 betaUAvg9 betaUAvg10];

%% Calculating average values - Dataset 2

QoAvg11 = mean([Qo2(1:100)]);
QoAvg12 = mean([Qo2(110:257)]);
QoAvg13 = mean([Qo2(285:440)]);
QoAvg14 = mean([Qo2(468:662)]);
QoAvg15 = mean([Qo2(690:837)]);

QoAvgListData2 = [QoAvg11 QoAvg12 QoAvg13 QoAvg14 QoAvg15];

betaUAvg11 = mean([betaU2(1:100)]);
betaUAvg12 = mean([betaU2(110:257)]);
betaUAvg13 = mean([betaU2(285:440)]);
betaUAvg14 = mean([betaU2(468:662)]);
betaUAvg15 = mean([betaU2(690:837)]);

betaUListData2 = [betaUAvg11 betaUAvg12 betaUAvg13 betaUAvg14 betaUAvg15];

%% Calculating separation efficiency

etaList = [];
for i = 1:length(betaUAvgList)
    eta = 1 - (betaUAvgList(i)/betaIn);
    etaList = [etaList, eta];
end

% etaList2 = [];
% for i = 1:length(betaUAvgList2)
%     eta2 = 1 - (betaUAvgList2(i)/betaIn2);
%     etaList2 = [etaList2, eta2];
% end


etaListData2 = [];
for i = 1:length(betaUListData2)
    eta3 = 1 - (betaUListData2(i)/betaIn);  % Betain is 300 for this aswell
    etaListData2 = [etaListData2, eta3];
end


%% Polynomial fitting


QoAvgList3 = sort([QoAvgList QoAvgListData2]);

etaList3 = sort([etaList etaListData2]);

p2 = polyfit(QoAvgList3, etaList3, 2);  % 2nd degree polynomial

xlist = linspace(min(QoAvgList3), max(QoAvgList3), 10000);

pol2 = polyval(p2,xlist);

p3 = polyfit(QoAvgList3, etaList3, 3);  % 3rd degree polynomial

pol3 = polyval(p3,xlist);   % 

%% Plotting

figure
scatter(QoAvgList3, etaList3, 'filled', 'r')
hold on
plot(xlist, pol2, 'b')

xlabel('$Q_O$ [m$^3$/h]')
ylabel('$\eta$ [-]')
legend('Exp data','Polynomial')
legend('Location', 'northwest')
% title('Polynomial Expression')
saveas(gcf, [pwd '/Figures/PolynomialApproximation'],'epsc')