"""f21_srk_psat -- predicted reduced vapour pressure of propane vs T_r from
RK, SRK and PR, plotted as relative deviation (%) from the SRK prediction
used as the reference curve (SRK sits identically at 0%, RK and PR deviate
from it). Not a comparison to experiment; annotated explicitly.

Uses the shared isofugacity solver in figB_common.py (same algorithm as
figA_16.py). One value is checked against a direct Maxwell equal-area
integration.
"""
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)
Pr_srk, Pr_rk, Pr_pr = [], [], []
for Tr in Trs:
    T = Tr * d["Tc"]
    Pr_srk.append(Psat_solve("srk", sp, T) / d["Pc"])
    Pr_rk.append(Psat_solve("rk", sp, T) / d["Pc"])
    Pr_pr.append(Psat_solve("pr", sp, T) / d["Pc"])
Pr_srk = np.array(Pr_srk); Pr_rk = np.array(Pr_rk); Pr_pr = np.array(Pr_pr)

dev_rk = 100 * (Pr_rk - Pr_srk) / Pr_srk
dev_pr = 100 * (Pr_pr - Pr_srk) / Pr_srk

# numerical check: verify one SRK point by Maxwell equal-area, and confirm
# the SRK reference deviation is exactly zero by construction
Tchk = 0.85 * d["Tc"]
Psat_chk = Psat_solve("srk", sp, Tchk)
maxwell_equal_area_check("srk", sp, Tchk, Psat_chk)
print(f"[check] SRK reference deviation from itself at all Tr: "
      f"max|dev|={np.nanmax(np.abs(100*(Pr_srk-Pr_srk)/Pr_srk)):.2e}% (should be 0 by construction; "
      f"{np.sum(np.isnan(Pr_srk))} of {len(Pr_srk)} points near Tr->1 excluded, isofugacity "
      f"solver loses the bracket a few K below Tc)")

fig, ax = newfig(6.8, 4.6)

ax.axhline(0.0, color=NAVY, lw=2.1, zorder=3)
ax.plot(Trs, dev_rk, color=SKYBLUE, lw=2.0)
ax.plot(Trs, dev_pr, color=AMBER, lw=2.0)

# the SRK reference label is dropped into the empty band opened below the
# zero line (ylim), clear of the descending RK curve that used to cut it
ax.text(0.99, -2.2, "SRK (reference curve,\nnot experiment)", color=NAVY,
        fontsize=10.5, ha="right", va="top", linespacing=1.3)
# RK label parked below its own curve, off the peak; PR label lifted off the
# navy zero line into the gap between the two curves
i_rk = np.argmin(np.abs(Trs - 0.66))
ax.text(0.66, dev_rk[i_rk] - 2.4, "RK", color=SKYBLUE,
        fontsize=10.5, ha="center", va="top")
i_pr = np.argmin(np.abs(Trs - 0.555))
ax.text(0.555, dev_pr[i_pr] + 1.8, "PR", color=AMBER,
        fontsize=10.5, ha="center", va="bottom")

ax.set_xlabel(r"Reduced temperature, $T_r$")
# set over two lines so it can be typeset larger (the stacked sub/superscripts
# were unreadable at the old single-line size) without running off the figure
ax.set_ylabel(r"$100\times\left(P_r^{sat,\,model}-P_r^{sat,\,SRK}\right)$"
              "\n" r"$/\;P_r^{sat,\,SRK}$   [%]",
              fontsize=14.5, labelpad=6, linespacing=1.6)
ax.set_xlim(0.5, 1.0)
ax.set_ylim(-6.8, 37.0)

save(fig, "f21_srk_psat")
