"""Original symmetric lattice-model illustration; no measured polymer data.

Run from the mirrored media directory with Python 3.14, NumPy 2.3.5 and
Matplotlib 3.10.7. Both chains have 100 equal-volume segments.
"""

from pathlib import Path

import matplotlib

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


N_SEGMENTS = 100
OUTPUT = "polymer-blend-miscibility-landscape.png"


def free_energy(phi, chi):
    return (phi * np.log(phi) + (1 - phi) * np.log1p(-phi)) / N_SEGMENTS + chi * phi * (1 - phi)


def derivative(phi, chi):
    return np.log(phi / (1 - phi)) / N_SEGMENTS + chi * (1 - 2 * phi)


def curvature(phi, chi):
    return 1 / (N_SEGMENTS * phi) + 1 / (N_SEGMENTS * (1 - phi)) - 2 * chi


def coexistence_left(chi):
    low, high = 1e-8, 0.25
    assert derivative(low, chi) < 0 < derivative(high, chi)
    for _ in range(80):
        mid = (low + high) / 2
        if derivative(mid, chi) < 0:
            low = mid
        else:
            high = mid
    return (low + high) / 2


def main():
    phi = np.linspace(0.0001, 0.9999, 1600)
    coexist = coexistence_left(0.03)
    coexist_other = 1 - coexist
    common_height = free_energy(coexist, 0.03)
    np.testing.assert_allclose(derivative(coexist, 0.03), 0, atol=1e-12)
    np.testing.assert_allclose(free_energy(coexist_other, 0.03), common_height, atol=1e-12)
    assert curvature(0.5, 0.01) > 0 > curvature(0.5, 0.03)
    assert np.min(free_energy(phi, 0.03)) >= common_height - 1e-8

    plt.rcParams.update({"font.size": 10, "axes.spines.top": False, "axes.spines.right": False})
    fig, ax = plt.subplots(figsize=(7.6, 4.8), dpi=100, facecolor="white")
    ax.plot(phi, 1000 * free_energy(phi, 0.01), color="#007c91", lw=2.2, label=r"$\chi=0.01$: one minimum")
    ax.plot(phi, 1000 * free_energy(phi, 0.03), color="#ba4d00", lw=2.2, label=r"$\chi=0.03$: two minima")
    ax.plot([coexist, coexist_other], [1000 * common_height] * 2, color="#555555", ls="--", lw=1.5, label="Two-phase common tangent")
    ax.scatter([coexist, coexist_other], [1000 * common_height] * 2, color="#ba4d00", zorder=4)
    ax.annotate(f"Coexisting compositions\n{coexist:.4f} and {coexist_other:.4f}", xy=(coexist_other, 1000 * common_height), xytext=(0.62, -2.6), arrowprops={"arrowstyle": "->", "color": "#555555"}, color="#555555", ha="center", fontsize=9)
    ax.set(xlim=(0, 1), ylim=(-5.1, 1.1), xlabel=r"First-polymer volume fraction $\phi$", ylabel=r"$10^3 f(\phi)$ (dimensionless)")
    ax.grid(alpha=0.18)
    ax.legend(loc="upper center", bbox_to_anchor=(0.5, 1.12), ncol=3, frameon=False, fontsize=8)
    ax.set_title("Polymer blend mixing: one minimum or phase coexistence", fontsize=12, pad=43)
    fig.text(0.5, 0.035, "Original equal-segment lattice model, N₁ = N₂ = 100. No specific polymer pair or data fitted.", ha="center", fontsize=9)
    fig.subplots_adjust(left=0.105, right=0.975, top=0.77, bottom=0.18)
    fig.savefig(Path.cwd() / OUTPUT, facecolor="white", transparent=False)
    plt.close(fig)


if __name__ == "__main__":
    main()