"""f23_pvt_envelope_compare -- density-temperature coexistence envelope
(rho_liq and rho_vap vs T) for propane from SRK and PR on the same axes,
each model's own critical point marked (T_c fixed at the SPEC.md value;
critical molar volume V_c = Z_c R T_c / P_c uses each model's own Z_c:
SRK Z_c = 1/3, PR Z_c = 0.3074, SPEC.md section 8).
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from figB_common import *

sp = "propane"
d = SPECIES[sp]
MW = 44.10  # g/mol, propane molar mass (standard value, used only to convert to kg/m3)

Trs = np.linspace(0.5, 0.999, 90)

def envelope(model):
    T_, rl, rv = [], [], []
    for Tr in Trs:
        T = Tr * d["Tc"]
        P = Psat_solve(model, sp, T)
        if np.isnan(P):
            continue
        r = Zroots(model, sp, T, P)
        if len(r) < 2:
            continue
        Vl = r[0] * Rb * T / P
        Vv = r[-1] * Rb * T / P
        T_.append(T)
        rl.append(1000 * MW / Vl)   # kg/m3 (V in cm3/mol)
        rv.append(1000 * MW / Vv)
    return np.array(T_), np.array(rl), np.array(rv)

fig, ax = newfig(7.4, 5.0)

specs = {"srk": (NAVY, "SRK", 1/3), "pr": (AMBER, "PR", 0.3074)}
for model, (col, lbl, Zc) in specs.items():
    T_, rl, rv = envelope(model)
    ax.plot(T_, rl, color=col, lw=2.1)
    ax.plot(T_, rv, color=col, lw=2.1)
    Vc = Zc * Rb * d["Tc"] / d["Pc"]
    rc = 1000 * MW / Vc
    ax.plot([d["Tc"]], [rc], marker="o", ms=7, color=col, zorder=5,
             markeredgecolor="white", markeredgewidth=1.0)
    # both critical-point labels go into the right-hand gutter (see set_xlim),
    # stacked well apart, each joined to its own marker by a short leader; on
    # top of the markers they collided with each other and with the grey
    # vapour-branch label
    y_lab = 250 if model == "pr" else 130
    ax.plot([d["Tc"] + 1.5, d["Tc"] + 6.5], [rc, y_lab], color=col, lw=0.9,
            alpha=0.8, zorder=4)
    ax.text(d["Tc"] + 8.5, y_lab, f"{lbl} critical point", color=col,
            fontsize=10.5, ha="left", va="center")

# direct labels on the liquid branch, placed inside the axes near Tr=0.62
Tlbl = 0.62 * d["Tc"]
for model, (col, lbl, Zc) in specs.items():
    T_, rl, rv = envelope(model)
    idx = np.argmin(np.abs(T_ - Tlbl))
    ax.text(T_[idx], rl[idx] + 18, lbl + " liquid", color=col, fontsize=10.5, ha="center", va="bottom")

# vapour-branch label (curves nearly coincide there; single shared label)
Tlbl_v = 0.90 * d["Tc"]
idxv = np.argmin(np.abs(envelope("pr")[0] - Tlbl_v))
Tv, rlv, rvv = envelope("pr")
ax.annotate("SRK / PR vapour branch\n(nearly coincide)", xy=(Tv[idxv], rvv[idxv]),
            xytext=(Tv[idxv] - 88, 152), fontsize=10.5, color=GRAPHITE, linespacing=1.3,
            ha="left", va="center",
            arrowprops=dict(arrowstyle="-", color=GRAPHITE, lw=0.9))

ax.set_xlabel(r"Temperature, $T$  [K]")
ax.set_ylabel(r"Molar density $\times M$, $\rho$  [kg m$^{-3}$]")
ax.text(0.01, 1.01, "propane, $M=44.10$ g mol$^{-1}$", transform=ax.transAxes,
        fontsize=10.5, color=SLATE, ha="left", va="bottom")
ax.set_xlim(150, 465)
ax.set_xticks(np.arange(150, 376, 50))   # ticks pinned to the data range
ax.set_ylim(0, 720)
fig.subplots_adjust(left=0.16)

save(fig, "f23_pvt_envelope_compare")
