# =====================================================================
# Companion code for "The Trap Number T: A Dimensionless Criterion for
# Boulder-Boulder Gravitational Dynamics on Near-Earth Rubble-Pile Asteroid
# Surfaces" (prepared for submission).
#
# THIS FILE: generates Figure 6 (stability of the trap point). Three panels:
#   (left)   illustration of a stable minimum (both surface-Hessian eigenvalues
#            > 0), the bowl z = x^2 + y^2 -- what a closed trap WOULD look like.
#   (middle) illustration of a saddle (eigenvalues of opposite sign),
#            z = x^2 - y^2 -- the Otohime case.
#   (right)  the REAL point-mass surface-Hessian eigenvalues vs grain distance:
#            radial lam_rad = -2Gm/r^3 (toward the boulder, NON-confining) and
#            contour lam_con = +Gm/r^3 (perpendicular, confining). They are of
#            opposite sign at every distance: a point-mass boulder potential is
#            trace-free (Laplace), so by Earnshaw's theorem a single boulder is
#            always a saddle and cannot form a minimum on its own.
#
# IMPORTANT: 𝒯 = 1 is shown as a NOMINAL order-unity scale, NOT a proven
# saddle-to-minimum bifurcation. The actual transition requires the curvature of
# the background geopotential and is set by that curvature versus the boulder
# term (manuscript Section 3.5; full treatment left to future work). The physical
# eigenvalues of the full Hessian at Otohime are -9.5e-5 and +1.0e-5 s^-2
# (a saddle), reproduced by code_hessian_stability.py.
# =====================================================================
import os
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import cm

TEAL  = "#16a085"   # lambda_1 (radial) line / "minimum" label
RED   = "#e8261e"   # lambda_2 (contour) line / "saddle" label
PINK_BG  = "#fce8e6"
GREEN_BG = "#e7f4ec"

T_OTO = 0.44        # Otohime trap number (long-axis catalogue value, Section 3.2)

def _style_3d(ax):
    """Clean, label-free 3D box matching the manuscript figure."""
    for pane in (ax.xaxis, ax.yaxis, ax.zaxis):
        pane.pane.set_facecolor("white")
        pane.pane.set_edgecolor("0.85")
        pane.pane.set_alpha(1.0)
    ax.grid(False)
    ax.set_xticks([]); ax.set_yticks([]); ax.set_zticks([])
    ax.view_init(elev=24, azim=-58)


def build_figure():
    fig = plt.figure(figsize=(15, 4.4))

    # ----- panel 1: stable minimum  z = x^2 + y^2 -----
    ax1 = fig.add_subplot(1, 3, 1, projection="3d")
    g = np.linspace(-2, 2, 60)
    X, Y = np.meshgrid(g, g)
    Zmin = X**2 + Y**2
    ax1.plot_surface(X, Y, Zmin, cmap="Blues", rstride=1, cstride=1,
                     linewidth=0, antialiased=True)
    ax1.set_title("Stable minimum\n(both eigenvalues > 0): a closed trap",
                  fontsize=11, pad=2)
    _style_3d(ax1)

    # ----- panel 2: saddle  z = x^2 - y^2 -----
    ax2 = fig.add_subplot(1, 3, 2, projection="3d")
    Zsad = X**2 - Y**2
    ax2.plot_surface(X, Y, Zsad, cmap="RdYlGn_r", rstride=1, cstride=1,
                     linewidth=0, antialiased=True)
    ax2.set_title("Saddle point\n(eigenvalues opposite sign): Otohime",
                  fontsize=11, pad=2)
    _style_3d(ax2)

    # ----- panel 3: REAL point-mass surface-Hessian eigenvalues vs grain distance -----
    ax3 = fig.add_subplot(1, 3, 3)
    Gc = 6.674e-11
    m_oto = (4.0 / 3.0) * np.pi * 1600.0 * (164.0 / 2.0)**3   # Otohime sphere mass
    r = np.linspace(10.0, 60.0, 300)                          # grain distance (m)
    lam_rad = -2.0 * Gc * m_oto / r**3                        # toward boulder: non-confining
    lam_con = +1.0 * Gc * m_oto / r**3                        # perpendicular: confining
    ax3.axhline(0, color="black", lw=1.0, zorder=1)
    ax3.plot(r, lam_rad, color=RED,  lw=2.6,
             label=r"$\lambda_{\rm radial}=-2Gm/r^3$ (non-confining)")
    ax3.plot(r, lam_con, color=TEAL, lw=2.6,
             label=r"$\lambda_{\rm contour}=+Gm/r^3$ (confining)")
    ax3.text(0.03, 0.05,
             "Opposite signs at every $r$:\npoint-mass field is trace-free\n"
             r"(Earnshaw) $\rightarrow$ always a saddle." + "\n"
             r"$\mathcal{T}=1$ nominal scale only.",
             transform=ax3.transAxes, ha="left", va="bottom", fontsize=8,
             bbox=dict(boxstyle="round", fc="white", ec="0.7", alpha=0.95))
    ax3.text(0.97, 0.97,
             "Full-Hessian at Otohime:\n"
             r"$(-9.5,\,+1.0)\times10^{-5}$ s$^{-2}$" + "\n= SADDLE",
             transform=ax3.transAxes, ha="right", va="top", fontsize=8,
             bbox=dict(boxstyle="round", fc="#fdf3e7", ec="0.7", alpha=0.95))
    ax3.set_xlim(10, 60); ax3.set_ylim(-5e-4, 2.5e-4)
    ax3.set_xlabel(r"grain distance from boulder centre  $r$ (m)", fontsize=10)
    ax3.set_ylabel(r"surface-Hessian eigenvalue (s$^{-2}$)", fontsize=10)
    ax3.set_title("A point-mass boulder is always a saddle", fontsize=10.5)
    ax3.legend(loc="lower right", fontsize=7.5, framealpha=0.9)

    fig.tight_layout()
    return fig


if __name__ == "__main__":
    here = os.path.dirname(os.path.abspath(__file__))
    out = os.path.join(here, "Fig6_stability_diagram.png")
    fig = build_figure()
    fig.savefig(out, dpi=200, bbox_inches="tight", pad_inches=0.1)
    print("saved:", out)
