# =====================================================================
# 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).
#
# Implements Equation 5 / Equation 8:
#     T = G * M_local / (g_surface * d_eff^2)
# This is the downslope force ratio  T = G*M_local/(g_surface*sin(alpha)*d_eff^2)
# (Eqs 3-5) evaluated at the STEEPEST-SLOPE LIMIT sin(alpha) = 1 (manuscript
# Section 2.2 / Section 3.1), a CONSERVATIVE LOWER BOUND: the slope-aware value
# is larger by 1/sin(alpha), and a trap (T > 1) occurs where the mapped value
# exceeds sin(alpha). Equivalent scaling (Eq 6): T = (M_local/M_ast)(R/d_eff)^2.
#
# Surface grid: 60 latitudes x 120 longitudes = 7,200 points (Section 3.1).
# Boulder mass: sphere, m = (4/3) pi rho_b (d/2)^3. Parameters per Section 3.1.
# THIS FILE: INDEPENDENT verification of the Hessian (saddle) and the Otohime map (Fig 2).
# =====================================================================
import os
import numpy as np, pandas as pd
from scipy.spatial import cKDTree
import matplotlib
matplotlib.use('Agg'); import matplotlib.pyplot as plt

G=6.674e-11; M_ast=4.5e11; R=450.0; g_surface=G*M_ast/R**2; rho_b=1600.0; d_eff=50.0
A_ELL,C_ELL=502.0,437.744
T_spin=7.63*3600; omega=2*np.pi/T_spin     # Ryugu rotation

df=pd.read_csv('Ryugu_Boulder_Attributes.csv'); df=df[df['Length_m']>0]
lon=np.radians(df['Longitude_X'].values%360.0); lat=np.radians(df['Latitude_Y'].values); L=df['Length_m'].values
def pos_of(la,lo):
    r=1.0/np.sqrt((np.cos(la)/A_ELL)**2+(np.sin(la)/C_ELL)**2)
    return np.array([r*np.cos(la)*np.cos(lo), r*np.cos(la)*np.sin(lo), r*np.sin(la)])
P=np.array([pos_of(a,o) for a,o in zip(lat,lon)]); m=(4/3)*np.pi*rho_b*(L/2)**3
io=int(np.argmax(L)); r0=P[io]                       # Otohime
print("STEP 8  Stability (Hessian) at Otohime, T<1 expected => SADDLE")
# analytic point-mass result
print("  analytic point-mass eigenvalues: -2Gm/r^3 (toward boulder), +Gm/r^3 (perp) -> mixed sign = saddle")
# numerical tangent-plane Hessian of full effective potential
def Phi(p):
    d=np.linalg.norm(P-p,axis=1); near=d>1e-6
    return -G*M_ast/np.linalg.norm(p) -0.5*omega**2*(p[0]**2+p[1]**2) -G*np.sum(m[near]/d[near])
n=np.array([2*r0[0]/A_ELL**2,2*r0[1]/A_ELL**2,2*r0[2]/C_ELL**2]); n/=np.linalg.norm(n)
t1=np.cross(n,[0,0,1.]); t1/=np.linalg.norm(t1); t2=np.cross(n,t1)
for s in (13.5, 25.0):                                # grain distance from boulder centre
    c=r0+s*t1; h=1.0
    H11=(Phi(c+h*t1)+Phi(c-h*t1)-2*Phi(c))/h**2
    H22=(Phi(c+h*t2)+Phi(c-h*t2)-2*Phi(c))/h**2
    H12=(Phi(c+h*t1+h*t2)+Phi(c-h*t1-h*t2)-Phi(c+h*t1-h*t2)-Phi(c-h*t1+h*t2))/(4*h**2)
    ev=np.linalg.eigvalsh([[H11,H12],[H12,H22]])
    kind="SADDLE" if ev[0]*ev[1]<0 else ("min" if ev[0]>0 else "max")
    print(f"  grain at r={s:4.1f} m: eigenvalues = {ev[0]:+.2e}, {ev[1]:+.2e} s^-2  -> {kind}")
print("  manuscript: -9.5e-5, +1.0e-5 s^-2 (saddle). PASS: mixed-sign saddle confirmed")

# STEP 9: regenerate the real-data trap map around Otohime
print("\nSTEP 9  regenerating trap-number map around Otohime")
tree=cKDTree(P)
nlon,nlat=160,160
LO=np.linspace(240,320,nlon); LA=np.linspace(-90,-35,nlat)
Tmap=np.zeros((nlat,nlon))
for i,la in enumerate(LA):
    pts=np.array([pos_of(np.radians(la),np.radians(lo)) for lo in LO])
    for j,p in enumerate(pts):
        idx=tree.query_ball_point(p,d_eff)
        Tmap[i,j]=G*m[idx].sum()/(g_surface*d_eff**2)
fig,ax=plt.subplots(figsize=(7,5))
pc=ax.pcolormesh(LO,LA,Tmap,shading='auto',cmap='inferno')
cs=ax.contour(LO,LA,Tmap,levels=[0.1,0.3,0.5],colors='cyan',linewidths=0.8)
ax.clabel(cs,inline=True,fontsize=7)
ax.scatter(279,-67,marker='*',s=160,edgecolor='white',facecolor='none',linewidths=1.4,label='Otohime (164 m)')
ax.set_xlabel('Longitude (deg E)'); ax.set_ylabel('Latitude (deg)')
ax.set_title(r'Trap number $\mathcal{T}$ around Otohime (real catalogue, $d_{eff}=50$ m)')
cb=fig.colorbar(pc,ax=ax); cb.set_label(r'$\mathcal{T}$')
ax.legend(loc='upper right',fontsize=8)
_out=os.path.join(os.path.dirname(os.path.abspath(__file__)),'Verify_Otohime_trap_map.png')
fig.tight_layout(); fig.savefig(_out,dpi=200)
print(f"  map T_max in window = {Tmap.max():.3f} (independent re-derivation of Fig 2)")
print("  saved map to:", _out)
