clear;close all;clc

% Gas properties
gases = struct(...
    'name',     {'H2',   'CO2', 'CH4'}, ...
    'C_p',      {14304,  839,   2191 }, ... % (J/kg/K)
    'hc_ratio', {1.41,   1.3,   1.31 }, ... % heat capacity ratio
    'p1',       {1,      35,    20   }, ... % (bar)
    'p2',       {100,    20,    5    }, ... % (bar)
    'T1',       {288,    288,   288  }, ... % (K)
    'eta_c',    {0.7,    0.7,   0.7  }  ... % mechanical efficiency
);

nGas = numel(gases);        % number of gases
p2 = linspace(1,100,100);   % Pressure range


% Compute compression work
Wadiab = zeros(nGas, numel(p2));
Wisot  = zeros(nGas, numel(p2));
for i = 1:nGas
    [Wadiab(i,:), Wisot(i,:)] = compr_Q_T(gases(i),p2);
end

[Wadiab_h2, Wisot_h2] = compr_Q_T(gases(2));
x=0.32;
W_tot = Wadiab_h2*x + Wisot_h2*(1-x);

% Plot
figure
hold on
for i = 1:nGas
    h = plot(p2, Wadiab(i,:), 'DisplayName', [gases(i).name '(adiabatic)'], 'LineWidth', 3);
    plot(p2, Wisot(i,:),'--','Color',h.Color,'DisplayName', [gases(i).name '(isothermal)'], 'LineWidth', 3);
end
legend('Location', 'northwest')
grid on
xlabel('Pressure (bar)');
ylabel('Compression cost (kWh/kg)');
ax = gca;
ax.FontSize = ax.FontSize * 2;

% Function
function [W_adiab, W_isot] = compr_Q_T(comp,p2)
    if nargin < 2
        p2 = comp.p2;
    end
    W_adiab = comp.C_p / 3.6e6 * comp.T1 / comp.eta_c ...
              * ((p2 ./ comp.p1) .^ ((comp.hc_ratio-1) / comp.hc_ratio) - 1);
    W_isot  = (comp.hc_ratio-1) / comp.hc_ratio * comp.C_p * comp.T1 / comp.eta_c ...
              / 3.6e6 * log(p2 / comp.p1);
end