% 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',15,...
'DefaultAxesFontName','Times','DefaultAxesFontSize',15,...
'DefaultLineLineWidth',5,'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];

bluecg       =[44,144,195]/256; 
movcg        =[198,102,177]/256; 
greencg      =[83,188,102]/256;
redcg        =[222,49,35]/256;
orangecg     =[255,135,81]/256; 
yellowcg     =[254,229,136]/256; 
greenlightcg =[146,205,91]/256;
greendarkcg  =[21,150,76]/256;

linetype  = {'-','-'};       % line type
col = [bluecg;greencg]; % colors

%% paramerters
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;
zV0      = 0.9;
mA0      = 0;
tsim     = 300;
%% Load data for floating pressure
  
filename = {'OpenLoopmG_Reg';'OpenLoopzV_reg'};
n = length(filename);
xax = [-50 150]; % limit of the xaxis
mG0  = 31.4018; 
figure
for i = 1:n

load(filename{i})
% Inputs: umG - fuel used
Pm   = v_save(:,end-3);                            %[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-4);                   % [bar] pressure at turbine inlet
TGE  = v_save(:,15);                      % [K] flue gas outlet temperature

t    = t_save;                                            % simulation time  


mG   = v_save(:,end-2);                                          % flue gas
zV   = v_save(:,end-1);                              % steam valve position
mA   = v_save(:,end);                                   % economizer bypass
Pm   = Pm/1e3;                                          % [kW] power output
Pm0  = Pm0/1e3;

pSh02 = pSh0;
zV02  = zV0;

figure(1)%power
plot([-50;0;t],[Pm0;Pm0;Pm],linetype{i},'Color',col(i,:))
hold on
xlim(xax)
ylim([16.54 16.7])
xlabel('Time, [s]')
ylabel('CV7 (Power), [MW]')


figure(2) % pressure
plot([-50;0;t],[pSh0;pSh0;pSh],linetype{i},'Color',col(i,:))
hold on
xlim(xax)
ylim([87.2 89])
xlabel('Time, [s]')
ylabel('CV8 (Pressure), [bar]')

figure(3) % fuel
plot([-50;0;t],[mG0;mG0;mG],linetype{i},'Color',col(i,:)) 
hold on
xlim(xax)
ylim([31.3 31.8])
xlabel('Time, [s]')
ylabel('MV1 (Fuel), [kg s$^{-1}$]')

figure(4) % valve
plot([-50;0;t],[zV0;zV0;zV],linetype{i},'Color',col(i,:)) 
hold on
xlim(xax)
ylim([0.89 1.01])
xlabel('Time, [s]')
ylabel('MV2 (Steam Valve), [-]')

end

figname = {'OL_power.eps','OL_pressure.eps','OL_fuel.eps','OL_valve.eps'};
lgd     = {'Step on MV1','Step on MV2'};
for i=1:4    
    figure(i)
    legend(lgd,'Location','east')
    set(gcf,'Position',[100 100 500 250]) 
    legend boxoff
    print('-depsc2','-r600',figname{i}) 
end

