from pyctp import saftvrmie, pcsaft, saftvrqmie
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.cm import get_cmap
from matplotlib.colors import Normalize

comps = 'XE,AR'
eos = saftvrmie.saftvrmie()
eos.init(comps)

temps = np.linspace(115, 125, 5)
norm = Normalize(vmin=min(temps), vmax=max(temps))
cmap = get_cmap('cool')
for T in temps:
    _, LL1V, _= eos.get_binary_pxy(T)

    plt.plot(LL1V[0], LL1V[2]*1e-5, color=cmap(norm(T)), label=int(round(T,0)))
    plt.plot(LL1V[1], LL1V[2]*1e-5, color=cmap(norm(T)), linestyle='--')

c1, c2 = comps.split(',')
plt.title(comps)
plt.xlabel(r'$x_{'+c1+'}$')
plt.ylabel('p [bar]')
plt.legend(title='T [K]')
plt.savefig('Twophase_'+comps)

exit(0)
comps = 'AR,C1'
max_p = 1e8
min_T = 20

s = 1.1
e = 1

eos = saftvrmie.saftvrmie()
eos.init(comps)
#params = eos.get_pure_fluid_param(1)
#eos.set_pure_fluid_param(2, params[0], params[1], s * params[2], e * params[3], params[4])
#eos.redefine_critical_parameters()

T = np.linspace(30,110)
p = np.array([35, 100, 200]) * 1e5
rho_saft = np.array([[1/eos.specific_volume(Ti, pi, [1], 1)[0] for Ti in T] for pi in p])

T_saft, p_saft = eos.get_envelope_twophase(1e5, [1], maximum_pressure=max_p, minimum_temperature=min_T)

eos = saftvrqmie.saftvrqmie()
eos.init(comps)
T_pc_saft, p_pc_saft = eos.get_envelope_twophase(1e5, [1], maximum_pressure=max_p, minimum_temperature=min_T)
rho_qu = np.array([[1/eos.specific_volume(Ti, pi, [1], 1)[0] for Ti in T] for pi in p])

'''

colors = ['r', 'g', 'b']
for i in range(3):
    plt.plot(T, rho_saft[i]/1e3, color=colors[i], label=str(p[i]/1e5) + ' [bar]')
    plt.plot(T, rho_qu[i]/1e3, color=colors[i], linestyle='--')
plt.xlabel('T [k]')
plt.ylabel(r'$\rho$ [kmol/m$^3$]')
plt.legend()
plt.savefig('Ne_supercrit_density')

'''

print('T (saft) [K]: ', T_saft)
print('p (saft) [bar]: ', p_saft/1e5)

print('T (pc-saft) [K]: ', T_pc_saft)
print('p (pc-saft) [bar]: ', p_pc_saft/1e5)

plt.plot(T_saft, p_saft/1e5, label='SAFT-VR-MIE')
plt.plot(T_pc_saft, p_pc_saft/1e5, label='SAFT-VRQ-MIE')
plt.legend()
plt.xlabel('T [K]')
plt.ylabel('p [bar]')
plt.yscale('log')
plt.title(comps.replace(',','-')+', (0.5, 0.5)')
#plt.title('s = '+str(s)+', e = '+str(e))
#plt.title('Neon')
#plt.savefig('Tp_phase_envelope_'+'s_'+str(s).replace('.','_')+'_e_'+str(e))
#plt.savefig('Tp_phase_envelope_Ne')
plt.savefig('Tp_phase_envelope_'+comps.replace(',','-'))
plt.clf()

exit(0)

rho = 1 / np.array([eos.specific_volume(T_saft[i], p_saft[i], [1], 1) for i in range(len(p_saft))])

plt.plot(rho/1e3, T_saft)
plt.xlabel(r'$\rho$ [kmol/m$^3$]')
plt.ylabel('T [K]')
plt.savefig('rho-T_phase_envelope_Ne')