#!/usr/bin/env python3
"""
================================================================================
hessian_stability.py
Stability analysis of the trap point via the surface-tangent Hessian
of the effective potential.

A critical point of the effective potential is a stable trap only if both
eigenvalues of the surface Hessian are positive. This code evaluates the
Hessian numerically at the Otohime critical point on Ryugu.

Result: the eigenvalues have opposite signs (a saddle point) at the real Otohime
value 𝒯 ~= 0.5. This is the synthetic-model computation; the headline Ryugu result
uses the real catalogue (see code_real_ryugu_trap_number.py).
Note: because each boulder is modelled as a point mass, its potential satisfies
Laplace's equation away from the source, so its Hessian is trace-free and a single
boulder cannot by itself produce a potential minimum (Earnshaw's theorem). A stable
trap therefore requires the curvature of the background geopotential; the transition
is governed by that background curvature versus the boulder term, not pinned exactly
at 𝒯 = 1. 𝒯 = 1 is reported as the natural order-unity scale, not a rigorously
derived bifurcation point (manuscript Section 3.5)."

Companion code for "The Trap Number 𝒯" (prepared for submission).

Usage:
    python hessian_stability.py
================================================================================
"""

import numpy as np
from code_trap_number import ASTEROIDS, generate_boulders, latlon_to_xyz

G = 6.674e-11


def effective_potential(lat, lon, params, bpos, masses, radii, omega):
    """
    Effective potential per unit mass at surface point (lat, lon):
    bulk gravity + centrifugal + N-body boulder gravity.
    The minimum distance to a boulder is its radius (a test particle cannot
    be inside a boulder).
    """
    R, M = params['R'], params['M']
    P = latlon_to_xyz(lat, lon, R)
    x, y, z = P
    phi = -G * M / R                          # bulk gravity (point-mass approx)
    phi += -0.5 * omega**2 * (x**2 + y**2)    # centrifugal
    dist = np.sqrt(np.sum((bpos - P)**2, axis=1))
    dist_eff = np.maximum(dist, radii)        # cannot be inside a boulder
    nearby = dist < 250                       # only nearby boulders matter
    phi += -np.sum(G * masses[nearby] / dist_eff[nearby])
    return phi


def surface_hessian(lat0, lon0, params, bpos, masses, radii, omega, h=0.1):
    """
    Numerical 2x2 surface-tangent Hessian (in physical metres) at (lat0, lon0).
    Returns the two eigenvalues (units 1/s^2).
    """
    R = params['R']
    m_per_deg = np.pi * R / 180.0
    cos_lat = np.cos(np.radians(lat0))

    def phi(i, j):
        return effective_potential(lat0 + i * h, lon0 + j * h,
                                   params, bpos, masses, radii, omega)

    d2_lat = (phi(1, 0) - 2 * phi(0, 0) + phi(-1, 0)) / (h * m_per_deg)**2
    d2_lon = (phi(0, 1) - 2 * phi(0, 0) + phi(0, -1)) / (h * m_per_deg * cos_lat)**2
    d2_cross = (phi(1, 1) - phi(1, -1) - phi(-1, 1) + phi(-1, -1)) / \
               (4 * h * h * m_per_deg**2 * cos_lat)
    H = np.array([[d2_lat, d2_cross], [d2_cross, d2_lon]])
    return np.linalg.eigvalsh(H)


if __name__ == '__main__':
    params = ASTEROIDS['Ryugu']
    omega = 2 * np.pi / (7.6326 * 3600)  # Ryugu spin rate (rad/s)

    d, masses, blat, blon = generate_boulders(params, seed=0)
    radii = d / 2
    bpos = np.array([latlon_to_xyz(blat[i], blon[i], params['R'])
                     for i in range(len(blat))])

    # Otohime is at lat = -80, lon = 180
    eigs = surface_hessian(-80.0, 180.0, params, bpos, masses, radii, omega)

    print("Hessian stability analysis at Otohime (Ryugu)")
    print(f"  Eigenvalue 1: {eigs[0]:+.3e} s^-2")
    print(f"  Eigenvalue 2: {eigs[1]:+.3e} s^-2")
    if (eigs > 0).all():
        print("  -> STABLE MINIMUM (true 2-D trap)")
    elif eigs[0] * eigs[1] < 0:
        print("  -> SADDLE POINT (confining in one direction only)")
    else:
        print("  -> local MAXIMUM")

    # Expected: eigenvalues approximately -9.5e-5 and +1.0e-5 -> SADDLE
