%% Extract data
clear
clc
close all
load('Actcons.mat','ActCon')
%load('AC.mat','ActCon')
X = [];
Y = [];
Z = [];
A = [];
C = [];
c=[];
m = 1;
l = ones(1,6);
w = 1;
ceq=[];
FsX =[]; 
FsY =[]; 
for i= 1:1:size(ActCon,2)
    if ActCon(i).Flag == 1 
        X(m) = ActCon(i).QEnv;
        Y(m) = ActCon(i).QBOG;
        Z(m) = ActCon(i).Objective;
        A{m} = ActCon(i).Act;
        I(m) = ActCon(i).T- 273.15;
        J(m) = ActCon(i).MBOG;
        P1(m) = ActCon(i).HighP*20;
        P2(m) = ActCon(i).LowP*3;
        F(m) = ActCon(i).Flow;
        Sf(m) = ActCon(i).SplitFlow*65;
        T(m) = ActCon(i).T3 - 273.15;
        [~,~,ceq(i,:)]=MainEO(ActCon(i).w,(ActCon(i).T-273.15+32),...
            ActCon(m).MBOG/610,ActCon(m).QEnv/80.5852);
        % Active constraints
        a = A{m};
        if a == [0] % No active constraints
            n = 1;
            QE_{n}(l(n)) = X(m);
            QBOG_{n}(l(n)) = Y(m);
            c(n,:) = [1 0 0]; % Red
            l(n)=l(n)+1;
        end
        if size(a,2) == 2
            if a == [0,35] % Lower bound High P
                n = 6;
                QE_{n}(l(n)) = X(m);
                QBOG_{n}(l(n)) = Y(m);
                c(n,:)=[0 0 1];
                l(n)=l(n)+1;% Blue
            elseif a == [0,36] % Lower bound Active Charge
                n = 3;
                QE_{n}(l(n)) = X(m);
                QBOG_{n}(l(n)) = Y(m);
                c(n,:)=[0 1 0]; % Green
                l(n)=l(n)+1;
            elseif a == [0,39] % Upper bound High P
                n = 5;
                QE_{n}(l(n)) = X(m);
                QBOG_{n}(l(n)) = Y(m);
                c(n,:) = [0 1 1];% Yellow

                l(n)=l(n)+1;

            end
        elseif size(a,2) == 3
            if a == [0,36,40] % Lower bound Active Charge  Lower bound High P
                n = 2;
                QE_{n}(l(n)) = X(m);
                QBOG_{n}(l(n)) = Y(m);
                c(n,:) = [1 0 0]; % Red
                l(n)=l(n)+1;
            elseif a == [0,36,39] % Lower bound Active Charge % Upper bound High P
                n = 4;
                QE_{n}(l(n)) = X(m);
                QBOG_{n}(l(n)) = Y(m);
                c(n,:) = [1 0 1]; % Magenta

                l(n)=l(n)+1;
            end
        end
        m = m+1;
    else
        FsX(w) = ActCon(i).QEnv;
        FsY(w) = ActCon(i).QBOG;
        w = w+1;
    end
end


%% Plot 1
np = 20;
plotres( X, Y, Z, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'W_{Compressor}[kj/s]',1,'Wc'  )

%% Plot 2

ilin = linspace(min(I),max(I),np);
jlin = linspace(min(J),max(J),np);
[i, j] = meshgrid(ilin,jlin);

k = griddata(I',J', Y', i,j,'cubic');

figure(2)

axis tight
hold on
grid on

xlabel('Temperature_{BOG,In}[^\circC]', 'FontSize',12)
ylabel('Flow_{BOG}[kmol/s]', 'FontSize',12)
zlabel('Q_{BOG}[kJ/s]')

contour(i,j,k, 'ShowText','on')
clabel([] ,[] ,'FontSize',12)
printfig(figure(2),'Qbog')


%% Plot 3

figure(3)
axis tight;
hold on
grid on
xlabel('Q_{Env}[kJ/s]', 'FontSize',12)
ylabel('Q_{BOG}[kJ/s]', 'FontSize',12)
% Active constraints
for n = 1:size(QE_,2)
    scatter(QE_{n}',QBOG_{n}',60,c(n,:),'filled')
end
scatter(FsX',FsY',60,[0 0 0],'filled')
legend('No active constraints', 'Failed states','Location', 'best', 'FontSize',12)
legend('show') 
printfig(figure(3),'Actcon')

%% Plot 4
plotres( X, Y, P2, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'P_{Low}[bar]',4,'LowP'  )
%% Plot 5
plotres( X, Y, P1, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'P_{High}[bar]',5,'HighP'  )
%% Plot 6
plotres( X, Y, F, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'Refrigerant flow (dimless)',6,'Flow'  )
%% Plot 7
plotres( X, Y, T, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'T_{3}[^\circ C]',7,'Temp3'  )
%% Plot 8
plotres( X, Y, Sf, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'Split fraction',8,'Splitf'  )

%% plot 
Eff = (P1)./(P2);
p = [0.0015   -0.0314    0.1911    0.4217]; %3rd degree polynomial. 
Eff = polyval(p,Eff);

plotres( X, Y, Eff, np,'Q_{Env}[kJ/s]' , 'Q_{BOG}[kJ/s]', 'Compressor efficiency',9,'Compeff'  )

