clear ;
addpath(genpath('.'))

% Compression ratio & critic ratio
r =  18 ; % [-]
rc =  2 ; % [-]

% Units
MPa = 1e6 ;    % [Pa]
cm  = 1/100 ;  % [m]
g   = 1/1000 ; % [kg]

% 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)]

% Constant cp and cv
gas.k  = 1.4 ; % [-]
gas.cv = gas.R / ( gas.k - 1 ) ;
gas.cp = gas.k * gas.cv ;

% States
state1.T = 300 ;         % [K]
state1.p = 0.1 * MPa ;   % [Pa]
state1.v = getVariableFromStateEquation("v", state1, gas) ;

% 1-2 (isoentropic compression)
state2.v = state1.v / r ;
state2   = getStateAfterIsoentropicProcessConst(state1, state2, gas) ;

% 2-3 (p=const)
state3.v = state2.v * rc ;
state3 = getStateAfterIsobaricProcess(state2, state3, gas) ;

% 3-4 (isoentropic expansion)
state4.v = state1.v ;
state4   = getStateAfterIsoentropicProcessConst(state3, state4, gas) ;

% Work per unit mass per cycle [J/kg]
qIn  = gas.cp * ( state3.T - state2.T ) ;
qOut = gas.cv * ( state4.T - state1.T ) ;
l = qIn - qOut ;

% Thermal efficiency
eta = l / qIn 

% Mean effective pressure [Pa]
Dv = state1.v - state2.v ;
pme = l / Dv