#!/usr/bin/env python3
from pathlib import Path
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from PIL import Image
import numpy as np
import csv
import shutil

ROOT = Path("/home/yyq/OpenFOAM/yyq-v2206/run/meshRobustness_ESM2")
RAWO = ROOT / "raw_overlay_v2"
TABLE = ROOT / "tables" / "Table_S5.csv"
PDF = ROOT / "esm2" / "ESM_2.pdf"
SUBMIT = ROOT / "submission_ready"

cases = {
    "coarse":   {"label": "256×24×24", "color": "black",   "marker": "o"},
    "baseline": {"label": "320×30×30", "color": "#1f77b4", "marker": "o"},
    "fine":     {"label": "384×36×36", "color": "#d62728", "marker": "o"},
}

plt.rcParams["font.family"] = "sans-serif"
plt.rcParams["font.sans-serif"] = ["Arial", "Helvetica", "Liberation Sans", "DejaVu Sans"]
plt.rcParams["font.size"] = 8

def read_csv_numeric(path):
    path = Path(path)
    if not path.exists() or path.stat().st_size <= 1:
        return None, None
    with open(path, newline="") as f:
        lines = [ln for ln in f.read().splitlines() if ln.strip()]
    if not lines:
        return None, None
    rr = list(csv.reader(lines))
    if not rr or not rr[0]:
        return None, None
    header = rr[0]
    cols = {h: [] for h in header}
    for row in rr[1:]:
        if len(row) < len(header):
            row = row + [""] * (len(header) - len(row))
        for h, v in zip(header, row):
            try:
                cols[h].append(float(v))
            except:
                cols[h].append(np.nan)
    return header, {k: np.array(v, dtype=float) for k, v in cols.items()}

def pick_key(header, keys):
    fmap = {f: f.strip().lower() for f in header}
    for key in keys:
        kl = key.lower()
        for f, fl in fmap.items():
            if fl == kl or kl in fl:
                return f
    raise RuntimeError(f"Cannot find {keys} in {header}")

def get_points(path):
    header, data = read_csv_numeric(path)
    if header is None:
        return None
    xk = pick_key(header, ["Points:0", "Points_0", "x"])
    yk = pick_key(header, ["Points:1", "Points_1", "y"])
    x = data[xk] * 1000.0
    y = data[yk] * 1000.0
    if len(x) == 0 or np.all(np.isnan(x)):
        return None
    return x, y

# -------------------------------------------------
# Rebuild Fig. S15 as 2x2 with z = 0.013 added
# -------------------------------------------------
z_panels = [
    ("0.003", "z = 0.003 m"),
    ("0.010", "z = 0.010 m"),
    ("0.013", "z = 0.013 m"),
    ("0.017", "z = 0.017 m"),
]

fig, axs = plt.subplots(2, 2, figsize=(174/25.4, 5.8))
letters = ["(a)", "(b)", "(c)", "(d)"]

for ax, (zkey, ztitle), letter in zip(axs.ravel(), z_panels, letters):
    plotted = False
    for ckey, meta in cases.items():
        path = RAWO / ckey / "contour_span" / f"T1200_z{zkey}_t0.00062000.csv"
        pts = get_points(path)
        if pts is None:
            continue
        x, y = pts
        ax.scatter(x, y, s=3.0, color=meta["color"], marker=meta["marker"], linewidths=0)
        plotted = True

    ax.set_xlim(8, 28)
    ax.set_ylim(4, 16)
    ax.set_title(ztitle, fontsize=9, pad=3)
    ax.text(0.02, 0.96, letter, transform=ax.transAxes,
            ha="left", va="top", fontsize=8, fontweight="bold")

    ax.set_xlabel("x (mm)", fontsize=8)
    ax.set_ylabel("y (mm)", fontsize=8)
    ax.tick_params(axis="both", labelsize=7, length=2)

    if not plotted:
        ax.text(0.5, 0.5, "no resolved\nT = 1200 K\nsplit contour",
                transform=ax.transAxes, ha="center", va="center",
                fontsize=8, color="0.45")

    for sp in ax.spines.values():
        sp.set_linewidth(0.6)

handles = [
    Line2D([0], [0], marker="o", linestyle="None", markersize=4,
           markerfacecolor=cases["coarse"]["color"], markeredgewidth=0, label=cases["coarse"]["label"]),
    Line2D([0], [0], marker="o", linestyle="None", markersize=4,
           markerfacecolor=cases["baseline"]["color"], markeredgewidth=0, label=cases["baseline"]["label"]),
    Line2D([0], [0], marker="o", linestyle="None", markersize=4,
           markerfacecolor=cases["fine"]["color"], markeredgewidth=0, label=cases["fine"]["label"]),
]
fig.legend(handles=handles, loc="lower center", ncol=3, frameon=False, fontsize=8, bbox_to_anchor=(0.5, 0.01))

plt.subplots_adjust(left=0.10, right=0.98, bottom=0.13, top=0.92, wspace=0.22, hspace=0.35)
fig.savefig(ROOT / "plots" / "Fig_S15.tif", dpi=600, pil_kwargs={"compression": "tiff_lzw"})
fig.savefig(ROOT / "plots" / "Fig_S15.png", dpi=300)
plt.close(fig)

# -------------------------------------------------
# Repack ESM_2.pdf
# -------------------------------------------------
plt.rcParams["font.size"] = 10

article_title = "Three-Dimensional Flame-Head Bifurcation and Compression-Hydrodynamic Coupling in LES of Premixed Hydrogen-Air Propagation in a Half-Open Channel"
journal_name = "Flow, Turbulence and Combustion"
authors = "Youquan Yuan; Jiandong He*"
affiliation = "College of Mechanics and Aerospace, Inner Mongolia University of Technology, Hohhot, Inner Mongolia, China"
corr_email = "Corresponding author: hjd@imut.edu.cn"

captions = [
(
"Fig. S12. Center-plane T = 1200 K contour overlays in a local flame-head window for the coarse, baseline,\n"
"and fine meshes at t = 0.520, 0.580, and 0.620 ms. The same local window is used in all panels.\n"
"The late common-window bifurcation topology can be compared directly across the three meshes.",
ROOT / "plots" / "Fig_S12.tif"
),
(
"Fig. S13. Centerline overpressure profiles in the baseline phi = 1.0 case at t = 0.520, 0.580, and 0.620 ms.\n"
"The three meshes preserve the same qualitative buildup of the ahead-of-flame compression layer,\n"
"although the profile amplitude remains mesh-sensitive.",
ROOT / "plots" / "Fig_S13.tif"
),
(
"Fig. S14. Mesh-resolution comparison of the precursor metrics within the common stable interval (t <= 0.620 ms):\n"
"(a) Db(t), (b) dpAhead(t), and (c) trajectories in the (Db, dpAhead) plane.",
ROOT / "plots" / "Fig_S14.tif"
),
(
"Fig. S15. T = 1200 K contour overlays at t = 0.620 ms on four spanwise planes, z = 0.003, 0.010,\n"
"0.013, and 0.017 m, for the coarse, baseline, and fine meshes. The near-center off-plane section\n"
"(z = 0.013 m) helps show that the three-dimensional flame-head structure is not confined to a single center slice.",
ROOT / "plots" / "Fig_S15.tif"
),
]

def trim_bg(img, tol=6, pad=10):
    arr = np.asarray(img.convert("RGB"))
    bg = arr[0, 0, :].astype(int)
    diff = np.max(np.abs(arr.astype(int) - bg), axis=2)
    mask = diff > tol
    ys, xs = np.where(mask)
    if len(xs) == 0:
        return img
    x0 = max(0, int(xs.min()) - pad)
    x1 = min(arr.shape[1], int(xs.max()) + pad + 1)
    y0 = max(0, int(ys.min()) - pad)
    y1 = min(arr.shape[0], int(ys.max()) + pad + 1)
    return img.crop((x0, y0, x1, y1))

with PdfPages(PDF) as pdf:
    # Cover page
    fig = plt.figure(figsize=(8.27, 11.69))
    fig.text(0.08, 0.945, "Online Resource 2 (ESM_2.pdf)", fontsize=13)
    fig.text(0.08, 0.900, "Supplementary mesh-resolution robustness checks for", fontsize=11)
    fig.text(
        0.08, 0.852,
        "“Three-Dimensional Flame-Head Bifurcation and Compression-Hydrodynamic\n"
        "Coupling in LES of Premixed Hydrogen-Air Propagation in a Half-Open\n"
        "Channel”.",
        fontsize=11
    )
    fig.text(0.08, 0.790, journal_name, fontsize=10.5)
    fig.text(0.08, 0.755, authors, fontsize=10.5)
    fig.text(0.08, 0.720, affiliation, fontsize=10.5)
    fig.text(0.08, 0.685, corr_email, fontsize=10.5)
    fig.text(
        0.08, 0.610,
        "This supplementary file contains Fig. S12-S15 and Table S5. The added materials document the\n"
        "mesh-resolution study for the baseline phi = 1.0 case within the common stable interval t <= 0.620 ms.",
        fontsize=10.5, va="top"
    )
    bullets = [
        "Fig. S12: center-plane T = 1200 K contour overlays in a local flame-head window",
        "Fig. S13: centerline overpressure profiles at matched late common-window times",
        "Fig. S14: mesh-resolution comparison of Db and dpAhead trends",
        "Fig. S15: spanwise T = 1200 K contour overlays at t = 0.620 ms",
        "Table S5: tabulated mesh-comparison values at 0.520, 0.580, and 0.620 ms",
    ]
    y = 0.515
    for b in bullets:
        fig.text(0.10, y, u"\u2022 " + b, fontsize=10.5)
        y -= 0.045
    fig.text(0.50, 0.04, "1", ha="center", va="center", fontsize=9, color="0.35")
    plt.axis("off")
    pdf.savefig(fig, dpi=300)
    plt.close(fig)

    # Figure pages
    page_no = 2
    for cap, img_path in captions:
        fig = plt.figure(figsize=(8.27, 11.69))
        fig.text(0.06, 0.965, "Online Resource 2 · Supplementary figures and robustness checks",
                 fontsize=9.5, color="0.35")

        img = trim_bg(Image.open(img_path).convert("RGB"))
        ax_img = fig.add_axes([0.06, 0.43, 0.88, 0.28])
        ax_img.imshow(np.asarray(img))
        ax_img.axis("off")

        fig.text(0.06, 0.22, cap, fontsize=11.2, va="top", ha="left")
        fig.text(0.50, 0.04, str(page_no), ha="center", va="center", fontsize=9, color="0.35")
        pdf.savefig(fig, dpi=300)
        plt.close(fig)
        page_no += 1

    # Table page
    with open(TABLE, newline="") as f:
        rr = list(csv.reader(f))
    header = rr[0]
    body = rr[1:]

    wrapped_header = [
        "Mesh",
        "Db at\n0.520 ms\n(mm)",
        "Db at\n0.580 ms\n(mm)",
        "Db at\n0.620 ms\n(mm)",
        "dpAhead at\n0.520 ms\n(kPa)",
        "dpAhead at\n0.580 ms\n(kPa)",
        "dpAhead at\n0.620 ms\n(kPa)",
        "t(dpAhead >=\n3 kPa)\n(ms)",
    ]

    fig = plt.figure(figsize=(11.69, 8.27))
    fig.text(0.05, 0.95, "Online Resource 2 · Supplementary figures and robustness checks",
             fontsize=9.5, color="0.35")
    fig.text(
        0.05, 0.90,
        "Table S5. Mesh-resolution comparison for the baseline phi = 1.0 case within the common stable interval.",
        fontsize=11.2
    )

    ax = fig.add_axes([0.04, 0.30, 0.92, 0.32])
    ax.axis("off")
    tbl = ax.table(cellText=body, colLabels=wrapped_header, cellLoc="center", loc="center")
    tbl.auto_set_font_size(False)
    tbl.set_fontsize(8.9)

    widths = [0.14, 0.12, 0.12, 0.12, 0.14, 0.14, 0.14, 0.16]
    for (r, c), cell in tbl.get_celld().items():
        if c < len(widths):
            cell.set_width(widths[c])
        cell.set_edgecolor("0.35")
        cell.set_linewidth(0.5)
        if r == 0:
            cell.set_height(0.24)
            cell.set_facecolor("#f2f2f2")
            cell.set_text_props(weight="bold")
            cell.get_text().set_wrap(True)
        else:
            cell.set_height(0.18)
            if r % 2 == 1:
                cell.set_facecolor("#fbfbfb")

    footer = (
        "Values are reported at matched late common-window times (0.520, 0.580, and 0.620 ms).\n"
        "All metric columns are sampled at the same matched times for coarse, baseline, and fine meshes.\n"
        "Absolute Db values remain mesh-sensitive at matched times; the table is used to support topology/trend robustness rather than strict quantitative convergence."
    )
    fig.text(0.05, 0.12, footer, fontsize=9.2, va="top")
    fig.text(0.50, 0.03, str(page_no), ha="center", va="center", fontsize=9, color="0.35")
    pdf.savefig(fig, dpi=300)
    plt.close(fig)

SUBMIT.mkdir(exist_ok=True)
shutil.copy2(ROOT / "plots" / "Fig_S15.tif", SUBMIT / "Fig_S15.tif")
shutil.copy2(TABLE, SUBMIT / "Table_S5.csv")
shutil.copy2(PDF, SUBMIT / "ESM_2.pdf")

print("rebuilt:", ROOT / "plots" / "Fig_S15.tif")
print("rebuilt:", PDF)
