"""f07_maxwell -- enlarged subcritical vdW loop (CO2, T_r = 0.90) with the
equal-area Maxwell construction; mechanically unstable branch (dP/dV > 0)
in rust, metastable branches dashed. View is zoomed to the loop region only
(the near-b divergence of the true isotherm is cropped by the axis limits).
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from eosstyle import *
from scipy.optimize import brentq
from scipy.integrate import quad

sp = "CO2"
d = SPECIES[sp]
Tc, Pc = d["Tc"], d["Pc"]
a, b, _, _ = cubic_params("vdw", Tc, Pc, d["w"], Tc, Rb)

Tr = 0.90
T = Tr*Tc

def Pvdw(V):
    return Rb*T/(V - b) - a/V**2

def dPdV(V):
    return -Rb*T/(V-b)**2 + 2*a/V**3

def three_roots_V(P):
    coeffs = [P, -(P*b + Rb*T), a, -a*b]
    r = np.roots(coeffs)
    r = np.sort(np.real(r[np.abs(np.imag(r)) < 1e-7*max(1, abs(r).max())]))
    return r[r > b*1.001]

def area_diff(P):
    roots = three_roots_V(P)
    if len(roots) < 3:
        return np.nan
    Vl, Vm, Vv = roots
    val, _ = quad(lambda V: Pvdw(V) - P, Vl, Vv, limit=200)
    return val

Ps = np.linspace(0.01*Pc, 0.999*Pc, 400)
vals = np.array([area_diff(P) for P in Ps])
idx = np.isfinite(vals)
Ps, vals = Ps[idx], vals[idx]
sc = np.where(np.diff(np.sign(vals)) != 0)[0][0]
Psat = brentq(area_diff, Ps[sc], Ps[sc+1], xtol=1e-6)

Vl, Vm, Vv = three_roots_V(Psat)
Vsp1 = brentq(dPdV, Vl, Vm)
Vsp2 = brentq(dPdV, Vm, Vv)
# Vl=77.5, Vsp1=92.3, Vm=140.1, Vsp2=196.3, Vv=301.7, Psat=47.8, Pc=73.8 (printed for record)
print(f"Vl={Vl:.1f} Vsp1={Vsp1:.1f} Vm={Vm:.1f} Vsp2={Vsp2:.1f} Vv={Vv:.1f} Psat={Psat:.2f}")

fig, ax = newfig(6.8, 4.6)

XMIN, XMAX = 65.0, 330.0
# the y-axis is carried a little below the loop minimum so that the two
# bottom captions get their own clear band of white space
YMIN, YMAX = 16.0, 65.0

Vgrid = np.linspace(XMIN, XMAX, 3000)
Pgrid = Pvdw(Vgrid)

m_stab_liq = (Vgrid <= Vl)
m_meta_liq = (Vgrid >= Vl) & (Vgrid <= Vsp1)
m_unstab   = (Vgrid >= Vsp1) & (Vgrid <= Vsp2)
m_meta_vap = (Vgrid >= Vsp2) & (Vgrid <= Vv)
m_stab_vap = (Vgrid >= Vv)

ax.plot(Vgrid[m_stab_liq], Pgrid[m_stab_liq], color=NAVY, lw=2.1)
ax.plot(Vgrid[m_stab_vap], Pgrid[m_stab_vap], color=NAVY, lw=2.1)
ax.plot(Vgrid[m_meta_liq], Pgrid[m_meta_liq], color=NAVY, lw=1.8, linestyle="--", alpha=0.9)
ax.plot(Vgrid[m_meta_vap], Pgrid[m_meta_vap], color=NAVY, lw=1.8, linestyle="--", alpha=0.9)
ax.plot(Vgrid[m_unstab], Pgrid[m_unstab], color=RUST, lw=2.2)

# equal-area shading: lower lobe (Vl..Vm, below Psat) teal, upper lobe (Vm..Vv, above Psat) amber
Vloop = Vgrid[(Vgrid >= Vl) & (Vgrid <= Vv)]
Ploop = Pvdw(Vloop)
ax.fill_between(Vloop, Ploop, Psat, where=(Ploop <= Psat), color=TEAL, alpha=0.25, lw=0)
ax.fill_between(Vloop, Ploop, Psat, where=(Ploop >= Psat), color=AMBER, alpha=0.28, lw=0)

ax.axhline(Psat, color=SLATE, lw=1.6, linestyle=":", zorder=2)
ax.plot([Vl, Vv], [Psat, Psat], color=GRAPHITE, lw=2.0, zorder=4)

for V in (Vl, Vm, Vv):
    ax.plot([V], [Psat], marker="o", ms=6.5, color=GRAPHITE, zorder=5, markeredgecolor="white", markeredgewidth=0.8)
ax.plot([Vm], [Psat], marker="o", ms=6.5, color=RUST, zorder=6, markeredgecolor="white", markeredgewidth=0.8)

# --- labels -----------------------------------------------------------
# below the saturation line, in the empty strip between V* and V^vap
ax.text(200.0, Psat - 1.2, r"$P_{sat}(T_r{=}0.90)$", fontsize=10.5, color=GRAPHITE,
        ha="center", va="top")

ax.annotate(r"$V^{liq}$", xy=(Vl, Psat), xytext=(Vl + 33, YMAX - 2.0),
            fontsize=11.0, color=GRAPHITE, ha="center", va="top",
            arrowprops=dict(arrowstyle="-", color=GRAPHITE, lw=1.0))
ax.annotate(r"$V^{vap}$", xy=(Vv, Psat), xytext=(Vv, YMAX - 2.0),
            fontsize=11.0, color=GRAPHITE, ha="center", va="top",
            arrowprops=dict(arrowstyle="-", color=GRAPHITE, lw=1.0))
# V* label placed off the (steep, rust) spinodal branch, in the empty band
# above the saturation line and to the left of that branch
ax.annotate(r"$V^{*}$", xy=(Vm, Psat), xytext=(Vm - 13, Psat + 4.6),
            fontsize=11.0, color=RUST, ha="right", va="center",
            arrowprops=dict(arrowstyle="-", color=RUST, lw=1.0))

# unstable branch label (bottom band, right half, clear of the lower-lobe caption)
ax.annotate("$dP/dV>0$\nmechanically unstable", xy=(122, Pvdw(122)),
            xytext=(152, YMIN + 1.0), fontsize=10.5, color=RUST, ha="left", va="bottom",
            linespacing=1.25, arrowprops=dict(arrowstyle="-", color=RUST, lw=1.0))

# equal-area lobe labels
ax.annotate("equal area\n(lower lobe)", xy=(100, 33), xytext=(76, YMIN + 1.0),
            fontsize=10.5, color="#0b5654", ha="left", va="bottom", linespacing=1.25,
            arrowprops=dict(arrowstyle="-", color=TEAL, lw=1.0))
ax.annotate("equal area\n(upper lobe)", xy=(220, 51.5), xytext=(222, YMAX - 2.0),
            fontsize=10.5, color="#8a6b16", ha="left", va="top", linespacing=1.25,
            arrowprops=dict(arrowstyle="-", color=AMBER, lw=1.0))

# metastable branch notes: label in confirmed-empty corners with a leader
# line to a point on the actual dashed segment
ax.annotate("metastable branches\n(dashed)", xy=(240, Pvdw(240)), xytext=(248, 26.5),
            fontsize=10.5, color=SLATE, ha="left", va="center", linespacing=1.2,
            arrowprops=dict(arrowstyle="-", color=SLATE, lw=0.9))

ax.set_xlim(XMIN, XMAX)
ax.set_ylim(YMIN, YMAX)
ax.set_xlabel(r"Molar volume, $V$ / (cm$^3$ mol$^{-1}$)")
ax.set_ylabel("Pressure, $P$ / bar")

save(fig, "f07_maxwell")
