clear ;
addpath(genpath('.'))

% Compression ratio
r = 10 ; % [-]

% Units
kPa = 1e3 ;    % [Pa]
g   = 1/1000 ; % [kg]
kW  = 1000 ;   % [W]
kJ  = 1000 ;   % [J]

% Ideal gas, air
gas.Rbar = 8.314 ;            % [J/(mol*K)]
gas.M    = 28.97 * g ;        % [kg/mol]
gas.R    = gas.Rbar / gas.M ; % [J/(kg*K)]

% States
state1.T = 300 ;       % [K]
state1.p = 100 * kPa ; % [Pa]
state1.v = getVariableFromStateEquation("v", state1, gas) ;
state1.pr = getValueFromTable("pr", "T", state1.T) ;      % [-]
state1.h  = getValueFromTable("h",  "T", state1.T) * kJ ; % [J/kg]

% 1-2 (isoentropic compression)
state2.p  = state1.p  * r ;
state2.pr = state1.pr * r ;
state2.T  = getValueFromTable("T", "pr", state2.pr) ; % [-]
state2.h  = getValueFromTable("h",  "T", state2.T) * kJ ; % [J/kg]
state2.v  = getVariableFromStateEquation("v", state2, gas) ;

% 2-3 (p=const)
state3.T  = 1400 ; % [K]
state3    = getStateAfterIsobaricProcess(state2, state3, gas) ;
state3.pr = getValueFromTable("pr", "T", state3.T) ;      % [-]
state3.h  = getValueFromTable("h",  "T", state3.T) * kJ ; % [J/kg]

% 3-4 (isoentropic expansion)
state4.p  = state3.p  / r ;
state4.pr = state3.pr / r ;
state4.T  = getValueFromTable("T", "pr", state4.pr) ;     % [-]
state4.h  = getValueFromTable("h",  "T", state4.T) * kJ ; % [J/kg]

% Work per unit mass [J/kg]
qIn  = state3.h - state2.h ;
qOut = state4.h - state1.h ;
l = qIn - qOut ;

% Compression (1-2) & expansion (3-4, turbine) work per unit mass & rlr
lc  = state2.h - state1.h ;
lt  = state3.h - state4.h ;
rlr = lc / lt

% Thermal efficiency & restitution ratio
eta = l / qIn 

% Power
volumetricFlowRate1 = 5 ; % [m^3/s]
rho = 1 / state1.v ;
massFlowRate = rho * volumetricFlowRate1 ;
L = massFlowRate * l