# Density of carbon dioxide at 100 bar and 25C in lbm/in^3
rho = props('D', 'CO2', {T:25 degC, P:100 bar}) in lbm/in^3

# Saturated vapor enthalpy [J/kg] of R134a at 25C
enthalpy = props('H', 'R134a', {T:25 celsius, Q:100%})

# Enthalpy (J per kg dry air) as a function of temperature, pressure, 
#    and relative humidity at STP
enthalpyDry = HAprops('H', {T:298.15 K, P:101325 Pa, R:0.5})

# Create an array with empty objects
cycle = [{},{}];

# Temperature of saturated air at the previous enthalpy
cycle[1].T = HAprops('T', {P:1 atm, H:enthalpyDry, R:1.0})

# Temperature of saturated air in Fahrenheit
cycle[2].T = HAprops('T', {H:enthalpyDry, R:1.0, P:1 atm}) to degF

# Saturation temperature of Water at 1 atm in K
props('T', 'Water', {P:101325 Pa, Q: 0})

# Saturated vapor enthalpy of Water at 1 atm in J/kg
H_V = props('H', 'Water', {P:101325 Pa, Q: 1})

# Saturated liquid enthalpy of Water at 1 atm in J/kg
H_L = props('H', 'Water', {P:101325 Pa, Q: 0})

# Latent heat of vaporization of Water at 1 atm in J/kg
H_V - H_L
# Get the density of Water at T = 461.1 K and P = 5.0e6 Pa, imposing the liquid phase
props('D', 'Water', {'T|liquid':461.1 K, P: 5e6 Pa})

# Get the density of Water at T = 597.9 K and P = 5.0e6 Pa, imposing the gas phase
props('D', 'Water', {T:597.9 K, 'P|gas': 5e6 Pa})

# You can call the props function directly using an empty object:
props('Tcrit', 'Water', {})

# It can be useful to know what the phase of a given state point is
phase('Water', {P:101325 Pa, Q:0})

# The phase index (as floating point number) can also be obtained using the props function:
props('Phase', 'Water', {P:101325 Pa, Q: 0})

# c_p using c_p
props('C', 'Water', {P:101325 Pa, T: 300 K})

# c_p using derivate
props('d(Hmass)/d(T)|P', 'Water', {P:101325 Pa, T: 300 K})

# c_p using second partial derivative
props('d(d(Hmass)/d(T)|P)/d(Hmass)|P', 'Water', {P:101325 Pa, T: 300 K})

# # Vapor compression cycle
# Specifies all states of a typical vapor compression cycle

# ## Inputs
# 
# Fluid and mass flow
fluid = 'R134a'
mDot  = 1 kg/minute

# Component data
evap  = {T: -20 degC, P_drop: 0 Pa, superHeating: 10 K}
cond  = {T:  40 degC, P_drop: 0 Pa, subCooling  : 10 K}
etaS  = 0.75

# ## Process

#- Define an empty array to add the cycle steps
c = [];

#- Make a function for this specific fluid
p(desiredProperty, currentState) = props(desiredProperty, fluid, currentState);

#- Define low and high pressure
pLow  = p('P', {'T|gas': evap.T, Q: 100%});
pHigh = p('P', {'T|liquid': cond.T, Q: 0%  });

# 4 to 1 Evaporation
c[1] = {};
c[1].P = pLow;
c[1].T = evap.T+ evap.superHeating;
c[1].D = p('D', {'T|gas':c[1].T, P:c[1].P});
c[1].H = p('H', {'T|gas':c[1].T, P:c[1].P});
c[1].S = p('S', {'T|gas':c[1].T, P:c[1].P});
c[1]

# 1 to 2 Compression of vapor
c[2] = {};
c[2].P = pHigh;
H_i        = p('H',{P:c[2].P, S:c[1].S});
c[2].H = (H_i-c[1].H)/etaS + c[1].H;
c[2].T = p('T', c[2]);
c[2].D = p('D', c[2]);
c[2].S = p('S', c[2]);
c[2]

# 2 to 3 Condensation
c[3] = {};
c[3].P = c[2].P - cond.P_drop;
c[3].T = cond.T-cond.subCooling;
c[3].D = p('D', {'T|liquid':c[3].T, P:c[3].P});
c[3].H = p('H', {'T|liquid':c[3].T, P:c[3].P});
c[3].S = p('S', {'T|liquid':c[3].T, P:c[3].P});
c[3]

# 3 to 4 Expansion
c[4] = {};
c[4].H = c[3].H;
c[4].P = c[1].P + evap.P_drop;
c[4].T = p('T', c[4]);
c[4].D = p('D', c[4]);
c[4].S = p('S', c[4]);
c[4]

# ## Work, Energy and Performance
# $W_{comp}$
W_comp   = mDot*(c[2].H - c[1].H)

# $Q_h$
Q_h      = mDot*(c[2].H - c[3].H)
# $Q_c$
Q_c      = mDot*(c[1].H - c[4].H)

# $COP_{evaporator}$
evap_COP = Q_c/W_comp

# $COP_{condenser}$
cond_COP = Q_h/W_comp

# # Final results

print('Compressor power   : $1 \t$2\t$3', W_comp to [W, BTU/h, TR], 4)
print('Condenser heat out : $1 \t$2\t$3', Q_h    to [W, BTU/h, TR], 4)
print('Evaporator heat in : $1 \t$2\t$3', Q_c    to [W, BTU/h, TR], 4)

print('COP(cooling)       : $1', [evap_COP], 3)
print('COP(heating)       : $1', [cond_COP], 3)