r"""
Andrews and van der Waals: continuity of the liquid and gas states
==================================================================

Andrews's CO2 isotherms showed a flat liquid-vapor plateau below a
critical temperature and none above it. The van der Waals equation
(:class:`~chemistrykit.thermo.systems.equations_of_state.VanDerWaals`),
built from CO2's critical constants with
:meth:`~chemistrykit.thermo.systems.equations_of_state.VanDerWaals.from_critical_constants`,
reproduces the whole family: loops with three volume roots below
:math:`T_c`, an inflection at :math:`T_c`, and single-valued isotherms
above it, so gas can be turned into liquid continuously by going around
the critical point.
"""

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

from chemistrykit.thermo.systems.equations_of_state import VanDerWaals

Tc, Pc = 304.13, 7.3773e6  # CO2
vdw = VanDerWaals.from_critical_constants(Tc, Pc)
Vc = 3.0 * vdw.b

Vm = np.linspace(1.4 * vdw.b, 8.0 * Vc, 600)
fig, ax = plt.subplots(figsize=(7, 5))
for Tr in [0.85, 0.9, 0.95, 1.0, 1.1, 1.2]:
    style = "k-" if Tr == 1.0 else "-"
    ax.plot(Vm / Vc, vdw.pressure(Vm, Tr * Tc) / Pc, style, label=rf"$T/T_c$ = {Tr}")
ax.plot([1.0], [1.0], "ro", label="critical point")
ax.set_xlim(0.4, 8.0)
ax.set_ylim(0.0, 2.0)
ax.set_xlabel(r"$V_m / V_c$")
ax.set_ylabel(r"$P / P_c$")
ax.set_title("van der Waals isotherms of CO2 around the critical point")
ax.legend()
fig.tight_layout()

# %%
# Below :math:`T_c` the cubic has three positive roots at a given
# pressure (liquid, unstable, vapor); above it only one. At the critical
# point the universal compressibility factor is exactly 3/8.

for Tr, Pr in [(0.9, 0.6), (1.2, 0.6)]:
    coeffs = [1.0, -(vdw.b + vdw.R * Tr * Tc / (Pr * Pc)), vdw.a / (Pr * Pc), -vdw.a * vdw.b / (Pr * Pc)]
    roots = np.sort(np.roots(coeffs).real[np.abs(np.roots(coeffs).imag) < 1e-12])
    print(f"T/Tc = {Tr}, P/Pc = {Pr}: {roots.size} real root(s), V/Vc = {np.round(roots / Vc, 3)}")
print(f"Zc = Pc Vc / (R Tc) = {Pc * Vc / (vdw.R * Tc):.4f}")

plt.show()