load All_data.mat

% Removing dummy values in points with no convergence
for i = 1:9
    [u,v] = find(X(:,:,i) == -1);
    for j = 1:length(u)
        X(u(j),v(j),i) = NaN;
    end
end

% Plotting
for i = 1:9
    figure(i)
    set(gca,'FontSize',14)
    X(:,:,i) = X(:,:,i)*scx(i);
    surf(Tci_ref,w_V_ref.*3.6,X(:,:,i),'EdgeColor','none')
    xlabel('\itd\rm_1 [°C]')
    ylabel('\itd\rm_2 [t/h]')
end

% Adding z labels
figure(1)
zlabel('\itu\rm_1')
axis([65 90 500 900 0 1])
figure(2)
zlabel('\itu\rm_2')
axis([65 90 500 900 0 1])
figure(3)
zlabel('\itu\rm_3')
axis([65 90 500 900 0 1])
figure(4)
zlabel('\itu\rm_5')
axis([65 90 500 900 0 1])
figure(5)
zlabel('\itu\rm_4')
axis([65 90 500 900 0 1])
figure(6)
zlabel('\itu\rm_6 [W]')
axis([65 90 500 900 5000 10000])
figure(7)
zlabel('\itu\rm_8 [W]')
axis([65 90 500 900 0 1000])
figure(9)
zlabel('\itu\rm_7 [W]')
axis([65 90 500 900 0 20000])

load All_data

% Invoking sub-plot for calculation of additional variables
zzplt = -1*ones(length(w_V_ref),length(Tci_ref),4);
for i = 1:length(w_V_ref)
    for j = 1:length(Tci_ref)
        for k = 1:4
            zzplt(i,j,k) = Brobekk_calc(X(i,j,:),Prob,w_V_ref(i),Tci_ref(j),k);
        end
    end
end

% Removing dummy values in points with no convergence
[u,v] = find(zzplt(:,:,1) < 88);
for j = 1:length(u)
    zzplt(u(j),v(j),1) = NaN;
end
[u,v] = find(zzplt(:,:,2) == 0);
for j = 1:length(u)
    zzplt(u(j),v(j),2) = NaN;
end
[u,v] = find(zzplt(:,:,3) == 0);
for j = 1:length(u)
    zzplt(u(j),v(j),3) = NaN;
end
[u,v] = find(zzplt(:,:,4) == 0);
for j = 1:length(u)
    zzplt(u(j),v(j),4) = NaN;
end
[u,v] = find(cst(:,:) == 0);
for j = 1:length(u)
    cst(u(j),v(j)) = NaN;
end

% Additional plots
figure(o-1)
set(gca,'FontSize',14)
surf(Tci_ref,w_V_ref.*3.6,zzplt(:,:,1),'EdgeColor','none')
xlabel('\itd\rm_1 [°C]')
ylabel('\itd\rm_2 [t/h]')
zlabel('\ity\rm_5 [°C]')
axis([65 90 500 900 90 140])

figure(o)
set(gca,'FontSize',14)
surf(Tci_ref,w_V_ref.*3.6,zzplt(:,:,2),'EdgeColor','none')
xlabel('\itd\rm_1 [°C]')
ylabel('\itd\rm_2 [t/h]')
zlabel('\ity\rm_2 [°C]')
axis([65 90 500 900 110 130])

figure(o+1)
set(gca,'FontSize',14)
surf(Tci_ref,w_V_ref.*3.6,zzplt(:,:,3),'EdgeColor','none')
xlabel('\itd\rm_1 [°C]')
ylabel('\itd\rm_2 [t/h]')
zlabel('\ity\rm_3 [°C]')
axis([65 90 500 900 90 150])

figure(o+2)
set(gca,'FontSize',14)
surf(Tci_ref,w_V_ref.*3.6,cst,'EdgeColor','none')
xlabel('\itd\rm_1 [°C]')
ylabel('\itd\rm_2 [t/h]')
zlabel('Cost')
axis([65 90 500 900 0 50000])