"""Original phase portrait. Python 3.14; NumPy and Matplotlib from pyproject.toml.

Run from the intended PNG output directory. No source-paper image is used.
"""

from pathlib import Path

import matplotlib

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


def rhs(z):
    x, y = z
    return np.array([x * (1 - y) * (1 - x - y), y * (2 - x * y)])


def trajectory(start, dt=0.002, steps=10000):
    z = np.array(start, dtype=float)
    values = [z.copy()]
    for _ in range(steps):
        k1 = rhs(z)
        k2 = rhs(z + dt * k1 / 2)
        k3 = rhs(z + dt * k2 / 2)
        k4 = rhs(z + dt * k3)
        z = z + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
        if not np.all(np.isfinite(z)) or np.max(np.abs(z)) > 50:
            raise RuntimeError("Trajectory left the numerical drawing range")
        values.append(z.copy())
    return np.array(values)


def main():
    fig, ax = plt.subplots(figsize=(9.0, 6.2), facecolor="white")
    ax.set_facecolor("white")
    xx, yy = np.meshgrid(np.arange(0.18, 5.45, 0.38),
                         np.arange(0.16, 3.4, 0.29))
    vx = xx * (1 - yy) * (1 - xx - yy)
    vy = yy * (2 - xx * yy)
    speed = np.hypot(vx, vy)
    ax.quiver(xx, yy, vx / speed, vy / speed, angles="xy", scale_units="xy",
              scale=8, color="#999999", alpha=0.65, width=0.0023, zorder=1)
    blue, orange = "#2567a5", "#cc7816"
    x = np.linspace(0.04, 5.5, 1000)
    ax.axvline(0, color=blue, lw=2.1)
    ax.axhline(1, color=blue, lw=1.7, label=r"$\dot{x}=0$")
    ax.plot([0, 1], [1, 0], color=blue, lw=1.7)
    ax.axhline(0, color=orange, lw=2.1)
    ax.plot(x, 2 / x, color=orange, lw=1.7, label=r"$\dot{y}=0$")
    starts = [(0.22, 0.12), (0.85, 0.04), (4.7, 0.15), (4.8, 2.8),
              (0.2, 2.0), (1.3, 2.8), (3.1, 0.8)]
    summaries = []
    for i, start in enumerate(starts):
        points = trajectory(start)
        ax.plot(points[:, 0], points[:, 1], color="#26724b", lw=1.2,
                alpha=0.9, label="Forward trajectories" if i == 0 else None)
        lengths = np.concatenate([[0], np.cumsum(np.linalg.norm(np.diff(points, axis=0), axis=1))])
        for fraction in [0.23, 0.63]:
            j = min(np.searchsorted(lengths, fraction * lengths[-1]), len(points) - 15)
            if 0 < points[j, 0] < 5.45 and 0 < points[j, 1] < 3.35:
                direction = rhs(points[j])
                direction = direction / np.linalg.norm(direction)
                end = points[j] + 0.095 * direction
                ax.annotate("", xy=end, xytext=points[j],
                            arrowprops={"arrowstyle": "->", "color": "#26724b", "lw": 1.4})
        summaries.append((start, points[-1].tolist()))
    ax.scatter([0, 1, 2], [0, 0, 1], s=[48, 48, 55], color="black", zorder=5,
                clip_on=False)
    ax.annotate("Unstable node", (0, 0), xytext=(0.12, 0.47),
                arrowprops={"arrowstyle": "-", "color": "black"}, fontsize=10)
    ax.annotate("Saddle", (1, 0), xytext=(1.15, 0.3),
                arrowprops={"arrowstyle": "-", "color": "black"}, fontsize=10)
    ax.annotate("Clockwise spiral sink", (2, 1), xytext=(2.65, 1.35),
                arrowprops={"arrowstyle": "-", "color": "black"}, fontsize=10)
    ax.text(0.25, 0.84, r"$x+y=1$", color=blue, fontsize=10)
    ax.text(4.75, 1.07, r"$y=1$", color=blue, fontsize=10)
    ax.text(0.75, 2.85, r"$xy=2$", color=orange, fontsize=10)
    for start, end in [((0.05, 0.6), (0.05, 0.95)),
                       ((0.25, 0.03), (0.55, 0.03)),
                       ((1.7, 0.03), (1.4, 0.03))]:
        ax.annotate("", xy=end, xytext=start,
                    arrowprops={"arrowstyle": "->", "color": "black", "lw": 1.3})
    ax.set(xlim=(-0.02, 5.5), ylim=(-0.02, 3.4), xlabel="$x$", ylabel="$y$")
    ax.set_title("Polynomial nullclines and the attracting focus", fontsize=14)
    ax.legend(loc="upper right", framealpha=0.95, fontsize=10)
    ax.spines[["top", "right"]].set_visible(False)
    fig.tight_layout()
    fig.savefig(Path.cwd() / "2015-polynomial-phase-plane.png", dpi=140,
                facecolor="white")
    plt.close(fig)
    print(summaries)


if __name__ == "__main__":
    main()