%% Plot effect of changing air flow
close all
%create lines for limits
rpm_max= ones(length(a.t),1)*N_cmax;
rpm_min= ones(length(a.t),1)*N_cmin;

z_max= ones(length(a.t),1);
z_min= zeros(length(a.t),1);

xMeOH_max_line= ones(length(a.t),1)*x_MeOHmax;
xMeOH_min_line= ones(length(a.t),1)*x_MeOHmin;

vmaxPI_line  = ones(length(a.t),1)*vmaxPI;
vminPI_line  =  ones(length(a.t),1)*vminPI;

vmaxPI_ax =vmaxPI+0.2;
vminPI_ax =vminPI-0.2;

%% Main figure
figure
subplot(3,2,1)
plot(a.t, a.x_MeOH)
hold on
plot(a.t, a.x_MeOHsp,'k--')
plot(a.t, xMeOH_max_line,'r--')
plot(a.t, xMeOH_min_line,'r--')
title('MeOH mol fraction')
ylabel('Mol fraction (-)')
%legend('x_{MeOH}','x_{MeOH}^{sp}',
legend('$x_{MeOH}$','$x_{MeOH}^{sp}$','Interpreter','latex')
ylim([x_MeOHmin_ax x_MeOHmax_ax])
ytickformat('%.2f')
% ylim([0.09 .11])
legend boxoff 

subplot(3,2,2)
plot(a.t, a.m_tot)
hold on
plot(a.t, a.m_tot_sp,'k--')
title('Total mass flow')
ylabel('Mass flow (kg/h)')
legend('$\dot{m}_{tot}$','$\dot{m}_{tot}^{sp}$', 'Location', 'Southeast','Interpreter','latex')
%ylim([mtot_min mtot_max])
ylim([mtot_min_ax mtot_max_ax])
ytickformat('%.1f')
legend boxoff 

subplot(3,2,3)
plot(a.t,a.m_air)
hold on
plot(a.t,a.m_air_sp,'k--')
title('Air mass flow')
ylabel('Mass flow (kg/h)')
legend('$\dot{m}_{air}$','$\dot{m}_{air}^{sp}$', 'Location', 'Southeast','Interpreter','latex')
ylim([mair_min_ax mair_max_ax])
%ylim([23000 25000])
legend boxoff 

subplot(3,2,4)
plot(a.t,a.m_MeOH)
hold on
plot(a.t,a.m_MEOH_sp,'k--')
title('MeOH mass flow')
ylabel('Mass flow (kg/h)')
%legend('m_{MeOH}','m_{MeOH}^{sp}', 'Location', 'Southeast')
legend('$\dot{m}_{MeOH}$','$\dot{m}_{MeOH}^{sp}$', 'Location', 'Southeast','Interpreter','latex')
%ylim([2800 3100])
ylim([m_MeOHmin_ax m_MeOHmax_ax])
xlim([0 tend])
legend boxoff 
%This changes that the y axis is in exponential format
ax = gca;
ax.YRuler.Exponent = 3;
ax.YRuler.TickLabelFormat = '%.1f';


subplot(3,2,5)
plot(a.t,a.comp_rpm)
hold on
plot(a.t,rpm_max,'r--')
plot(a.t,rpm_min,'r--')
title('Air compressor')
xlabel('Time (s)')
ylabel('Rotational speed (rpm)') %Rotational speed, \omega_c
ylim([N_cmin-100 N_cmax+100])
ax = gca;
ax.YRuler.Exponent = 3;
ax.YRuler.TickLabelFormat = '%.1f';


subplot(3,2,6)
plot(a.t,a.z_MeOH)
hold on
plot(a.t,z_max,'r--')
plot(a.t,z_min,'r--')
title('MeOH valve')
ylim([-0.05 1.05])
xlabel('Time (s)')
ylabel('Opening (-)')
ytickformat('%.2f')
yticks(0:0.25:1)


%%
figure

subplot(2,2,1)
plot(a.t, a.u_air_scaled)
hold on
plot(a.t, a.u_MeOH_scaled)
plot(a.t,a.v, 'k--')
legend('$u_{air}$', '$u_{MeOH}$', 'v','Interpreter','latex', 'Location', 'Southeast')
ylim([vminPI_ax vmaxPI_ax])
ylabel('Signal')
legend boxoff 
ytickformat('%.2f')
title('Signals')


subplot(2,2,2)
plot(a.t, a.x_MeOH_error)
title ('x_{MeOH} error')
ylabel('error (mol fracc)')
ytickformat('%.3f')


subplot(2,2,3)
plot(a.t, a.m_air_sp)
title('Air mass flow set point')
ylim([mair_min_ax mair_max_ax])
%ylim([23000 25000])
%legend boxoff 
xlabel('Time (s)')
ylabel('Mass flow (kg/h)')

subplot(2,2,4)
plot(a.t, a.m_MEOH_sp_SR)
hold on
plot(a.t, a.m_MeOH_sp_m_tot)
plot(a.t, a.m_MEOH_sp,'k--')
legend('SR', '$\dot{m}_{tot}$', '$\dot{m}_{MeOH}^{sp}$','Interpreter','latex', 'Location','Southeast')
ylim([m_MeOHmin_ax m_MeOHmax_ax])
xlim([0 tend])
title('MeOH set point')
xlabel('Time (s)')
ylabel('Mass flow (kg/h)')
ax = gca;
ax.YRuler.Exponent = 3;
ax.YRuler.TickLabelFormat = '%.1f';
legend boxoff 

%%
% figure
%  plot(a.t, a.Nc_act)
% 
% plot(a.t,a.I_air)
% title('controller output signal')
% legend('I_{0}')

% plot(a.t,a.z_actuator)
% hold on
% plot(a.t,a.wcomp)
% title('integrator outputs')
% legend('z_{actuator}','w')

% figure
% plot(a.t,a.Nc_act)
% % plot(a.t, a.I_air0)
% hold on
% plot(a.t, a.I_air)
% legend('from PI', '+ bias')


