#eos = cubic.cubic()
eos = cpa.cpa()
eos.init('PROP1OL,H2O', 'SRK')

x_line = np.linspace(0.0001,0.9999,100)
A1_list = np.zeros(len(x_line))
A2_list = np.zeros(len(x_line))
mu1 = np.zeros(len(x_line))
mu2 = np.zeros(len(x_line))

mu1g = np.zeros(len(x_line))
mu2g = np.zeros(len(x_line))

dmudn = np.zeros(len(x_line))
T = 400
for i, x in enumerate(x_line):
    v = eos.specific_volume(T, 1e5, [x, 1-x], 1)[0]
    vg = eos.specific_volume(T, 1e5, [x, 1-x], 2)[0]

    A1_list[i] = eos.helmholtz_tv(T, v, [x, 1 - x], 1)[0]
    A2_list[i] = eos.helmholtz_tv(T, v, [x, 1 -x], 2)[0]

    mu = eos.chemical_potential_tv(T, v, [x, 1-x], 1)[0]
    mu1[i], mu2[i] = mu

    mug = eos.chemical_potential_tv(T, v, [x, 1 - x], 2)[0]
    mu1g[i], mu2g[i] = mug

dmix_A = (A1_list - (x_line*A1_list[-1] - x_line*A1_list[0]) - A1_list[0])
dmix_mu1 = (mu1 - (x_line*mu1[-1] - x_line*mu1[0]) - mu1[0])
dmix_mu2 = (mu2 - (x_line*mu1[-1] - x_line*mu2[0]) - mu2[0])
dmix_mu = x_line * dmix_mu1 + (1 - x_line) * dmix_mu2

fig, ax = plt.subplots()
twin = ax.twinx()

a_line, = ax.plot(x_line,  dmix_A* 1e-3,
         label=r'$\Delta_{mix}A$')
# twin.plot(x_line, dmix_mu1 * 1e-3,
#           label=r'$\Delta_{mix}\mu1$', color='red')
# twin.plot(x_line,  dmix_mu2 * 1e-3,
#           label=r'$\Delta_{mix}\mu2$', color='green')
# twin.plot(x_line,  dmix_mu * 1e-3,
#           label=r'$\Delta_{mix}\mu$', color='orange')
dadx_line, = twin.plot(x_line[:-1], np.diff(dmix_A)/np.diff(x_line), label=r'$\frac{dA_{mix}}{dx}$', color='red')
#twin.plot(x_line[:-2], np.diff(np.diff(dmix_A)/np.diff(x_line))/np.diff(x_line[:-1]),
#          label=r'$\frac{d^2A_{mix}}{dx^2}$', color='green')

twin.plot()

ax.set_xlim(0,1)
ax.set_ylabel('A')
twin.set_ylabel(r'$\mu$')

plt.legend(handles=[a_line, dadx_line], loc='upper center')
plt.title('Mixing of H2O and 1-propanol')
plt.savefig('Mixing')