clear ;
addpath(genpath('.'))

% Compression ratio
r = 10 ; % [-]

% Isoentropic efficiencies ;
etaCompression = 0.8 ;
etaExpansion   = 0.8 ;

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

% 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 = 100 * kPa ; % [Pa]
state1.v = getVariableFromStateEquation("v", state1, gas) ;

% 1-2 (isoentropic compression)
state2s.p   = state1.p * r ;
state2s     = getStateAfterIsoentropicProcessConst(state1, state2s, gas) ;
deltaT_isoC = state2s.T - state1.T ;

% non-isoentropic compression
state2.p  = state2s.p ;
deltaT_12 = deltaT_isoC / etaCompression ;
state2.T  = state1.T + deltaT_12 ;
state2.v  = getVariableFromStateEquation("v", state2, gas) ;

% 2-3 (p=const)
state3.T = 1400 ; % [K]
state3 = getStateAfterIsobaricProcess(state2, state3, gas) ;

% 3-4 (isoentropic expansion)
state4s.p   = state1.p ;
state4s     = getStateAfterIsoentropicProcessConst(state3, state4s, gas) ;
deltaT_isoE = state3.T - state4s.T ;

% non-isoentropic expansion
state4.p  = state4s.p ;
deltaT_43 = deltaT_isoE * etaExpansion ;
state4.T  = state3.T - deltaT_43 ;
state4.v  = getVariableFromStateEquation("v", state4, gas) ;

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

% Compression (1-2) & expansion (3-4, turbine) work per unit mass & rlr
lc  = gas.cp * ( state2.T - state1.T ) ;
lt  = gas.cp * ( state3.T - state4.T ) ;
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