fix: hysteresis shading, label collision, axis limits; close reference loop

Author: Hj1308

README.md

@@ -27,7 +27,7 @@ Developed and validated for **graphene-like carbon nitride (C₃N₄)**, MOFs, z
 
 *Output of `bet_analysis.py` on the bundled reference dataset — a synthetic
 Type IV isotherm whose true surface area is known by construction
-(S_BET = n_m × 4.353 = 177.00 m² g⁻¹; the tool recovers 177.01, R² = 0.999999).
+(S_BET = n_m × 4.353 = 177.28 m² g⁻¹; the tool recovers 177.28, R² = 0.999999).
 Reproduce it with:*
 
 ```bash
@@ -202,6 +202,10 @@ Cychosz & Thommes 2018 §3). §6.1 also recommends argon at 87 K over nitrogen
 at 77 K where surface functional groups interact with the N2 quadrupole.
 ```
 
+*This example message is from a measurement whose lowest adsorption point sits
+above p/p₀ = 0.015. The bundled reference dataset is purely mesoporous and does
+not trigger this gate.*
+
 | Check | Behaviour |
 |-------|-----------|
 | **BET C constant** | `UserWarning` raised if C < 0 — invalid p/p₀ range; adjust `start_pt`/`end_pt` to 0.05 ≤ p/p₀ ≤ 0.35 |
@@ -317,6 +321,9 @@ It contains no measured data. All four sheets (AdsDes, BET, BJH, Summary)
 derive from the same monolayer capacity and pore-size distribution, so the file
 is internally self-consistent: S_BET and S_BJH agree to within 1.2 %.
 
+Every figure in this README was produced from this synthetic file. No measured
+instrument data is included in this repository.
+
 ---
 
 ## 📚 References

assets/bet_analysis_example.png

@@ -775,10 +775,19 @@ def plot_all(data: dict, iso_cls: dict, hyst_cls: dict,
         sort_d = np.argsort(des[:, 0])[::-1]
         ax.plot(des[sort_d, 0], des[sort_d, 1], "s--",
                 color=C_DES, ms=4, lw=1.4, label="Desorption")
-        pp0_d_s, Va_d_s = des[sort_d, 0], des[sort_d, 1]
-        pp0_all = np.concatenate([ads[:, 0], pp0_d_s[::-1]])
-        Va_all  = np.concatenate([ads[:, 1], Va_d_s[::-1]])
-        ax.fill(pp0_all, Va_all, alpha=0.10, color=C_ADS)
+        # Shade only the hysteresis loop: the p/p0 interval where both branches
+        # exist, bounded by the two interpolated branches (not a single polygon
+        # with straight closing edges, which produced a full-width wedge).
+        p_lo = max(ads[:, 0].min(), des[:, 0].min())
+        p_hi = min(ads[:, 0].max(), des[:, 0].max())
+        if p_hi > p_lo:
+            p_grid = np.linspace(p_lo, p_hi, 200)
+            s_a = np.argsort(ads[:, 0])
+            s_d = np.argsort(des[:, 0])
+            Va_a_g = np.interp(p_grid, ads[s_a, 0], ads[s_a, 1])
+            Va_d_g = np.interp(p_grid, des[s_d, 0], des[s_d, 1])
+            ax.fill_between(p_grid, Va_a_g, Va_d_g, alpha=0.10,
+                            color=C_ADS, linewidth=0)
 
     ax.set_xlabel(r"Relative Pressure ($p/p_0$)")
     ax.set_ylabel(r"Volume Adsorbed (cm$^3$ g$^{-1}$ STP)")
@@ -791,7 +800,7 @@ def plot_all(data: dict, iso_cls: dict, hyst_cls: dict,
     iso_label  = iso_cls["type"]
     hyst_label = hyst_cls["type"] if hyst_cls["type"] != "None" else ""
     ax_ann = iso_label + (f" / {hyst_label}" if hyst_label else "")
-    ax.text(0.03, 0.96, ax_ann, transform=ax.transAxes,
+    ax.text(0.03, 0.78, ax_ann, transform=ax.transAxes,
             va="top", ha="left", fontsize=8.5,
             bbox=dict(boxstyle="round,pad=0.3", fc="white",
                       ec="0.7", lw=0.7, alpha=0.9))
@@ -833,6 +842,16 @@ def plot_all(data: dict, iso_cls: dict, hyst_cls: dict,
     rp    = bjh[:, 0] * 2          # radius (nm) -> diameter (nm)
     dVdd  = bjh[:, 1] / 2.0        # dVp/drp -> dVp/ddp
 
+    # Upper x-limit from the data: the smallest diameter where the cumulative
+    # pore volume reaches 99 % of its final value (drops the mostly-empty
+    # high-diameter tail), with a floor so microporous samples aren't cramped.
+    x_max = float(rp[-1])
+    if bjh[-1, 2] > 0:
+        done = np.where(bjh[:, 2] >= 0.99 * bjh[-1, 2])[0]
+        if len(done):
+            x_max = float(rp[done[0]])
+    x_max = max(x_max, 5.0)
+
     ax.plot(rp, dVdd, "-", color=C_BJH, lw=1.5)
     ax.fill_between(rp, dVdd, alpha=0.15, color=C_BJH)
 
@@ -843,15 +862,15 @@ def plot_all(data: dict, iso_cls: dict, hyst_cls: dict,
 
     ax.set_xlabel(r"Pore Diameter (nm)")
     ax.set_ylabel(r"d$V_p$/d$d_p$  (cm$^3$ g$^{-1}$ nm$^{-1}$)")
-    ax.set_xlim(left=0)
+    ax.set_xlim(left=0, right=x_max)
     ax.set_ylim(bottom=0)
     ax.xaxis.set_minor_locator(AutoMinorLocator())
     ax.yaxis.set_minor_locator(AutoMinorLocator())
 
     ax.axvline(N2_CAVITATION_NM, ls=":", lw=0.8, color="0.6", alpha=0.7)
-    ax.text(N2_CAVITATION_NM + 0.2,
-            ax.get_ylim()[1] * 0.01 if ax.get_ylim()[1] > 0 else 0.001,
-            "cavitation\n(~3.4 nm)", fontsize=6.5, color="0.5", va="bottom")
+    ax.text(N2_CAVITATION_NM, ax.get_ylim()[1] * 0.9,
+            "cavitation\n(~3.4 nm)", fontsize=6.5, color="0.5",
+            va="top", ha="left")
     _label_panel(ax, "C")
 
     # ── [D] Cumulative Pore Volume ────────────────────────────
@@ -876,7 +895,7 @@ def plot_all(data: dict, iso_cls: dict, hyst_cls: dict,
     ax2.set_ylabel(r"Cum. Surface Area (m$^2$ g$^{-1}$)", color=C_BJH)
     ax.tick_params(axis="y", colors=C_CUM)
     ax2.tick_params(axis="y", colors=C_BJH)
-    ax.set_xlim(left=0)
+    ax.set_xlim(left=0, right=x_max)
     ax.set_ylim(bottom=0)
 
     lines1, lbl1 = ax.get_legend_handles_labels()

bet_analysis.py

@@ -72,7 +72,7 @@
 # Reuse the project's own constants and physics generators (no second BET
 # implementation, no duplicated physical constants).
 from bet_analysis import N2_BET_FACTOR, N2_STP_TO_LIQUID  # noqa: E402
-from synthetic_isotherms import bet, step, desorption     # noqa: E402
+from synthetic_isotherms import bet, step, desorption_closed  # noqa: E402
 
 _trapezoid = getattr(np, "trapezoid", None) or np.trapz
 
@@ -94,10 +94,12 @@
 STEP_W = 0.08
 STEP_H = 200.0                   # cm³(STP) g⁻¹ added by capillary condensation
 
-# Desorption branch (H1 loop): adsorption pressure shifted down by SHIFT above
-# CLOSE_AT (mirrors the TypeIV_H1 fixture in synthetic_isotherms).
-SHIFT    = 0.15
-CLOSE_AT = 0.45
+# Desorption branch (H1 loop, closed at both ends): a volume-domain bump
+# vd = v_ads + AMPLITUDE·sin(pi·u) that vanishes at CLOSE_AT and p/p0 = 1
+# (desorption_closed). AMPLITUDE = 25 keeps the closure gap < 0.02, the branch
+# monotonic, slope_ratio_max <= 1.10, and a visible H1 loop (area_norm ~ 0.024).
+AMPLITUDE = 25.0
+CLOSE_AT  = 0.45
 
 # BJH mesopore peak (radius, nm) → 12 nm diameter peak.
 R_PEAK  = 6.0
@@ -120,7 +122,7 @@ def make_isotherm(n_m):
     x = _grid()
     va = bet(x, n_m, C, n=N_LAYERS) + step(x, X_STEP, STEP_W, STEP_H)
     ads = np.column_stack([x, va])
-    des = desorption(x, va, SHIFT, CLOSE_AT)
+    des = desorption_closed(x, va, AMPLITUDE, CLOSE_AT)
     return x, ads, des
 
 
@@ -132,7 +134,7 @@ def make_bjh(vp_total):
     trapezoidal integrals of ``dV`` over the same radius grid (cylindrical
     pore geometry: dS = 2 dV / r, with 1 cm³→1e-6 m³ and 1 nm→1e-9 m).
     """
-    rp = np.logspace(np.log10(1.0), np.log10(100.0), 60)
+    rp = np.logspace(np.log10(1.0), np.log10(25.0), 60)
     dV = np.exp(-0.5 * ((np.log(rp) - np.log(R_PEAK)) / R_SIGMA) ** 2)
     dV /= _trapezoid(dV, rp)
     dV *= vp_total

examples/make_reference_data.py

@@ -88,6 +88,22 @@ def desorption(x, v_ads, shift, close_at):
     vd = np.maximum(vd, np.interp(xd, x, v_ads))
     return np.column_stack([xd, vd])
 
+
+def desorption_closed(x, v_ads, amplitude, close_at):
+    """H1 desorption branch that closes at close_at.
+
+    Separation is added in the *volume* domain (vd = v_ads + A·sin(pi·u))
+    rather than the *pressure* domain, so the two branches stay nearly
+    parallel: condensation and evaporation keep the same steepness, as
+    expected for a single pore-size distribution.  The bump vanishes at both
+    ends, so the loop is closed at the bottom and the top.
+    """
+    m = x >= close_at
+    xd = x[m]
+    u = (xd - close_at) / (1.0 - close_at)
+    vd = v_ads[m] + amplitude * np.sin(np.pi * u)
+    return np.column_stack([xd, vd])
+
 CASES = {
     "TypeIa_noHyst":  (type_Ia,  None),
     "TypeIb_noHyst":  (type_Ib,  None),