"""
.. _pah_phase_diagram_example:

PAH size-dependent melting vs. mechanism solid-carbon thresholds.
------------------------------------------------------------------

Every Cantera mechanism Bloc supports resolves its own solid-carbon
carbon-count threshold, ``n_C_min``, via
:func:`~bloc.chem.solid_carbon_threshold.resolve_n_c_min` (see
:mod:`bloc.chem.solid_carbon_threshold` for why the right value depends on
the mechanism). This example places each *registered* mechanism's threshold
on an (approximate) PAH melting-point-vs-size phase diagram, so it's visually
clear which mechanisms treat a given PAH size as "solid" well past the
liquid/solid transition -- and which are more aggressive.

The phase diagram (``Liquid`` / ``Size-dependent`` / ``Solid`` regions, five
labeled bulk-PAH melting points, and the two dashed boundary lines
separating the regions) is digitized **by eye, approximately** from Fig. 5 of
Chen et al. (2015), "Solid-liquid transitions in homogenous ovalene,
hexabenzocoronene and circumcoronene clusters: A molecular dynamics study",
Combustion and Flame, https://doi.org/10.1016/j.combustflame.2014.07.025
(see ``data/literature/README.md`` for provenance and caveats). This is a
qualitative, order-of-magnitude comparison, not a quantitative one --
refer to the original publication for exact values.

That figure compares melting points across three cluster sizes of the same
PAH molecule:

* **PAH50** -- A cluster consisting of 50 molecules of a given peri-condensed
  PAH. Its melting point is treated as a lower bound for the size-dependent
  region, and a PAH50 cluster is about 4 nm in diameter, making it an
  analogue for nascent soot particles.
* **PAH200** -- A cluster consisting of 200 molecules of a given PAH. Its
  size is approximately 1.5 times that of PAH50. It's used to show how the
  cluster melting point changes with cluster size.
* **PAHbulk** -- The bulk (macroscopic) solid of the same PAH. Its melting
  point is treated as an upper bound for the size-dependent region, and it
  serves as an approximation for mature soot particles.
"""  # noqa: D205, D400

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from bloc.chem.solid_carbon_threshold import registered_mechanisms
from bloc.paths import ROOT_FOLDER_PATH

# %%
# Load the digitized bulk-PAH melting-point data and phase-boundary points
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

literature_dir = ROOT_FOLDER_PATH / "data" / "literature"
pah_data = pd.read_csv(
    literature_dir / "chen2015_pah_bulk_melting_point.csv", comment="#"
)
boundary_points = pd.read_csv(
    literature_dir / "chen2015_phase_boundary_points.csv", comment="#"
)
pah_data

# %%
# Map a carbon count to an approximate ring count / mass
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#
# A linear fit through the 5 digitized (n_C, rings) / (n_C, mass) points --
# good enough to *place* a mechanism's n_C_min on the same axes as the
# figure, not a physical model of PAH structure.

ring_fit = np.polyfit(pah_data["n_C"], pah_data["n_aromatic_rings"], 1)
mass_fit = np.polyfit(pah_data["n_C"], pah_data["MW_amu"], 1)


def n_c_to_mass(n_c: float) -> float:
    return float(np.polyval(mass_fit, n_c))


# %%
# Group registered mechanisms by their resolved n_C_min
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#
# Several mechanisms share the same threshold (e.g. the whole CRECK
# sectional-soot family uses n_C_min=80) -- group them so the plot shows one
# line per distinct threshold instead of a dozen overlapping ones.
# Mechanisms with n_c_min=None (no PAH continuum, e.g. Fincke_GRC: only
# named C(s)/C(soot) ever count as solid) have no carbon-count threshold to
# place on this axis at all -- listed separately below instead.

by_n_c_min: dict[int, list[str]] = {}
species_names_only: list[str] = []
for basename, result in sorted(registered_mechanisms().items()):
    if result.n_c_min is None:
        species_names_only.append(basename)
        continue
    by_n_c_min.setdefault(result.n_c_min, []).append(basename)

print("Mechanisms with no carbon-count threshold (named species only):")
for basename in species_names_only:
    result = registered_mechanisms()[basename]
    print(f"  {basename}: solid species = {result.solid_species_names}")

# %%
# Plot
# ^^^^

fig, ax = plt.subplots(figsize=(9, 6.5))

# Phase-diagram boundaries (Liquid <-> Size-dependent, Size-dependent <->
# Solid): each fit through its two digitized anchor points (A+C, B+D --
# see data/literature/README.md). A straight-line fit, not a physical model.
mass_range = np.linspace(0, 1500, 200)


def _boundary_fit(name: str):
    pts = boundary_points[boundary_points["boundary"] == name]
    slope, intercept = np.polyfit(pts["MW_amu"], pts["melting_point_K"], 1)
    return slope * mass_range + intercept


upper_boundary = _boundary_fit("liquid_size_dependent")  # Size-dependent <-> Liquid
lower_boundary = _boundary_fit("size_dependent_solid")  # Solid <-> Size-dependent

ax.fill_between(
    mass_range, 0, lower_boundary, color="tab:blue", alpha=0.20, label="Solid"
)
ax.fill_between(
    mass_range,
    lower_boundary,
    upper_boundary,
    color="tab:red",
    alpha=0.18,
    label="Size-dependent",
)
ax.fill_between(
    mass_range, upper_boundary, 2700, color="tab:green", alpha=0.18, label="Liquid"
)

# "Flame temperature" band, as in the source figure: a shaded horizontal
# strip between the two reference temperatures the boundary points (A/B and
# C/D) were digitized at (2000 K and 1500 K), with solid horizontal lines at
# each edge and a label.
ax.axhspan(1500, 2000, color="gray", alpha=0.15, zorder=0)
ax.axhline(2000, color="k", linewidth=0.8)
ax.axhline(1500, color="k", linewidth=0.8)
ax.text(
    30, 1750, "Flame temperature", fontsize=10, va="center", ha="left", color="dimgray"
)

# Region labels, as in the source figure (italic text inside each region,
# rather than relying on the legend alone). "Size Dependent" is labeled twice
# (center and upper-right) since that band spans the full width of the plot.
ax.text(250, 2350, "Liquid", fontsize=16, fontstyle="italic", ha="left")
ax.text(650, 1300, "Size\nDependent", fontsize=13, fontstyle="italic", ha="center")
ax.text(1300, 2600, "Size Dependent", fontsize=11, fontstyle="italic", ha="right")
ax.text(550, 100, "Solid", fontsize=16, fontstyle="italic", ha="left")

# Three cluster sizes of the same PAH, matching the source figure's own
# marker convention: squares (PAHbulk), crosses (PAH200), triangles (PAH50).
# Markers only (no connecting line) -- the digitized points aren't dense
# enough for a line between them to mean anything beyond the marker itself.
ax.plot(
    pah_data["MW_amu"],
    pah_data["melting_point_K"],
    "ks",
    label=r"PAH$_{bulk}$ melting point (digitized)",
)
ax.plot(
    pah_data["MW_amu"],
    pah_data["melting_point_K_pah200"],
    "kx",
    label=r"PAH$_{200}$ melting point (digitized)",
)
ax.plot(
    pah_data["MW_amu"],
    pah_data["melting_point_K_pah50"],
    "k^",
    label=r"PAH$_{50}$ melting point (digitized)",
)
for _, row in pah_data.iterrows():
    ax.annotate(
        row["formula"],
        (row["MW_amu"], row["melting_point_K"]),
        textcoords="offset points",
        xytext=(6, -4),
        fontsize=11,
    )

# The 4 digitized boundary-line anchor points themselves (A, B, C, D in the
# source figure), so the fitted boundaries above can be visually checked
# against the points they were fit through.
ax.scatter(
    boundary_points["MW_amu"],
    boundary_points["melting_point_K"],
    marker="D",
    facecolors="none",
    edgecolors="k",
    s=40,
    label="Digitized phase-boundary points (A, B, C, D)",
)
for _, row in boundary_points.iterrows():
    ax.annotate(
        row["label"],
        (row["MW_amu"], row["melting_point_K"]),
        textcoords="offset points",
        xytext=(6, 4),
        fontsize=8,
    )

for n_c_min, basenames in sorted(by_n_c_min.items()):
    mass = n_c_to_mass(n_c_min)
    ax.axvline(mass, color="k", linestyle="--", linewidth=0.8, alpha=0.6)
    label = f"n_C_min={n_c_min}\n" + "\n".join(basenames)
    ax.text(mass, 2550, label, rotation=90, va="top", ha="right", fontsize=6.5)

ax.set_xlim(0, 1500)
ax.set_ylim(0, 2700)
ax.set_xlabel(
    "PAH mass, amu (approximate; converted from n_C_min via the linear fit above)"
)
ax.set_ylabel("Melting point, K")
ax.set_title("PAH bulk melting point vs. size, with each mechanism's n_C_min overlaid")
ax.legend(loc="lower right", fontsize=8)
plt.tight_layout()
plt.show()