r"""
UNIQUAC: separating molecular size from molecular energy
========================================================

Abrams and Prausnitz's UNIQUAC model
(:class:`~chemistrykit.thermo.UNIQUACSolution`) writes
:math:`\ln\gamma_i` as the sum of two physically distinct parts. The
*combinatorial* part depends only on each molecule's van der Waals
volume :math:`r` and surface area :math:`q`, and arises because big and
small molecules pack differently. The *residual* part depends on the
interaction energies. Below, for ethanol-water with tabulated
:math:`r, q` values, the two parts are plotted separately: the size
mismatch alone contributes a small negative (Flory-Huggins-like) term,
while the energetic part drives the large positive deviations.
"""

# %%
import matplotlib.pyplot as plt
import numpy as np

from chemistrykit.thermo import UNIQUACSolution

r = dict(ethanol=2.1055, water=0.92)
q = dict(ethanol=1.972, water=1.40)
full = UNIQUACSolution(r["ethanol"], r["water"], q["ethanol"], q["water"], tau12=0.70, tau21=0.95, P1_star=72.2, P2_star=31.2)
size_only = UNIQUACSolution(r["ethanol"], r["water"], q["ethanol"], q["water"], tau12=1.0, tau21=1.0, P1_star=72.2, P2_star=31.2)

x = np.linspace(0.001, 0.999, 400)
g1, g2 = full.activity_coefficients(x)
c1, c2 = size_only.activity_coefficients(x)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))
for ax, total, comb, name in ((axes[0], g1, c1, "ethanol"), (axes[1], g2, c2, "water")):
    ax.plot(x, np.log(total), color="black", linewidth=2, label=r"total $\ln\gamma$")
    ax.plot(x, np.log(comb), "--", color="steelblue", label="combinatorial (size and shape)")
    ax.plot(x, np.log(total) - np.log(comb), ":", color="darkorange", label="residual (energy)")
    ax.axhline(0.0, color="gray", linewidth=0.8)
    ax.set_xlabel("mole fraction of ethanol")
    ax.set_ylabel(rf"$\ln\gamma_{{{name}}}$")
    ax.set_title(f"UNIQUAC decomposition for {name}")
    ax.legend()
fig.tight_layout()

# %%
# With all :math:`\tau=1` only the combinatorial part is left, and it
# vanishes too when the molecules are the same size:

same = UNIQUACSolution(1.5, 1.5, 1.2, 1.2, 1.0, 1.0, 1.0, 1.0)
print("same-size, athermal gamma:", [round(float(g), 12) for g in same.activity_coefficients(0.3)])
x_az, P_az = full.azeotrope()
print(f"full model azeotrope: x_ethanol = {x_az:.3f} at {P_az:.1f} kPa")

plt.show()