% Script to plot the saved data from a dynamic simulation of a H2P cycle
% 16.10.2019.
% Figure(1):power; figure(2):valve; figure(3):fuel
% setpoint step

clc;clear;close all
%% Set default options for plotting. You can change it to your preferences
set(0,'DefaultTextFontName','Times','DefaultTextFontSize',20,...
'DefaultAxesFontName','Times','DefaultAxesFontSize',20,...
'DefaultLineLineWidth',6,'DefaultLineMarkerSize',7.75,...
'DefaultTextInterpreter','latex')
%% Define new colors

green    = [0       0.5      0];
gold     = [0.9805, 0.7773,  0.25];
hotpink  = [0.8945, 0.4102,  0.5898];
lilac    = [0.7383, 0.5898,  0.7930];
indigo   = [0.293 0 0.5078]; 
cerulean = [0.1094, 0.6641,  0.8359];

R        = 8.314462618*1e-5;            % [m3*bar/K/mol] ideal gas constant
Mwt      = 18*1e-3;                             % [kg/mol] water molar mass
MSh0     = 100;                                   % [kg] superheater holdup
pSh0     = 8.805290778575733e+01;            % [bar] superheater pressure
pShsp    = pSh0;
pT0      = pSh0-5;                % [bar] nominal pressure at turbine inlet
TSh0     = 868;
VSh      = MSh0/Mwt*R*TSh0/pSh0;  
TGE0     = 150+273;                  % [K] flue gas outlet temperature. nominal
zV0      = 0.9;
tsim     = 300;

%% Load data for floating pressure
%filename = 'FloatingPressure_zV1' ;    
filename = 'FloatingPressure' ;    

load(filename)
% Inputs: umG - fuel used
Pm   = v_save(:,end-2);                            %[kW] - mechanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure
pT   = v_save(:,end-3);                   % [bar] pressure at turbine inlet
TGE  = v_save(:,15);                      % [K] flue gas outlet temperature

t    = t_save;                                            % simulation time  

mG0  = ctrlPm(1); 
mG   = v_save(:,end-1);                                          % flue gas
zV   = v_save(:,end);                                % steam valve position

Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output

pSh02 = pSh0;87.5522;
zV02  = zV0;1;
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',indigo)
hold on

figure(2)%pressure
plot([-50;0;t],[pSh02;pSh02;pSh],'Color',indigo)
hold on

figure(3)% valve 
plot([-50;0;t],[zV02;zV02;zV],'Color',indigo) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',indigo) 
hold on

figure(5) % pressure at turbine inlet
plot([-50;0;t],[pT0;pT0;pT],'Color',indigo)
hold on

figure(6) % flue gas outlet temperature
plot([-50;0;t],[TGE0;TGE0; TGE],'Color',indigo)
hold on
%% Load data for boiler driven
clear t mG zV Pm pSh KcPm tauIPm errPm interPm Pmsp pT
filename = 'BoilerDriven' ;    
load(filename)
% Inputs: mG - fuel used; zV valve position
Pm   = v_save(:,end-2);                             %[kW] - mchanical power
MSh  = v_save(:,5);                                  % superheater holdup
TSh  = v_save(:,3);                             % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                % [bar] superheater pressure
pT   = v_save(:,end-3);                 % [bar] pressure at turbine inlet
TGE  = v_save(:,17);                      % [K] flue gas outlet temperature

mG   = v_save(:,end-1);
zV   = v_save(:,end);

t    = t_save;                                            % simulation time  
 
Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',green)
hold on
xlabel('Time, [s]')
ylabel('Power, [kW]')
xlim([-50,max(t)])

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',green)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zV],'Color',green) 
hold on


figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',green) 
hold on

figure(5) % pressure at turbine inlet
plot([-50;0;t],[pT0;pT0;pT],'Color',green)

figure(6) % flue gas outlet temperature
plot([-50;0;t],[TGE0;TGE0; TGE],'Color',green)

%% Load data for turbine driven
clear t umG uzV Pm pSh
filename = 'TurbineDriven' ;    
load(filename)
% Inputs: mG - fuel used; zV - steam valve
Pm   = v_save(:,end-2);                             %[kW] - mchanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure
pT   = v_save(:,end-3);                   % [bar] pressure at turbine inlet
TGE  = v_save(:,17);                      % [K] flue gas outlet temperature

mG   = v_save(:,end-1);
zV   = v_save(:,end);

t = t_save;


Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',hotpink)
hold on

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',hotpink)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zV],'Color',hotpink) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',hotpink) 
hold on

figure(5) % pressure at turbine inlet
plot([-50;0;t],[pT0;pT0;pT],'Color',hotpink)

figure(6) % flue gas outlet temperature
plot([-50;0;t],[TGE0;TGE0; TGE],'Color',hotpink)
%% Load data for valve position controller
clear t umG uzV Pm pSh pT
filename = 'VPC' ;    
load(filename)
% Inputs: mG - fuel used; zV-steam valve
Pm   = v_save(:,end-2);                             %[kW] - mchanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure
pT   = v_save(:,end-3);                   % [bar] pressure at turbine inlet
TGE  = v_save(:,17);                      % [K] flue gas outlet temperature

mG   = v_save(:,end-1);
zV   = v_save(:,end);

t    = t_save;                                            % simulation time  

Pm0 = Pm0/1e3;                                  % [MW] nominal power output
Pm  = Pm/1e3;                                           % [MW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'--','Color',gold)
hold on,

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'--','Color',gold)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zV],'--','Color',gold) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'--','Color',gold) 
hold on

figure(5) % pressure at turbine inlet
plot([-50;0;t],[pT0;pT0;pT],'--','Color',gold)

figure(6) % flue gas outlet temperature
plot([-50;0;t],[TGE0;TGE0; TGE],'--','Color',gold)

%% Load data for 2 controllers
clear t umG uzV Pm pSh pT
filename = '2Controllers' ;    
load(filename)
% Inputs: mG - fuel used; zV - steam valve position
Pm   = v_save(:,end-2);                             %[kW] - mchanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure
pT   = v_save(:,end-3);                   % [bar] pressure at turbine inlet
mG   = v_save(:,end-1);
zV   = v_save(:,end);
TGE  = v_save(:,15);                      % [K] flue gas outlet temperature 

t    = t_save;                                            % simulation time  

Pm0  = Pm0/1e3;                                 % [MW] nominal power output
Pm   = Pm/1e3;                                          % [MW] power output                                  
Pmsp = Pmsp/1e3;                                      % [MW] power setpoint
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'--','Color',cerulean)
plot([-50;0;t],[Pm0;Pm0;ones(length(t),1)*Pmsp],'k:') % plot setpoint
xlabel('Time, [s]')
ylabel('CV7 (Power), [MW]')
xlim([-50,300])
ylim([16 17.5])
textlegend = {'Floating Pressure','Boiler Driven','Turbine Driven',...
       'VPC','2 Controllers','Setpoint'};
legend(textlegend,'Location','northeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Power_dist.eps')                       % save eps 

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'--','Color',cerulean)
plot([-50;0;t],[pSh0;pSh0;ones(length(t),1)*pSh0],'k:') % plot setpoint
xlim([-50,300])
ylim([86 93])
xlabel('Time, [s]')
ylabel('CV8 (Pressure), [bar]')
textlegend = {'Floating Pressure','Boiler Driven','Turbine Driven',...
       'VPC','2 Controllers','Setpoint'};
legend(textlegend,'Location','northeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Pressure_dist.eps')                      % save eps 


figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zV],'--','Color',cerulean) 
hold on
xlabel('Time, [s]')
ylabel('MV2 (Steam valve), [-]')
xlim([-50,300])
ylim([0 1.05])
textu = {'Floating Pressure','Boiler Driven','Turbine Driven','VPC','2 Controllers'};
legend(textu,'Location','southeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Valve_dist.eps')                       % save eps 


figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'--','Color',cerulean) 
hold on
xlabel('Time, [s]')
ylabel('MV1 (Fuel), [kg s$^{-1}$]')
xlim([-50,300])
ylim([28; 32])
legend(textu,'Location','northeast','NumColumns',1)
set(gcf,'Position',[100 100 750 500]) 
legend boxoff% set the figure size                 
print('-depsc2','-r600','Fuel_dist.eps')                        % save eps 

figure(5) % pressure at turbine inlet
plot([-50;0;t],[pT0;pT0;pT],'--','Color',cerulean)

figure(6) % flue gas outlet temperature
plot([-50;0;t],[TGE0;TGE0; TGE],'--','Color',cerulean)
hold on
plot([-50;0;t],[TGE0;TGE0;ones(length(t),1)*TGE0],'k:') % plot setpoint
xlabel('Time, [s]')
ylabel('Flue gas outlet temperature, [K]')
xlim([-50,300])
%ylim([30.8 31.6])
textlegend = {'Floating Pressure','Boiler Driven','Turbine Driven',...
       'VPC','2 Controllers','T_E^{g,min}'};
legend(textlegend,'Location','southeast','NumColumns',1)
set(gcf,'Position',[100 100 750 500]) 
legend boxoff% set the figure size                 
print('-depsc2','-r600','TGE_dist.eps')   