
% This function is used for plotting the values obtained which where saved.
% It plots the temperature of the first and final stage as well as the 
% concentration of the final stage.

% Definition of the current path folder
FileName = mfilename('fullpath');
[directory,~,~] = fileparts(FileName);
[dirmain,~,~] = fileparts(directory);

% Import of CASADI
addpath([dirmain '\Util1. Casadi3.1.0'])
addpath([dirmain '\Util2. Functions'])
import casadi.*

% Clearing of the workspace and closing of all windows
% clear;
clc; close all;

% Loading of the desired file
% load('eNMPC_DistRej_pQ_V2')
load([dirmain '\Data\SS_opt.mat'])

% Clearance of the figure
clear fig 

% Calculation of the inlet temperature to bed one
if size(sol_Dyn.z,1) < 31
    mSHEx  = (1-sum(sol_Dyn.u(:,[1 1:end]))).*in_int(1,:);
    CStar  = mSHEx./in_int(1,:);
    NTU    = 151688./(mSHEx*3500);
    eps    = (1-exp(-NTU.*(1-CStar)))./(1-CStar.*exp(-NTU.*(1-CStar)));
    TexHEx = (eps.*sol_Dyn.x(end,:)+(1-eps).*in_int(3,:));

    mS1   = sol_Dyn.u(1,[1 1:end]).*in_int(1,:);
    m1    = mS1+mSHEx;
    Tin   = (mSHEx.*TexHEx + mS1.*in_int(3,:))./m1;
else
    Tin   = sol_Dyn.z(31,:);
    inl(1) = 252/3.6;         % Feed flow (kg/s)
    inl(2) = 200;             % Pressure [bar]
    inl(3) = 250;             % Temperature [C]
    inl(4) = 0.08;            % Mass fraction ammonia [-]
    inl(5) = 3;               % Hydrogen/nitrogen molar ratio [-]
    inl(6) = 1;               % Catalyst activity 
    in_int = in_int+inl'*ones(1,t.N+1);
end

% Calculation of the step changes and the optimal value at these points
opt.u = SS_opt{2}.u(:,51);
opt.x = SS_opt{2}.x(:,51);
opt.z = SS_opt{2}.z(:,51);
opt.t = 0;
tic
for k = 2:length(in_int)
    if any(in_int(:,k)-in_int(:,k-1))
        opt.t = [opt.t; t.samp(k); t.samp(k)];
        for l = 1:4
            if in_int(l,k)-in_int(l,k-1) ~= 0 && in_int(l,k)-in_int(l,1) ~=0
                opt.u = [opt.u opt.u(:,end) SS_opt{l}.u(:,SS_opt{l}.Var == in_int(l,k))];
                opt.x = [opt.x opt.x(:,end) SS_opt{l}.x(:,SS_opt{l}.Var == in_int(l,k))];
                opt.z = [opt.z opt.z(:,end) SS_opt{l}.z(:,SS_opt{l}.Var == in_int(l,k))];
            end
        end
    end
end
toc
opt.t = [opt.t; t.samp(end)];
opt.u = [opt.u opt.u(:,end)];
opt.x = [opt.x opt.x(:,end)];
opt.z = [opt.z opt.z(:,end)];

% Definition of integer values
itg.fig = 1;
itg.lab = 1;
itg.arr = 1;
itg.fda = 1;

% Definition of figure and values
fig.dimw = 12;
fig.dimh = 8;
fig.Font = [16 14];
fig.label = [];

% Plotting of the temperature of the first and the final stage
fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters');
set(gcf, 'position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
yyaxis left
hold on
fig.axes(itg.fig) = gca;
fig.axes(itg.fig).YAxis(1).Color = 'k';
fig.axes(itg.fig).YAxis(2).Color = 'r';
fig.dat(itg.fda) = plot(t.samp,Tin,'k'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(opt.t,opt.z(31,:),'k--'); itg.fda = itg.fda+1;
fig.label(itg.lab) = ylabel('Inlet temperature \itT\rm_0 [^{\circ}C]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = xlabel('Time \itt\rm [min]'); itg.lab = itg.lab+1;
yyaxis right
fig.dat(itg.fda) = plot(t.samp,sol_Dyn.x(30,:),'r'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(opt.t,opt.x(30,:),'r--'); itg.fda = itg.fda+1;
fig.label(itg.lab) = ylabel('Outlet temperature \itT\rm_{30} [^{\circ}C]'); itg.lab = itg.lab+1;
legend BOXOFF 
itg.fig = itg.fig + 1;  


% Plotting of the input movement
fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters');
set(gcf, 'position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
fig.axes(itg.fig) = gca;
fig.dat(itg.fda) = plot(t.samp(1:end),sol_Dyn.u(1,:)*100,'k'); itg.fda = itg.fda+1;
hold on
fig.dat(itg.fda) = plot(t.samp(1:end),sol_Dyn.u(2,:)*100,'r'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(t.samp(1:end),sol_Dyn.u(3,:)*100,'b'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(opt.t,opt.u(1,:)*100,'k--'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(opt.t,opt.u(2,:)*100,'r--'); itg.fda = itg.fda+1;
fig.dat(itg.fda) = plot(opt.t,opt.u(3,:)*100,'b--'); itg.fda = itg.fda+1;
fig.label(itg.lab) = ylabel('Split ratio [%]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = xlabel('Time \itt\rm [min]'); itg.lab = itg.lab+1;
% fig.label(itg.lab) = legend('Input \itu\rm_1','Input \itu\rm_2','Input \itu\rm_3','Location','NE'); itg.lab = itg.lab+1;
legend BOXOFF 
% ylim([0.15 0.40])
itg.fig = itg.fig + 1;


fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters');
set(gcf, 'position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
fig.axes(itg.fig) = gca;
fig.dat(itg.fda) = plot(t.samp,sol_Dyn.z(30,:)*100,'k'); itg.fda = itg.fda+1;
hold on
fig.dat(itg.fda) = plot(opt.t,opt.z(30,:)*100,'k--'); itg.fda = itg.fda+1;
fig.label(itg.lab) = ylabel('Mass fraction NH_3 \itw\rm_{30} [wt.%]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = xlabel('Time \itt\rm [min]'); itg.lab = itg.lab+1;
legend BOXOFF 
itg.fig = itg.fig + 1;
% ylim([28 35])

% set(fig.label, 'Interpreter', 'Latex');
% set(fig.axes,'TickLabelInterpreter', 'Latex');
% set(fig.label(6), 'Fontsize', 12);
set(fig.label, 'Fontsize', fig.Font(1));
set(fig.axes,'Fontsize', fig.Font(2));
set(fig.axes,'LineWidth', 1);
set(fig.dat,'LineWidth', 0.75);

fig.axes(1).XTick = [0 20 40 60 80 100];
