"""f22_pr_liquid_density -- saturated liquid molar volume of propane vs T_r
from SRK and PR (isofugacity liquid root V_liq = Z_liq RT/P at each
model's own Psat, SPEC.md section 8 constants). PR's smaller Zc (0.3074 vs
0.3333) makes its liquid-root branch systematically smaller than SRK's,
qualitatively mimicking the effect that a posteriori volume translation
(Peneloux, Rauzy & Freze, 1982, SPEC.md reference list) is designed to
correct; this figure shows the built-in PR effect, not an applied
translation (no translation constant is invented here).
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from figB_common import *

sp = "propane"
d = SPECIES[sp]

Trs = np.linspace(0.5, 0.99, 60)
Vl_srk, Vl_pr = [], []
Tr_used = []
for Tr in Trs:
    T = Tr * d["Tc"]
    Ps = Psat_solve("srk", sp, T)
    Pp = Psat_solve("pr", sp, T)
    if np.isnan(Ps) or np.isnan(Pp):
        continue
    rs = Zroots("srk", sp, T, Ps)
    rp = Zroots("pr", sp, T, Pp)
    Vl_srk.append(rs[0] * Rb * T / Ps)
    Vl_pr.append(rp[0] * Rb * T / Pp)
    Tr_used.append(Tr)
Tr_used = np.array(Tr_used); Vl_srk = np.array(Vl_srk); Vl_pr = np.array(Vl_pr)

print(f"[info] at Tr=0.7: V_liq SRK = {np.interp(0.7, Tr_used, Vl_srk):.2f} cm3/mol, "
      f"V_liq PR = {np.interp(0.7, Tr_used, Vl_pr):.2f} cm3/mol "
      f"(PR/SRK ratio = {np.interp(0.7, Tr_used, Vl_pr)/np.interp(0.7, Tr_used, Vl_srk):.3f})")

fig, ax = newfig(6.8, 4.6)

ax.plot(Tr_used, Vl_srk, color=NAVY, lw=2.1)
ax.plot(Tr_used, Vl_pr, color=AMBER, lw=2.1)
ax.fill_between(Tr_used, Vl_pr, Vl_srk, color=SAND, alpha=0.5, zorder=1)

ax.text(0.56, np.interp(0.56, Tr_used, Vl_srk) + 4, "SRK", color=NAVY, fontsize=10.5, va="bottom")
ax.text(0.56, np.interp(0.56, Tr_used, Vl_pr) - 4, "PR", color=AMBER, fontsize=10.5, va="top")

ax.annotate("PR gives the smaller, more realistic\nliquid molar volume\n"
            "(built-in effect of $Z_c^{PR}=0.307 < Z_c^{SRK}=0.333$,\n"
            "qualitatively like a volume translation)",
            xy=(0.80, (np.interp(0.80, Tr_used, Vl_srk)+np.interp(0.80, Tr_used, Vl_pr))/2),
            xytext=(0.60, 145), fontsize=10.5, color=GRAPHITE, linespacing=1.3,
            arrowprops=dict(arrowstyle="-", color=GRAPHITE, lw=0.9))

ax.set_xlabel(r"Reduced temperature, $T_r$")
# two lines: as a single line the label is taller than the axes and the
# closing unit bracket was cut off at the top of the figure
ax.set_ylabel("Saturated liquid molar volume,\n"
              r"$V_{liq}^{sat}$  [cm$^3$ mol$^{-1}$]", linespacing=1.6)
ax.set_xlim(0.5, 1.0)
ax.set_ylim(0, 210)

save(fig, "f22_pr_liquid_density")
