clear ;
addpath(genpath('.'))

% Compression ratio
r = 8 ; % [-]

% Units
bar = 1e5 ;    % [Pa]
cm  = 1/100 ;  % [m]
g   = 1/1000 ; % [kg]
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)]

% State 1
state1.T  = 300 ;         % [K]
state1.p  = 1.013 * bar ; % [Pa]
state1.V  = 566 * cm^3 ;  % [m^3]
state1.vr = getValueFromTable("vr", "T", state1.T) ;         % [-]
state1.u  = getValueFromTable("u",  "T", state1.T) * kJ ;    % [J/kg]
gas.m     = getVariableFromStateEquation("m", state1, gas) ; % [kg]

% 1-2 (isoentropic compression)
state2.V  = state1.V  / r ;
state2.vr = state1.vr / r ;
state2.T  = getValueFromTable("T", "vr", state2.vr) ; % [K]
state2.u  = getValueFromTable("u", "T",  state2.T) * kJ ;
state2.p  = getVariableFromStateEquation("p", state2, gas) ;

% 2-3 (V=const)
state3.T  = 2000 ; % [K]
state3    = getStateAfterIsochoricProcess(state2, state3, gas) ;
state3.vr = getValueFromTable("vr", "T", state3.T) ;
state3.u  = getValueFromTable("u",  "T", state3.T) * kJ ;

% 3-4 (isoentropic expansion)
state4.V  = state3.V  * r ;
state4.vr = state3.vr * r ;
state4.T  = getValueFromTable("T", "vr", state4.vr) ;
state4.u  = getValueFromTable("u", "T",  state4.T) * kJ ;

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

% Thermal efficiency
eta = l / qIn 

% Total work per cycle [J]
L = gas.m * l 

% Displacement volume [m^3]  and mean effective pressure [Pa]
V = state1.V - state2.V ;
pme = L / V