import CoolProp.CoolProp
from CoolProp.CoolProp import PropsSI


# By enabling the following section of the code, you can change the reference state of CoolProp to make it match the one from the psychrometric chart (s = 0 and h = 0 at 273.15K for air; s = 0 and h = 0 for saturated liquid water at 273.15K).
###
Dmolar = CoolProp.CoolProp.PropsSI('Dmolar', 'T', 273.15, 'P', 101325, 'Air')
CoolProp.CoolProp.set_reference_state('Air', 273.15, Dmolar, 0, 0)
Dmolar = CoolProp.CoolProp.PropsSI('Dmolar', 'T', 273.15, 'Q', 0, 'Water')
CoolProp.CoolProp.set_reference_state('Water', 273.15, Dmolar, 0, 0)

# Test for dry air at 0 degrees C:
#h_test1 = 0.001 * PropsSI('H','T', 273.15, 'P', 101325, 'Air') # kJ/kg
#print(f"htest1 = {h_test1:.3f} kJ/kg") # Should be near zero

# Test for saturated air at 0 degrees C:
#h_air = 0.001 * PropsSI('H','T', 267.45, 'P', 101325, 'Air') # kJ/kg
#h_vap = 0.001 * PropsSI('H','T', 267.45, 'Q', 1, 'Water') # kJ/kg

#h_test2 = 0.0023/1.0023 * h_vap + 1/1.0023*h_air
#print(f"htest2 = {h_test2:.3f} kJ/kg") # Should be near zero

###


# Given data
m_w1  = 1.3738      # kg/s
h_w1  = 0.001 * PropsSI('H','T', 323.15, 'P', 101325, 'Water') # kJ/kg  (h_f at 50°C)
print(f"h_w1 = {h_w1:.3f} kJ/kg")
m_a   = 1.9981      # kg/s
omega3 = 0.0073177  # kg_vapor / kg_dry_air
omega4 = 0.027581 # kg_vapor / kg_dry_air
P_fan = 8.0         # kW = kJ/s
m_w2  = 1.333       # kg/s

# Calculate enthalpies of air and water vapour
h_a3 = 0.001 * PropsSI('H','T', 290.15, 'P', 101325, 'Air') # kJ/kg  (h_a at 17°C)
h_a4 = 0.001 * PropsSI('H','T', 303.15, 'P', 101325, 'Air') # kJ/kg  (h_a at 30°C)

h_v3 = 0.001 * PropsSI('H','T', 290.15, 'P', 101325 * omega3 /(0.622+omega3), 'Water') # kJ/kg  (h_v at 17°C)
h_v4 = 0.001 * PropsSI('H','T', 303.15, 'Q', 1, 'Water') # kJ/kg  (h_v at 30°C)

# Energy balance:
# ( m_w1*h_w1 + m_a*(h_a3 + omega3*h_v3 - h_a4 - omega4*h_v4) + P_fan ) / m_w2 = h_f(T2)

numerator = (
    m_w1 * h_w1
    + m_a * (h_a3 + omega3*h_v3 - h_a4 - omega4*h_v4)
    + P_fan
)

h_w2_T2 = numerator / m_w2

T2 = PropsSI('T','H', h_w2_T2*1000, 'P', 101325, 'Water') - 273.15

print(f"h_w2_T2 = {h_w2_T2:.3f} kJ/kg")
print(f"h_a3 = {h_a3:.3f}, ha4 = {h_a4:.3f}, h_v3 = {h_v3:.3f}, h_v4 = {h_v4:.3f} kJ/kg")

print(f"T2 = {T2:.3f} C")