{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aee6052b",
   "metadata": {},
   "source": [
    "# Online Resource 1: ORC simulation, figures and tables\n",
    "\n",
    "Supplementary code for: Synergistic Oil–Gas–Geothermal Integration and Advanced Exergy Optimization of Organic Rankine Cycles for Refinery Waste Heat Recovery in Nigeria. Running all cells writes Figs. 1–10 to `figures/` and prints every number quoted in Sections 3 to 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74b97a52",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q coolprop numpy pandas matplotlib seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b16b703a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Oil-gas-geothermal ORC integration for Nigerian refineries - simulation and figures.\"\"\"\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "matplotlib.use(\"Agg\")\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from CoolProp.CoolProp import PropsSI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3d11b2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 0. Global settings and reported base-case data\n",
    "# ----------------------------------------------------------------------------\n",
    "T0, P0 = 298.15, 101_325.0                 # dead state (K, Pa)\n",
    "K = lambda c: c + 273.15\n",
    "OUT = \"figures\"; os.makedirs(OUT, exist_ok=True)\n",
    "sns.set_theme(style=\"whitegrid\", context=\"paper\", font_scale=1.0)\n",
    "plt.rcParams.update({\"font.family\": \"sans-serif\", \"font.sans-serif\": [\"Arial\", \"Helvetica\", \"Liberation Sans\", \"DejaVu Sans\"],\n",
    "                     \"font.size\": 9, \"axes.labelsize\": 9, \"xtick.labelsize\": 8, \"ytick.labelsize\": 8, \"legend.fontsize\": 8,\n",
    "                     \"axes.titlesize\": 9, \"savefig.dpi\": 600, \"figure.dpi\": 100})\n",
    "PAL = sns.color_palette(\"deep\", 3)\n",
    "\n",
    "COND_C, ETA_P, ETA_E, DT_PP = 90.0, 0.75, 0.80, 10.0\n",
    "# Reported optima (W, Ed in kW). Tsrc = effective hybrid-carrier inlet temperature (ASSUMED\n",
    "# for Warri/Kaduna; 335 C is the reported Port Harcourt flue-gas temperature).\n",
    "SITES = {\n",
    "    \"Warri\":         dict(fluid=\"IsoButane\", Tev=145.0, W=180.0, Ed=1320.0, psi=0.12, Tsrc=320.0),\n",
    "    \"Port Harcourt\": dict(fluid=\"Isopentane\", Tev=140.0, W=200.0, Ed=1338.0, psi=0.13, Tsrc=335.0),\n",
    "    \"Kaduna\":        dict(fluid=\"n-Pentane\",  Tev=155.0, W=170.0, Ed=1375.0, psi=0.11, Tsrc=300.0),\n",
    "}\n",
    "SITE_COLOR = dict(zip(SITES, PAL))\n",
    "LS = dict(zip(SITES, [\"-\", \"--\", \"-.\"])); MK = dict(zip(SITES, [\"o\", \"s\", \"^\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b556a829",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 1. Thermodynamic cycle model (CoolProp, Helmholtz EOS)\n",
    "# ----------------------------------------------------------------------------\n",
    "def ref_state(f):\n",
    "    return PropsSI(\"H\", \"T\", T0, \"P\", P0, f), PropsSI(\"S\", \"T\", T0, \"P\", P0, f)\n",
    "\n",
    "def cycle(f, Tev_C, pr=1.10, eta_e=ETA_E, eta_p=ETA_P, Tcond_C=COND_C):\n",
    "    \"\"\"Simple ORC: 1 sat. liquid -> 2 pump -> 3 evaporator exit -> 4 expander -> 1 condenser.\n",
    "    Subcritical: saturated vapour at Tev. Transcritical: P_high = pr * p_crit at Tev.\"\"\"\n",
    "    Tc, pc = PropsSI(\"Tcrit\", f), PropsSI(\"pcrit\", f)\n",
    "    Tev, Tcd = K(Tev_C), K(Tcond_C)\n",
    "    P_lo = PropsSI(\"P\", \"T\", Tcd, \"Q\", 0, f)\n",
    "    h1, s1 = PropsSI(\"H\", \"T\", Tcd, \"Q\", 0, f), PropsSI(\"S\", \"T\", Tcd, \"Q\", 0, f)\n",
    "    if Tev < Tc - 1.0:\n",
    "        regime, P_hi = \"subcritical\", PropsSI(\"P\", \"T\", Tev, \"Q\", 1, f)\n",
    "        h3, s3 = PropsSI(\"H\", \"T\", Tev, \"Q\", 1, f), PropsSI(\"S\", \"T\", Tev, \"Q\", 1, f)\n",
    "    else:\n",
    "        regime, P_hi = \"transcritical\", pr * pc\n",
    "        h3, s3 = PropsSI(\"H\", \"T\", Tev, \"P\", P_hi, f), PropsSI(\"S\", \"T\", Tev, \"P\", P_hi, f)\n",
    "    h2 = h1 + (PropsSI(\"H\", \"P\", P_hi, \"S\", s1, f) - h1) / eta_p\n",
    "    s2, T2 = PropsSI(\"S\", \"P\", P_hi, \"H\", h2, f), PropsSI(\"T\", \"P\", P_hi, \"H\", h2, f)\n",
    "    h4 = h3 - eta_e * (h3 - PropsSI(\"H\", \"P\", P_lo, \"S\", s3, f))\n",
    "    s4, T4 = PropsSI(\"S\", \"P\", P_lo, \"H\", h4, f), PropsSI(\"T\", \"P\", P_lo, \"H\", h4, f)\n",
    "    h0, s0 = ref_state(f)\n",
    "    ex = lambda h, s: (h - h0) - T0 * (s - s0)\n",
    "    c = dict(f=f, regime=regime, P_hi=P_hi, P_lo=P_lo, Tc=Tc, pc=pc, Tev=Tev, Tcd=Tcd, T2=T2, T4=T4,\n",
    "             h=[h1, h2, h3, h4], s=[s1, s2, s3, s4])\n",
    "    c.update(w_net=((h3 - h4) - (h2 - h1)) / 1e3, q_in=(h3 - h2) / 1e3, q_cond=(h4 - h1) / 1e3,\n",
    "             ex=[ex(h1, s1), ex(h2, s2), ex(h3, s3), ex(h4, s4)], w_p=(h2 - h1) / 1e3)\n",
    "    c[\"eta_th\"] = c[\"w_net\"] / c[\"q_in\"]\n",
    "    return c\n",
    "\n",
    "def pinch_margin(c, m, C_s, Tout, n=80):\n",
    "    \"\"\"Minimum source-WF temperature difference along the counter-current evaporator (K).\"\"\"\n",
    "    f, P = c[\"f\"], c[\"P_hi\"]\n",
    "    hs = np.linspace(c[\"h\"][1], c[\"h\"][2], n)\n",
    "    Twf = np.array([PropsSI(\"T\", \"P\", P, \"H\", h, f) for h in hs])\n",
    "    Tsrc = Tout + m * (hs - c[\"h\"][1]) / 1e3 / C_s\n",
    "    return float(np.min(Tsrc - Twf))\n",
    "\n",
    "def site_balance(f, Tev_C, Tsrc_C, C_s, W_kW=None, Q_kW=None, **kw):\n",
    "    \"\"\"Component-level exergy balance with the source-side stack loss made explicit (kW).\"\"\"\n",
    "    c = cycle(f, Tev_C, **kw)\n",
    "    m = W_kW / c[\"w_net\"] if W_kW is not None else Q_kW / c[\"q_in\"]           # kg/s\n",
    "    Q = m * c[\"q_in\"]; Tin = K(Tsrc_C); Tout = Tin - Q / C_s\n",
    "    exf = lambda T: C_s * ((T - T0) - T0 * np.log(T / T0))\n",
    "    Ex_src, Ex_out = exf(Tin), exf(Tout)\n",
    "    s, ex = c[\"s\"], c[\"ex\"]\n",
    "    # ex[] and s[] are in J/kg; m in kg/s -> divide by 1e3 to obtain kW.\n",
    "    comp = {\"Stack loss\": Ex_out,\n",
    "            \"Evaporator\": (Ex_src - Ex_out) - m * (ex[2] - ex[1]) / 1e3,\n",
    "            \"Expander\": m * T0 * (s[3] - s[2]) / 1e3,\n",
    "            \"Pump\": m * T0 * (s[1] - s[0]) / 1e3,\n",
    "            \"Condenser\": m * (ex[3] - ex[0]) / 1e3}\n",
    "    W = m * c[\"w_net\"]\n",
    "    return dict(c=c, m=m, Q=Q, W=W, Tin=Tin, Tout=Tout, Ex_src=Ex_src, comp=comp, Ed=sum(comp.values()),\n",
    "                Ed_cycle=sum(v for k_, v in comp.items() if k_ != \"Stack loss\"),\n",
    "                psi=W / Ex_src, pinch=pinch_margin(c, m, C_s, Tout), Q_cond=m * c[\"q_cond\"])\n",
    "\n",
    "def calibrate(name):\n",
    "    \"\"\"Closure: Ex_src = W + Ed; C_s from Ex_src; flow from W_net.\"\"\"\n",
    "    d = SITES[name]; Tin = K(d[\"Tsrc\"])\n",
    "    Ex_src = d[\"W\"] + d[\"Ed\"]\n",
    "    C_s = Ex_src / ((Tin - T0) - T0 * np.log(Tin / T0))                      # kW/K\n",
    "    return Ex_src, C_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21f8ee7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 2. Advanced exergy: endogenous/exogenous and avoidable/unavoidable\n",
    "# ----------------------------------------------------------------------------\n",
    "def advanced_exergy(name):\n",
    "    d = SITES[name]; _, C_s = calibrate(name)\n",
    "    real = site_balance(d[\"fluid\"], d[\"Tev\"], d[\"Tsrc\"], C_s, W_kW=d[\"W\"])\n",
    "    # Endogenous evaporator destruction: all other components ideal (eta = 1) at the same net product W.\n",
    "    ideal = site_balance(d[\"fluid\"], d[\"Tev\"], d[\"Tsrc\"], C_s, W_kW=real[\"W\"], eta_e=1.0, eta_p=1.0)\n",
    "    ed_evap, ed_endo = real[\"comp\"][\"Evaporator\"], ideal[\"comp\"][\"Evaporator\"]\n",
    "    # Unavoidable conditions: best-practice machines, Tev at the 160 C window limit, same duty.\n",
    "    un = site_balance(d[\"fluid\"], 160.0, d[\"Tsrc\"], C_s, Q_kW=real[\"Q\"], eta_e=0.90, eta_p=0.85)\n",
    "    ed_un = {k_: v / un[\"W\"] * real[\"W\"] for k_, v in un[\"comp\"].items()}     # (Ed/P)_UN * P_real\n",
    "    rows = [dict(Component=k_, Ed_kW=real[\"comp\"][k_], Ed_UN_kW=min(ed_un[k_], real[\"comp\"][k_]),\n",
    "                 Ed_AV_kW=max(real[\"comp\"][k_] - ed_un[k_], 0.0)) for k_ in real[\"comp\"]]\n",
    "    return real, dict(endo=ed_endo, exo=ed_evap - ed_endo), pd.DataFrame(rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51bbffc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 3. NSGA-II (compact implementation): max W_net, min cycle exergy destruction\n",
    "# ----------------------------------------------------------------------------\n",
    "def nds(F):\n",
    "    n = len(F); dom = [[] for _ in range(n)]; cnt = np.zeros(n, int); fronts = [[]]\n",
    "    for p in range(n):\n",
    "        for q in range(n):\n",
    "            if np.all(F[p] <= F[q]) and np.any(F[p] < F[q]): dom[p].append(q)\n",
    "            elif np.all(F[q] <= F[p]) and np.any(F[q] < F[p]): cnt[p] += 1\n",
    "        if cnt[p] == 0: fronts[0].append(p)\n",
    "    i = 0\n",
    "    while fronts[i]:\n",
    "        nxt = []\n",
    "        for p in fronts[i]:\n",
    "            for q in dom[p]:\n",
    "                cnt[q] -= 1\n",
    "                if cnt[q] == 0: nxt.append(q)\n",
    "        fronts.append(nxt); i += 1\n",
    "    return fronts[:-1]\n",
    "\n",
    "def crowd(F, fr):\n",
    "    d = np.zeros(len(F))\n",
    "    for j in range(F.shape[1]):\n",
    "        o = sorted(fr, key=lambda i: F[i, j]); d[o[0]] = d[o[-1]] = np.inf\n",
    "        rng = F[o[-1], j] - F[o[0], j] or 1.0\n",
    "        for a in range(1, len(o) - 1): d[o[a]] += (F[o[a + 1], j] - F[o[a - 1], j]) / rng\n",
    "    return d\n",
    "\n",
    "def nsga2(fun, lb, ub, pop=60, gen=60, seed=7, eta_c=15, eta_m=20):\n",
    "    rng = np.random.default_rng(seed); lb, ub = np.array(lb), np.array(ub)\n",
    "    X = lb + rng.random((pop, len(lb))) * (ub - lb); F = np.array([fun(x) for x in X]); archive = [(X.copy(), F.copy())]\n",
    "    for _ in range(gen):\n",
    "        fr = nds(F); rank = np.zeros(pop, int); cd = np.zeros(pop)\n",
    "        for r, f_ in enumerate(fr): rank[f_] = r; cd[f_] = crowd(F, f_)[f_]\n",
    "        pick = lambda: min(rng.integers(0, pop, 2), key=lambda i: (rank[i], -cd[i]))\n",
    "        kids = []\n",
    "        while len(kids) < pop:\n",
    "            p1, p2 = X[pick()], X[pick()]; c1, c2 = p1.copy(), p2.copy()\n",
    "            if rng.random() < 0.9:\n",
    "                u = rng.random(len(lb)); b = np.where(u <= .5, (2 * u) ** (1 / (eta_c + 1)), (1 / (2 * (1 - u))) ** (1 / (eta_c + 1)))\n",
    "                c1, c2 = .5 * ((1 + b) * p1 + (1 - b) * p2), .5 * ((1 - b) * p1 + (1 + b) * p2)\n",
    "            for c_ in (c1, c2):\n",
    "                m_ = rng.random(len(lb)) < 1 / len(lb); u = rng.random(len(lb))\n",
    "                dl = np.where(u < .5, (2 * u) ** (1 / (eta_m + 1)) - 1, 1 - (2 * (1 - u)) ** (1 / (eta_m + 1)))\n",
    "                c_[m_] += dl[m_] * (ub - lb)[m_]; kids.append(np.clip(c_, lb, ub))\n",
    "        Xk = np.array(kids[:pop]); Fk = np.array([fun(x) for x in Xk]); archive.append((Xk.copy(), Fk.copy()))\n",
    "        XA, FA = np.vstack([X, Xk]), np.vstack([F, Fk]); nxt = []\n",
    "        for f_ in nds(FA):\n",
    "            if len(nxt) + len(f_) <= pop: nxt += f_\n",
    "            else:\n",
    "                cdd = crowd(FA, f_)[f_]; nxt += [f_[i] for i in np.argsort(-cdd)[:pop - len(nxt)]]; break\n",
    "        X, F = XA[nxt], FA[nxt]\n",
    "    XA = np.vstack([a[0] for a in archive]); FA = np.vstack([a[1] for a in archive])\n",
    "    return X[nds(F)[0]], F[nds(F)[0]], XA, FA\n",
    "\n",
    "def run_optimisation(name, seed=7):\n",
    "    d = SITES[name]; _, C_s = calibrate(name)\n",
    "    ref = site_balance(d[\"fluid\"], d[\"Tev\"], d[\"Tsrc\"], C_s, W_kW=d[\"W\"])\n",
    "    def fun(x):                                   # x = [Tev (C), pr (P_hi/p_crit), phi (duty fraction)]\n",
    "        try:\n",
    "            r = site_balance(d[\"fluid\"], x[0], d[\"Tsrc\"], C_s, Q_kW=x[2] * ref[\"Q\"], pr=x[1])\n",
    "            if r[\"pinch\"] < DT_PP or r[\"Tout\"] < K(100.0): return np.array([1e3, 1e5])\n",
    "            return np.array([-r[\"W\"], r[\"Ed_cycle\"]])\n",
    "        except Exception:\n",
    "            return np.array([1e3, 1e5])\n",
    "    Xp, Fp, XA, FA = nsga2(fun, [140.0, 1.05, 0.5], [160.0, 1.40, 1.3], seed=seed)\n",
    "    return dict(ref=ref, Xp=Xp, Fp=Fp, XA=XA, FA=FA, C_s=C_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38dc3de3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 4. Techno-economic and environmental model\n",
    "# ----------------------------------------------------------------------------\n",
    "CEPCI, CEPCI_2001 = 800.0, 397.0\n",
    "I_DISC, N_LIFE, CF, CF_SA, OM_FRAC = 0.12, 20, 0.90, 0.75, 0.04   # CF_SA: stand-alone (no brine stabilisation)\n",
    "PRICE_BASE, CAPEX_WELL, EXP_USD_KW = 0.20, 0.30e6, 900.0                # USD/kWh avoided cost; well/brine-loop retrofit\n",
    "EF_ELEC, EF_HEAT, F_USE = 0.60, 0.238, 0.25           # tCO2/MWh_e displaced; tCO2/MWh_th; usable condenser-heat share\n",
    "GAS_USD_MWH = 10.0\n",
    "E_BASE = {\"Warri\": 10_500.0, \"Port Harcourt\": 10_000.0, \"Kaduna\": 9_500.0}   # t/yr (PLACEHOLDER baselines)\n",
    "UPLIFT = 0.25                                          # co-production uplift (20-30 % reported range)\n",
    "\n",
    "def turton(K1, K2, K3, X): return 10 ** (K1 + K2 * np.log10(X) + K3 * np.log10(X) ** 2)\n",
    "\n",
    "def capex(res):\n",
    "    lm = lambda a, b: (a - b) / np.log(a / b) if abs(a - b) > 1e-6 else a\n",
    "    A_ev = res[\"Q\"] / (0.12 * lm(res[\"Tin\"] - res[\"c\"][\"Tev\"], res[\"Tout\"] - res[\"c\"][\"T2\"]))  # gas-side U, kW/m2K\n",
    "    A_cd = res[\"Q_cond\"] / (0.45 * 25.0)                                                       # water-preheat duty\n",
    "    C_hx = (turton(4.3247, -0.3030, 0.1634, A_ev) * (1.63 + 1.66 * 2.73)\n",
    "            + turton(4.3247, -0.3030, 0.1634, A_cd) * (1.63 + 1.66)) * CEPCI / CEPCI_2001\n",
    "    C_pp = turton(3.3892, 0.0536, 0.1538, res[\"m\"] * res[\"c\"][\"w_p\"]) * 3.3 * CEPCI / CEPCI_2001\n",
    "    C_ex = EXP_USD_KW * res[\"W\"]                                                               # packaged expander-generator\n",
    "    orc = (C_hx + C_pp + C_ex) * 1.18                                                          # 18 % contingency and fee\n",
    "    return orc\n",
    "\n",
    "def econ(res, price=PRICE_BASE, uplift=UPLIFT, hybrid=True):\n",
    "    scale = 1.0 if hybrid else 1.0 / (1.0 + uplift)                    # stand-alone = hybrid/(1+u)\n",
    "    W = res[\"W\"] * scale\n",
    "    capex_tot = capex(res) * scale ** 0.7 + (CAPEX_WELL if hybrid else 0.0)\n",
    "    cf = CF if hybrid else CF_SA\n",
    "    E = W * 8760 * cf                                                  # kWh/yr\n",
    "    heat_MWh = res[\"Q_cond\"] * scale * 8760 * cf * F_USE / 1e3\n",
    "    rev = E * price + heat_MWh * GAS_USD_MWH\n",
    "    net = rev - OM_FRAC * capex_tot\n",
    "    crf = I_DISC * (1 + I_DISC) ** N_LIFE / ((1 + I_DISC) ** N_LIFE - 1)\n",
    "    return dict(capex=capex_tot, payback=capex_tot / net if net > 0 else np.inf,\n",
    "                lcoe=(capex_tot * crf + OM_FRAC * capex_tot) / E, E_MWh=E / 1e3, heat_MWh=heat_MWh)\n",
    "\n",
    "def co2_reduction(name, res, uplift=UPLIFT, hybrid=True):\n",
    "    scale = (1.0 + uplift) / (1.0 + UPLIFT) if hybrid else 1.0 / (1.0 + UPLIFT)\n",
    "    cf = CF if hybrid else CF_SA\n",
    "    mwh = res[\"W\"] * scale * 8760 * cf / 1e3\n",
    "    heat = res[\"Q_cond\"] * scale * 8760 * cf * F_USE / 1e3\n",
    "    return 100.0 * (mwh * EF_ELEC + heat * EF_HEAT) / E_BASE[name]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0069c8a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 5. Figures\n",
    "# ----------------------------------------------------------------------------\n",
    "def ts_diagram(name, num):\n",
    "    d = SITES[name]; f = d[\"fluid\"]; c = cycle(f, d[\"Tev\"], pr=1.10)\n",
    "    Tt = np.linspace(K(30.0), c[\"Tc\"] - 0.05, 250)\n",
    "    sL = [PropsSI(\"S\", \"T\", t, \"Q\", 0, f) / 1e3 for t in Tt]; sV = [PropsSI(\"S\", \"T\", t, \"Q\", 1, f) / 1e3 for t in Tt]\n",
    "    fig, ax = plt.subplots(figsize=(6.2, 4.6))\n",
    "    ax.plot(sL, Tt - 273.15, \"k\", lw=1.6); ax.plot(sV, Tt - 273.15, \"k\", lw=1.6, label=\"Saturation dome\")\n",
    "    s1, s2, s3, s4 = [x / 1e3 for x in c[\"s\"]]; h2, h3 = c[\"h\"][1], c[\"h\"][2]\n",
    "    hs = np.linspace(h2, h3, 60)\n",
    "    Ts = [PropsSI(\"T\", \"P\", c[\"P_hi\"], \"H\", h, f) for h in hs]; ss = [PropsSI(\"S\", \"P\", c[\"P_hi\"], \"H\", h, f) / 1e3 for h in hs]\n",
    "    ax.plot([s1, s2], [K(COND_C) - 273.15, c[\"T2\"] - 273.15], color=PAL[0], lw=2)\n",
    "    ax.plot(ss, np.array(Ts) - 273.15, color=\"firebrick\", lw=2.2, label=\"Heat addition\")\n",
    "    ax.plot([s3, s4], [c[\"Tev\"] - 273.15, c[\"T4\"] - 273.15], color=\"darkgreen\", lw=2.2, label=\"Expansion\")\n",
    "    sg = PropsSI(\"S\", \"P\", c[\"P_lo\"], \"Q\", 1, f) / 1e3\n",
    "    ax.plot([s4, sg, s1], [c[\"T4\"] - 273.15, COND_C, COND_C], color=\"navy\", lw=2.2, label=\"Condensation\")\n",
    "    for i, (s_, T_) in enumerate([(s1, COND_C), (s2, c[\"T2\"] - 273.15), (s3, c[\"Tev\"] - 273.15), (s4, c[\"T4\"] - 273.15)], 1):\n",
    "        ax.plot(s_, T_, \"o\", color=\"k\"); ax.annotate(str(i), (s_, T_), textcoords=\"offset points\", xytext=(6, 4))\n",
    "    ax.axhline(c[\"Tc\"] - 273.15, ls=\":\", color=\"grey\"); ax.text(ax.get_xlim()[0], c[\"Tc\"] - 273.15 + 2, f\"$T_c$ = {c['Tc']-273.15:.1f} °C\", fontsize=9, color=\"grey\")\n",
    "    ax.set(xlabel=\"Specific entropy s (kJ/kg K)\", ylabel=\"Temperature (°C)\")\n",
    "    ax.set_ylim(40, c[\"Tc\"] - 273.15 + 25); ax.legend(loc=\"lower left\", fontsize=8); fig.tight_layout(); fig.savefig(f\"{OUT}/fig{num:02d}_ts_{name.replace(' ', '_')}.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "def make_all():\n",
    "    results, adv, opt, tables = {}, {}, {}, {}\n",
    "    for i, n in enumerate(SITES, 1):\n",
    "        ts_diagram(n, i)\n",
    "        results[n], adv[n], tables[n] = advanced_exergy(n)\n",
    "        opt[n] = run_optimisation(n)\n",
    "\n",
    "    # Fig 4: base-case exergy efficiency\n",
    "    fig, ax = plt.subplots(figsize=(6, 4.2))\n",
    "    vals = [100 * SITES[n][\"psi\"] for n in SITES] + [28.6]; labs = list(SITES) + [\"Warri refinery\\n(plant base case)\"]\n",
    "    bars = ax.bar(labs, vals, color=PAL + [\"lightgrey\"], edgecolor=\"k\"); [b.set_hatch(h) for b, h in zip(bars, [\"\", \"..\", \"xx\", \"//\"])]\n",
    "    for b, v in zip(bars, vals): ax.text(b.get_x() + b.get_width() / 2, v + 0.6, f\"{v:.1f}\", ha=\"center\")\n",
    "    ax.set(ylabel=\"Exergy efficiency (%)\"); fig.tight_layout(); fig.savefig(f\"{OUT}/fig04_exergy_efficiency.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "    # Fig 5: component breakdown (Warri)\n",
    "    comp = results[\"Warri\"][\"comp\"]; tot = sum(comp.values()); fig, ax = plt.subplots(figsize=(6.4, 4.6))\n",
    "    wedges, _ = ax.pie(comp.values(), startangle=90, counterclock=False, colors=sns.color_palette(\"Set2\", len(comp)), wedgeprops=dict(edgecolor=\"w\"))\n",
    "    [w_.set_hatch(h) for w_, h in zip(wedges, [\"\", \"//\", \"..\", \"xx\", \"\\\\\\\\\"])]\n",
    "    ax.legend(wedges, [f\"{k_}: {v:.0f} kW ({100 * v / tot:.1f}%)\" for k_, v in comp.items()], loc=\"center left\", bbox_to_anchor=(1.0, 0.5), frameon=False)\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig05_warri_breakdown.png\", dpi=600, bbox_inches=\"tight\"); plt.close(fig)\n",
    "\n",
    "    # Fig 6: endogenous vs exogenous (evaporator, Warri)\n",
    "    a = adv[\"Warri\"]; fig, ax = plt.subplots(figsize=(5, 4.6))\n",
    "    ax.pie([a[\"endo\"], a[\"exo\"]], labels=[f\"Endogenous\\n{a['endo']:.0f} kW\", f\"Exogenous\\n{a['exo']:.0f} kW\"], autopct=\"%1.1f%%\",\n",
    "           colors=[\"#4c72b0\", \"#dd8452\"], startangle=90, wedgeprops=dict(edgecolor=\"w\"))\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig06_endo_exo.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "    # Fig 7: Pareto fronts\n",
    "    fig, ax = plt.subplots(figsize=(6.6, 4.8))\n",
    "    for n in SITES:\n",
    "        o = opt[n]; ok = o[\"FA\"][:, 0] < 0\n",
    "        ax.scatter(-o[\"FA\"][ok, 0], o[\"FA\"][ok, 1], s=6, alpha=.15, color=SITE_COLOR[n])\n",
    "        idx = np.argsort(o[\"Fp\"][:, 0]); ax.plot(-o[\"Fp\"][idx, 0], o[\"Fp\"][idx, 1], ls=LS[n], marker=MK[n], ms=3.5, lw=1.4, color=SITE_COLOR[n], label=f\"{n} front\")\n",
    "        r = o[\"ref\"]; ax.plot(r[\"W\"], r[\"Ed_cycle\"], \"*\", ms=14, mec=\"k\", color=SITE_COLOR[n])\n",
    "    ax.plot([], [], \"k*\", ms=11, label=\"Reported operating point\")\n",
    "    ax.set(xlabel=\"Net power output (kW)\", ylabel=\"Cycle exergy destruction and loss (kW)\"); ax.legend(fontsize=8)\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig07_pareto.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "    # Fig 8: parametric effect of evaporator temperature at constant heat duty\n",
    "    fig, ax = plt.subplots(figsize=(6.2, 4.4)); Tgrid = np.linspace(140, 160, 21); rows8 = []\n",
    "    for n, d in SITES.items():\n",
    "        _, C_s = calibrate(n); ref = results[n]; W = []\n",
    "        for T in Tgrid:\n",
    "            W.append(site_balance(d[\"fluid\"], T, d[\"Tsrc\"], C_s, Q_kW=ref[\"Q\"])[\"W\"])\n",
    "        ax.plot(Tgrid, W, lw=2, ls=LS[n], color=SITE_COLOR[n], label=f\"{n} ({d['fluid']})\"); ax.plot(d[\"Tev\"], ref[\"W\"], \"o\", color=SITE_COLOR[n], mec=\"k\")\n",
    "        rows8.append((n, W[0], W[-1]))\n",
    "    ax.set(xlabel=\"Evaporator temperature (°C)\", ylabel=\"Net power (kW)\"); ax.legend(fontsize=8)\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig08_temperature_sweep.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "    # Fig 9: CO2 reductions\n",
    "    fig, ax = plt.subplots(figsize=(6.2, 4.4)); x = np.arange(3); w = .36; rows9 = {}\n",
    "    for j, hyb in enumerate([False, True]):\n",
    "        mid = [co2_reduction(n, results[n], hybrid=hyb) for n in SITES]\n",
    "        lo = [co2_reduction(n, results[n], uplift=0.20, hybrid=hyb) for n in SITES]; hi = [co2_reduction(n, results[n], uplift=0.30, hybrid=hyb) for n in SITES]\n",
    "        err = np.array([np.array(mid) - np.array(lo), np.array(hi) - np.array(mid)]) if hyb else None\n",
    "        ax.bar(x + (j - .5) * w, mid, w, yerr=err, capsize=4, label=\"Hybrid co-production\" if hyb else \"Stand-alone waste heat\", edgecolor=\"k\", hatch=\"//\" if hyb else \"\", color=PAL[1] if hyb else \"lightgrey\")\n",
    "        rows9[hyb] = mid\n",
    "    ax.axhspan(15, 25, color=\"green\", alpha=.08); ax.set_xticks(x, list(SITES))\n",
    "    ax.set(ylabel=\"CO$_2$ emission reduction (%)\"); ax.legend(fontsize=8)\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig09_co2.png\", dpi=600); plt.close(fig)\n",
    "\n",
    "    # Fig 10: payback sensitivity\n",
    "    fig, ax = plt.subplots(figsize=(6.2, 4.4)); prices = np.linspace(0.06, 0.30, 25)\n",
    "    for n in SITES: ax.plot(prices, [econ(results[n], price=p)[\"payback\"] for p in prices], lw=2, ls=LS[n], color=SITE_COLOR[n], label=n)\n",
    "    ax.axhline(5, ls=\"--\", color=\"k\"); ax.set(ylim=(0, 15), xlabel=\"Avoided electricity price (USD/kWh)\", ylabel=\"Simple payback (years)\"); ax.legend(fontsize=8)\n",
    "    fig.tight_layout(); fig.savefig(f\"{OUT}/fig10_payback.png\", dpi=600); plt.close(fig)\n",
    "    return results, adv, tables, opt, rows8, rows9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b1ac401",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------------------------------------------------------\n",
    "# 6. Reproducible numerical report (every number quoted in the manuscript)\n",
    "# ----------------------------------------------------------------------------\n",
    "def brine_share(name, T_br_C, dT=DT_PP):\n",
    "    \"\"\"Largest share of the working-fluid heating duty that brine at T_br_C can supply (sensible preheat only).\"\"\"\n",
    "    d = SITES[name]; _, C_s = calibrate(name)\n",
    "    r = site_balance(d[\"fluid\"], d[\"Tev\"], d[\"Tsrc\"], C_s, W_kW=d[\"W\"]); c = r[\"c\"]; f = c[\"f\"]\n",
    "    T_pre = K(T_br_C) - dT\n",
    "    if c[\"regime\"] == \"subcritical\": T_pre = min(T_pre, PropsSI(\"T\", \"P\", c[\"P_hi\"], \"Q\", 0, f))\n",
    "    if T_pre <= c[\"T2\"]: return 0.0, 0.0, 0.0\n",
    "    h_pre = (PropsSI(\"H\", \"P\", c[\"P_hi\"], \"Q\", 0, f) if (c[\"regime\"] == \"subcritical\" and T_pre >= PropsSI(\"T\", \"P\", c[\"P_hi\"], \"Q\", 0, f) - 1e-6)\n",
    "             else PropsSI(\"H\", \"P\", c[\"P_hi\"], \"T\", T_pre, f))\n",
    "    share = (h_pre - c[\"h\"][1]) / (c[\"h\"][2] - c[\"h\"][1]); Q_br = share * r[\"Q\"]\n",
    "    m_br = Q_br / (4.0 * max(K(T_br_C) - (c[\"T2\"] + dT), 1.0))      # kg/s, brine cp = 4.0 kJ/kgK\n",
    "    return share, Q_br, m_br * 3.6                                   # m3/h at rho ~ 1000 kg/m3\n",
    "\n",
    "def report(results, adv, tables, opt):\n",
    "    pd.set_option(\"display.width\", 220); pd.set_option(\"display.max_columns\", 40)\n",
    "    print(\"\\n== Table 2: fluid properties ==\")\n",
    "    for f in [\"IsoButane\", \"Isopentane\", \"n-Pentane\"]:\n",
    "        print(f, round(PropsSI(\"Tcrit\", f) - 273.15, 2), round(PropsSI(\"pcrit\", f) / 1e6, 3),\n",
    "              round(PropsSI(\"P\", \"T\", K(COND_C), \"Q\", 0, f) / 1e6, 3), round(PropsSI(\"M\", f) * 1e3, 2))\n",
    "    print(\"Carnot limit 145/90 C:\", round(1 - K(COND_C) / K(145.0), 4), \"| Tr(145 C) isobutane:\", round(K(145.0) / PropsSI(\"Tcrit\", \"IsoButane\"), 4),\n",
    "          \"| Tr,cond:\", round(K(COND_C) / PropsSI(\"Tcrit\", \"IsoButane\"), 4))\n",
    "    print(\"\\n== Table 3: base case ==\")\n",
    "    rows = []\n",
    "    for n, r in results.items():\n",
    "        e_h, e_s = econ(r), econ(r, hybrid=False); Tin = r[\"Tin\"]\n",
    "        fac = ((Tin - T0) - T0 * np.log(Tin / T0)) / (Tin - T0)\n",
    "        rows.append(dict(Site=n, Regime=r[\"c\"][\"regime\"], W=r[\"W\"], Ex_F=r[\"Ex_src\"], Ed=r[\"Ed\"], psi=r[\"psi\"], Ed_check=r[\"Ex_src\"] - r[\"W\"],\n",
    "                         m_wf=r[\"m\"], Q_in=r[\"Q\"], Q_cond=r[\"Q_cond\"], eta_th=100 * r[\"c\"][\"eta_th\"], T_out=r[\"Tout\"] - 273.15, pinch=r[\"pinch\"],\n",
    "                         P_hi=r[\"c\"][\"P_hi\"] / 1e6, P_lo=r[\"c\"][\"P_lo\"] / 1e6, exp_ratio=r[\"c\"][\"P_hi\"] / r[\"c\"][\"P_lo\"],\n",
    "                         Ex_factor=fac, H_above_T0_MW=r[\"Ex_src\"] / fac / 1e3, C_s_kW_per_K=calibrate(n)[1],\n",
    "                         carnot_pct=100 * (1 - K(COND_C) / K(SITES[n][\"Tev\"])),\n",
    "                         CAPEX_MUSD=e_h[\"capex\"] / 1e6, USD_per_kW=e_h[\"capex\"] / r[\"W\"], payback=e_h[\"payback\"], payback_SA=e_s[\"payback\"],\n",
    "                         LCOE=e_h[\"lcoe\"], E_hyb_MWh=e_h[\"E_MWh\"], E_SA_MWh=e_s[\"E_MWh\"],\n",
    "                         CO2_hyb=co2_reduction(n, r), CO2_hyb_20=co2_reduction(n, r, 0.20), CO2_hyb_30=co2_reduction(n, r, 0.30),\n",
    "                         CO2_SA=co2_reduction(n, r, hybrid=False), CO2_t=co2_reduction(n, r) * E_BASE[n] / 100))\n",
    "    print(pd.DataFrame(rows).set_index(\"Site\").round(3).T)\n",
    "    print(\"\\n== Component breakdown (kW) and advanced exergy ==\")\n",
    "    for n in SITES:\n",
    "        comp = results[n][\"comp\"]; tot = sum(comp.values())\n",
    "        print(n, {k_: (round(v, 1), f\"{100 * v / tot:.1f}%\") for k_, v in comp.items()})\n",
    "        print(\"   endo/exo evaporator:\", {k_: round(v, 1) for k_, v in adv[n].items()}, \"endo share %:\", round(100 * adv[n][\"endo\"] / (adv[n][\"endo\"] + adv[n][\"exo\"]), 1))\n",
    "        t = tables[n].copy(); t[\"AV_share_%\"] = 100 * t[\"Ed_AV_kW\"] / t[\"Ed_kW\"]; print(t.round(1).to_string(index=False)); print(\"   totals:\", t[[\"Ed_kW\", \"Ed_UN_kW\", \"Ed_AV_kW\"]].sum().round(1).to_dict())\n",
    "        hx = t.set_index(\"Component\").loc[[\"Evaporator\", \"Condenser\"], [\"Ed_kW\", \"Ed_UN_kW\"]].sum(); print(\"   HX share of Ed / of UN %:\", round(100 * hx[\"Ed_kW\"] / t[\"Ed_kW\"].sum(), 1), round(100 * hx[\"Ed_UN_kW\"] / t[\"Ed_UN_kW\"].sum(), 1))\n",
    "        print(\"   LMTD evaporator (K):\", round(((r_ := results[n])[\"Tin\"] - r_[\"c\"][\"Tev\"] - (r_[\"Tout\"] - r_[\"c\"][\"T2\"])) / np.log((r_[\"Tin\"] - r_[\"c\"][\"Tev\"]) / (r_[\"Tout\"] - r_[\"c\"][\"T2\"])), 1))\n",
    "    print(\"\\n== Pareto fronts ==\")\n",
    "    for n in SITES:\n",
    "        o = opt[n]; W = -o[\"Fp\"][:, 0]; Ed = o[\"Fp\"][:, 1]; i = np.argsort(W); ref = o[\"ref\"]\n",
    "        gap = float(np.interp(ref[\"W\"], W[i], Ed[i]))\n",
    "        print(n, \"W range\", W.min().round(1), W.max().round(1), \"| Tev range\", o[\"Xp\"][:, 0].min().round(1), o[\"Xp\"][:, 0].max().round(1),\n",
    "              \"| pr range\", o[\"Xp\"][:, 1].min().round(2), o[\"Xp\"][:, 1].max().round(2), \"| ref Ed_cyc\", round(ref[\"Ed_cycle\"], 1), \"front Ed at ref W\", round(gap, 1),\n",
    "              \"gap %\", round(100 * (ref[\"Ed_cycle\"] - gap) / ref[\"Ed_cycle\"], 1), \"| n front\", len(W))\n",
    "    print(\"\\n== Fig 8 sweep (W at 140, 145, 160 C; constant duty) ==\")\n",
    "    for n, d in SITES.items():\n",
    "        _, C_s = calibrate(n); Q = results[n][\"Q\"]\n",
    "        print(n, [round(site_balance(d[\"fluid\"], T, d[\"Tsrc\"], C_s, Q_kW=Q)[\"W\"], 1) for T in (140.0, 145.0, 150.0, 155.0, 160.0)])\n",
    "    print(\"\\n== Table 5: fluid cross-check at reported Tev and duty (kW) ==\")\n",
    "    for n, d in SITES.items():\n",
    "        _, C_s = calibrate(n)\n",
    "        print(n, {f: round(site_balance(f, d[\"Tev\"], d[\"Tsrc\"], C_s, Q_kW=results[n][\"Q\"])[\"W\"], 1) for f in [\"IsoButane\", \"Isopentane\", \"n-Pentane\"]})\n",
    "    print(\"\\n== Payback sensitivity ==\")\n",
    "    for n in SITES:\n",
    "        r = results[n]; be = next((p for p in np.arange(0.05, 0.40, 0.001) if econ(r, price=p)[\"payback\"] <= 5.0), None)\n",
    "        print(n, \"payback @0.12/0.14/0.20/0.30:\", [round(econ(r, price=p)[\"payback\"], 2) for p in (0.12, 0.14, 0.20, 0.30)], \"| 5-yr break-even price:\", round(be, 3) if be else None)\n",
    "    print(\"\\n== Brine preheat share of working-fluid duty (sensible preheat, dT_pp = 10 K) ==\")\n",
    "    for n in SITES:\n",
    "        r = results[n]\n",
    "        for Tb in (100.0, 120.0, 140.0):\n",
    "            sh, Qb, mb = brine_share(n, Tb)\n",
    "            print(n, f\"brine {Tb:.0f} C: share {100*sh:.1f}%  Q_brine {Qb:.0f} kW  brine flow {mb:.0f} m3/h  flue-gas duty needed {r['Q'] - Qb:.0f} kW  (Q_in {r['Q']:.0f})\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    results, adv, tables, opt, rows8, rows9 = make_all()\n",
    "    report(results, adv, tables, opt)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}