"""figB_common -- small shared helpers for Set B (f17-f34), built on top of
eosstyle.py (not a figure script itself; imported by figB_*.py).

Contains the general cubic-EOS fugacity-coefficient expression and an
isofugacity Psat solver that work for vdW, RK, SRK and PR alike (the
routine mirrors the one written for f16 in figA_16.py, generalised to
delta1 != delta2 for RK/SRK/PR and to the vdW limit delta1=delta2=0
separately), plus a Maxwell equal-area check used to validate the
isofugacity solutions numerically (SPEC.md section 8 for all EOS constants).
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from eosstyle import *
from scipy.optimize import brentq
from scipy.integrate import quad


def ln_phi(model, sp, T, P, Z, Rg=Rb):
    """ln(fugacity coefficient) for a cubic-EOS root Z (general two-parameter
    cubic P = RT/(V-b) - a/((V+d1 b)(V+d2 b)); reduces correctly to the vdW
    form d1=d2=0)."""
    d = SPECIES[sp]
    a, b, d1, d2 = cubic_params(model, d["Tc"], d["Pc"], d["w"], T, Rg)
    A = a * P / (Rg * T) ** 2
    B = b * P / (Rg * T)
    if abs(d1 - d2) < 1e-12:
        # vdW limit (delta1=delta2=0): general term -> A/Z (verified by
        # taking the delta1,delta2 -> 0 limit of the general expression
        # below, and independently by direct integration of the vdW
        # residual Helmholtz energy; standard result ln phi = Z-1-ln(Z-B)-A/Z).
        term = A / Z
    else:
        term = A / (B * (d1 - d2)) * np.log((Z + d1 * B) / (Z + d2 * B))
    return (Z - 1) - np.log(Z - B) - term


def Psat_solve(model, sp, T, Rg=Rb):
    """Solve ln(phi_vap) = ln(phi_liq) (isofugacity) for the saturation
    pressure of species sp at temperature T using the given cubic model.
    Returns np.nan if no bracketing sign change is found (e.g. above Tc)."""
    d = SPECIES[sp]

    def f(P):
        r = Zroots(model, sp, T, P, Rg)
        if len(r) < 2:
            return np.nan
        Zl, Zv = r[0], r[-1]
        return ln_phi(model, sp, T, P, Zv) - ln_phi(model, sp, T, P, Zl)

    Ps = np.geomspace(1e-7 * d["Pc"], 0.999 * d["Pc"], 400)
    vals = np.array([f(P) for P in Ps])
    idx = np.isfinite(vals)
    Ps_f, vals_f = Ps[idx], vals[idx]
    sc = np.where(np.diff(np.sign(vals_f)) != 0)[0]
    if len(sc) == 0:
        return np.nan
    i = sc[0]
    return brentq(f, Ps_f[i], Ps_f[i + 1], xtol=1e-10)


def maxwell_equal_area_check(model, sp, T, Psat, Rg=Rb, verbose=True):
    """Independent check of an isofugacity Psat: integrate P_of_V(V) between
    the liquid and vapour roots at that Psat and confirm the net signed area
    between the EOS isotherm and the horizontal line P = Psat is ~0
    (Maxwell equal-area rule), which is equivalent to equal Gibbs energy.
    Returns the residual area (bar*cm3/mol), printed as the numerical check
    requested by the figure spec."""
    d = SPECIES[sp]
    a, b, d1, d2 = cubic_params(model, d["Tc"], d["Pc"], d["w"], T, Rg)
    r = Zroots(model, sp, T, Psat, Rg)
    Vl = r[0] * Rg * T / Psat
    Vv = r[-1] * Rg * T / Psat

    def integrand(V):
        return P_of_V(model, sp, T, V, Rg) - Psat

    area, err = quad(integrand, Vl, Vv, limit=200)
    if verbose:
        print(f"[check] {model.upper()} {sp} @ T={T:.2f} K: Psat={Psat:.4f} bar, "
              f"Maxwell equal-area residual = {area:.3e} bar*cm3/mol "
              f"(quad err est {err:.1e})")
    return area
