"""Regenerate the SA-only CRM chart from the published CSV.

Run from the repository root:
    python3 validation/cases/data/dbns-crm-dpw5-sa/plot_grid_convergence.py
"""

import csv
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt


HERE = Path(__file__).resolve().parent
DATA = HERE / "cl05_results_sa.csv"
BAND_DATA = HERE / "dpw5_participant_band.csv"
OUTPUT = HERE.parent.parent / "images" / "dbns-crm-dpw5-sa" / "grid_convergence.png"

with DATA.open(newline="", encoding="utf-8") as source:
    rows = list(csv.DictReader(source))
with BAND_DATA.open(newline="", encoding="utf-8") as source:
    band_rows = list(csv.DictReader(source))

series = [
    ("hex", "WakeCFD SA Hex", "#1864ab", "o"),
    ("prism", "WakeCFD SA Prism", "#d9480f", "s"),
    ("hybrid", "WakeCFD SA Hybrid", "#2b8a3e", "^"),
]
panels = [
    ("alpha_cl05_deg", "alpha", "Angle of attack at CL = 0.5 [deg]"),
    ("cd_cl05_counts", "cd_counts", "Drag at CL = 0.5 [counts]"),
    ("cmy_cl05", "cmy", "Pitching moment at CL = 0.5"),
]

fig, axes = plt.subplots(3, 1, figsize=(8.1, 10.4), sharex=True, constrained_layout=True)
band_x = [float(row["grdfac_N^-2/3"]) for row in band_rows]
for ax, (field, band_field, ylabel) in zip(axes, panels):
    band_lo = [float(row[f"{band_field}_lo"]) if row[f"{band_field}_lo"] else float("nan") for row in band_rows]
    band_mid = [float(row[f"{band_field}_med"]) if row[f"{band_field}_med"] else float("nan") for row in band_rows]
    band_hi = [float(row[f"{band_field}_hi"]) if row[f"{band_field}_hi"] else float("nan") for row in band_rows]
    ax.fill_between(band_x, band_lo, band_hi, color="#d9d9d9", alpha=0.8,
                    label="DPW-5 participants (digitized range)", zorder=0)
    ax.plot(band_x, band_mid, color="#777777", linestyle=":", linewidth=1.6,
            label="DPW-5 median", zorder=1)
    for family, label, color, marker in series:
        points = sorted((row for row in rows if row["family"] == family),
                        key=lambda row: int(row["nodes"]), reverse=True)
        ax.plot(
            [int(row["nodes"]) ** (-2 / 3) for row in points],
            [float(row[field]) for row in points],
            label=label,
            color=color,
            marker=marker,
            linewidth=2,
            markersize=5.5,
            zorder=2,
        )
    ax.set_ylabel(ylabel)
    ax.set_xlim(0, 1.4e-4)
    ax.grid(alpha=0.28)

axes[-1].set_xlabel("N$^{-2/3}$ (N = mesh nodes; finer meshes to the left)")
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="upper center", bbox_to_anchor=(0.5, 1.08), ncol=3, frameon=False)
OUTPUT.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(OUTPUT, format="png", dpi=160, bbox_inches="tight")