"""f17_acentric_definition -- same axes as f16 (log10 Pr^sat vs 1/Tr, PR EOS,
methane/propane/n-butane/water), with a vertical line at 1/Tr = 1/0.7, the
simple-fluid reference line log10 Pr^sat = -1.0, and amber arrows showing
omega = -1 - log10(Pr^sat)|_{Tr=0.7} for two fluids (propane, water).

Reuses the isofugacity saturation-pressure routine of figA_16.py (same
algorithm, copied into figB_common.py as Psat_solve/ln_phi so that this
script does not import/execute figA_16.py's module-level plotting code).
Axis limits, species set and colours are matched to figA_16.py exactly so
f16 and f17 read as one pair.
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from figB_common import *

species_list = ["methane", "propane", "n-butane", "water"]
colors = {"methane": NAVY, "propane": TEAL, "n-butane": SKYBLUE, "water": RUST}

fig, ax = newfig(6.8, 4.6)

ends = []
curves = {}
for sp in species_list:
    d = SPECIES[sp]
    Trs = np.linspace(0.5, 0.99, 60)
    invTr, logPr = [], []
    for Tr in Trs:
        T = Tr * d["Tc"]
        Ps = Psat_solve("pr", sp, T)
        if np.isnan(Ps):
            continue
        invTr.append(1.0 / Tr)
        logPr.append(np.log10(Ps / d["Pc"]))
    invTr = np.array(invTr); logPr = np.array(logPr)
    curves[sp] = (invTr, logPr)
    ax.plot(invTr, logPr, color=colors[sp], lw=2.1)
    ends.append((sp, invTr[0], logPr[0]))

# direct labels, staggered as in f16
ends_sorted = sorted(ends, key=lambda t: -t[2])
min_gap = 0.30
ytexts = []
for i, (sp, x, y) in enumerate(ends_sorted):
    yt = y if i == 0 else min(y, ytexts[-1] - min_gap)
    ytexts.append(yt)
for (sp, x, y), yt in zip(ends_sorted, ytexts):
    col = colors[sp]
    if abs(yt - y) > 1e-6:
        ax.plot([x, x + 0.03], [y, yt], color=col, lw=0.8, alpha=0.7)
    ax.text(x + 0.045, yt, sp, color=col, fontsize=10.5, ha="left", va="center")

ax.plot([1.0], [0.0], marker="o", ms=6, color=GRAPHITE, zorder=5,
         markeredgecolor="white", markeredgewidth=0.8)

# --- Pitzer definition construction --------------------------------------
invTr07 = 1.0 / 0.7
ax.axvline(invTr07, color=GRAPHITE, lw=1.2, linestyle=":", zorder=2)
ax.text(invTr07, 0.52, r"$T_r=0.7$", color=GRAPHITE, fontsize=10.5,
        ha="center", va="bottom")

ax.axhline(-1.0, color=NAVY, lw=1.3, linestyle="--", zorder=2)
# label parked in the right-hand gutter (see set_xlim below), just above the
# right end of the dashed line so the line does not strike through the text
ax.text(2.12, -0.94, "simple-fluid reference,\n" r"$\log_{10}P_r^{sat}=-1$",
        color=NAVY, fontsize=10.5, ha="left", va="bottom", linespacing=1.3)

# omega arrows for two fluids computed directly from the PR isofugacity
# solution at Tr = 0.7 (consistent with the curves already plotted, not
# re-fitted): omega = -1 - log10(Pr^sat(Tr=0.7)). Text placed well clear of
# the two amber points and of each other (propane above, water below).
# propane: open box above the dashed reference line and right of the Tr=0.7
# line; water: open triangle below all four curves and left of that line.
label_pos = {"propane": (1.50, -0.28), "water": (1.50, -3.75)}
label_ha = {"propane": "left", "water": "left"}
for sp in ["propane", "water"]:
    d = SPECIES[sp]
    T07 = 0.7 * d["Tc"]
    Ps07 = Psat_solve("pr", sp, T07)
    logPr07 = np.log10(Ps07 / d["Pc"])
    omega_pr = -1.0 - logPr07
    ax.annotate("", xy=(invTr07, logPr07), xytext=(invTr07, -1.0),
                arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=2.0,
                                 mutation_scale=14, shrinkA=0, shrinkB=0),
                zorder=6)
    ax.plot([invTr07], [logPr07], marker="o", ms=5.5, color=AMBER, zorder=7,
             markeredgecolor="white", markeredgewidth=0.6)
    xt, yt = label_pos[sp]
    ax.annotate(fr"$\omega_{{\rm {sp}}}={omega_pr:.3f}$" "\n"
                fr"(PR here; table {d['w']:.3f})",
                xy=(invTr07, logPr07), xytext=(xt, yt),
                fontsize=10.5, color=AMBER, ha=label_ha[sp], va="center",
                linespacing=1.3,
                arrowprops=dict(arrowstyle="-", color=AMBER, lw=0.9, alpha=0.85))
    print(f"[check] omega({sp}) from PR curve at Tr=0.7: {omega_pr:.4f}  "
          f"(SPEC.md table value {d['w']:.3f})")

ax.set_xlabel(r"$1/T_r$")
ax.set_ylabel(r"$\log_{10}\!\left(P_r^{sat}\right)$")
ax.set_xlim(1.0, 2.95)
ax.set_xticks(np.arange(1.0, 2.01, 0.2))   # ticks pinned to the data range
ax.set_ylim(-4.6, 0.6)

save(fig, "f17_acentric_definition")
