import CoolProp.CoolProp
from CoolProp.CoolProp import PropsSI
from scipy.optimize import brentq
m_H2 = 1; # kg
m_N2 = 4.667 # kg
T1_H2 = 300
T1_N2 = 600
u1_H2 = PropsSI('U','T', T1_H2, 'P', 600000, 'H2') # kJ/kg
u1_N2 = PropsSI('U','T', T1_N2, 'P', 600000, 'N2') # kJ/kg

P2 = 600000.0  # Pa  (pressure at state 2), assumed constant, as the change is small, as seen in part b)

# -------------------------------------------------------------
# DEFINE THE ENERGY BALANCE RESIDUAL
# -------------------------------------------------------------
def energy_balance(T2):

    # Residual of u1_H2*m1_H2 + u1_N2*m1_N2 - [u2_H2(T2)*m1_H2 + u2_N2(T2)*m1_N2] = 0

    u2_H2 = PropsSI('U', 'T', T2, 'P', P2, 'H2')
    u2_N2 = PropsSI('U', 'T', T2, 'P', P2, 'N2')

    lhs = u1_H2 * m_H2 + u1_N2 * m_N2
    rhs = u2_H2 * m_H2 + u2_N2 * m_N2

    return lhs - rhs

# SOLVE FOR T2

T_low = 300.0   # K
T_high = 400.0 # K

T2_solution = brentq(energy_balance, T_low, T_high)

print(f"T2 = {T2_solution:.2f} K")