# =====================================================================
# 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: recomputes the headline numbers (PASS/FAIL),
#   plus the Eq 5 = Eq 6 identity and a dimensional check.
# =====================================================================
"""
INDEPENDENT FROM-SCRATCH VERIFICATION of the Trap Number analysis.
Written fresh (not importing the paper's code) to check the math and simulation.

Physical model
--------------
Trap number at a surface point P:
    T(P) = G * M_local(P) / (g_surface * d_eff^2)
where M_local(P) = sum of masses of boulders within d_eff of P.

Each check prints PASS/FAIL against the manuscript's claimed values.
"""
import numpy as np, pandas as pd
from scipy.spatial import cKDTree

# ---- physical constants and adopted Ryugu parameters ----
G        = 6.674e-11          # gravitational constant  [m^3 kg^-1 s^-2]
M_ast    = 4.5e11             # Ryugu mass              [kg]
R        = 450.0              # volume-equivalent radius[m]
g_surface= G*M_ast/R**2       # surface gravity GM/R^2  [m s^-2]
rho_b    = 1600.0             # boulder density         [kg m^-3]
d_eff    = 50.0              # interaction radius      [m]
A_ELL, C_ELL = 502.0, 437.744 # Ryugu ellipsoid semi-axes [m]

print("="*64)
print("STEP 0  adopted parameters")
print(f"  g_surface = G*M/R^2 = {g_surface:.4e} m/s^2   (manuscript 1.48e-4)")
assert abs(g_surface-1.48e-4) < 2e-6, "g_surface mismatch"
print("  PASS: g_surface reproduces 1.48e-4")

# ---- load the real catalogue ----
df = pd.read_csv('Ryugu_Boulder_Attributes.csv')
df = df[df['Length_m'] > 0].copy()            # drop zero-length entries
lon = np.radians(df['Longitude_X'].values % 360.0)
lat = np.radians(df['Latitude_Y'].values)
L   = df['Length_m'].values                   # long-axis length [m]
print(f"\nSTEP 1  catalogue loaded: {len(L)} boulders, largest long axis = {L.max():.1f} m")

def positions(lat, lon):
    """3-D positions on the oblate ellipsoid (geocentric latitude)."""
    r = 1.0/np.sqrt((np.cos(lat)/A_ELL)**2 + (np.sin(lat)/C_ELL)**2)
    return np.column_stack([r*np.cos(lat)*np.cos(lon),
                            r*np.cos(lat)*np.sin(lon),
                            r*np.sin(lat)])

def mass(diam):
    """Sphere mass from diameter."""
    return (4.0/3.0)*np.pi*rho_b*(diam/2.0)**3

def trap_field(diam, d):
    """Trap number (sin a = 1) at every boulder's location for interaction radius d."""
    pos = positions(lat, lon)
    m   = mass(diam)
    tree= cKDTree(pos)
    Mloc= np.array([m[idx].sum() for idx in tree.query_ball_point(pos, d)])
    return G*Mloc/(g_surface*d**2), pos, m

# ---- STEP 2: headline T_max from long-axis diameters ----
T_long, pos, m = trap_field(L, d_eff)
imax = int(np.argmax(T_long))
print("\nSTEP 2  headline trap number (long-axis diameter, d_eff=50 m)")
print(f"  T_max = {T_long[imax]:.3f}  at lon={np.degrees(lon[imax])%360:.0f}E "
      f"lat={np.degrees(lat[imax]):.0f}  (boulder L={L[imax]:.1f} m)")
print(f"  manuscript: 0.67 long-axis, at Otohime (~279E, ~-67 lat, 164 m)")
ok = abs(T_long[imax]-0.67) < 0.03
print("  PASS" if ok else "  FAIL", ": T_max reproduces ~0.67")

# ---- STEP 3: shape-corrected (effective diameter = 0.8 * long axis) ----
T_shape,_,_ = trap_field(0.8*L, d_eff)
print(f"\nSTEP 3  shape-corrected (d_eff_diam = 0.8 x long axis)")
print(f"  T_max = {T_shape.max():.3f}   manuscript shape-corrected: 0.34")
print("  PASS" if abs(T_shape.max()-0.34) < 0.03 else "  FAIL", ": reproduces ~0.34")

# ---- STEP 4: d_eff sensitivity (expect 1.9, 0.67, 0.31, 0.19) ----
print("\nSTEP 4  d_eff sensitivity (long-axis)")
exp = {30:1.9, 50:0.67, 75:0.31, 100:0.19}
allok=True
for d in (30,50,75,100):
    Td,_,_ = trap_field(L, float(d))
    val=Td.max(); good=abs(val-exp[d])<0.06; allok&=good
    print(f"  d_eff={d:3d} m : T_max={val:.2f}  (manuscript {exp[d]})  {'ok' if good else 'OFF'}")
print("  PASS: d_eff sweep reproduced" if allok else "  FAIL")
# check inverse-square law: T*d^2 should be ~constant (Otohime-dominated)
prod=[ (trap_field(L,float(d))[0].max())*d**2 for d in (30,50,75,100)]
print(f"  T_max * d_eff^2 = {[f'{p:.0f}' for p in prod]}  (constant => exact 1/d_eff^2 scaling)")

# ---- STEP 5: analytic scaling Eq 6 vs direct Eq 5 at the max point ----
Mloc_max = m[ [i for i in cKDTree(pos).query_ball_point(pos[imax], d_eff)] ].sum()
T_eq5 = G*Mloc_max/(g_surface*d_eff**2)
T_eq6 = (Mloc_max/M_ast)*(R/d_eff)**2
print("\nSTEP 5  analytic scaling check (Eq 6 must equal Eq 5)")
print(f"  Eq 5 direct      : {T_eq5:.4f}")
print(f"  Eq 6 (Mloc/Mast)(R/deff)^2 : {T_eq6:.4f}")
print("  PASS" if abs(T_eq5-T_eq6)<1e-6 else "  FAIL", ": Eq 6 is an exact rearrangement of Eq 5")

# ---- STEP 6: dimensional check ----
# [G][M]/([g][d^2]) = (m^3 kg^-1 s^-2)(kg)/((m s^-2)(m^2)) = dimensionless
print("\nSTEP 6  dimensional analysis")
num = "m^3 kg^-1 s^-2 * kg = m^3 s^-2"
den = "m s^-2 * m^2 = m^3 s^-2"
print(f"  numerator G*M_local -> {num}")
print(f"  denominator g*d^2   -> {den}")
print("  PASS: ratio is dimensionless (m^3 s^-2 / m^3 s^-2)")

# ---- STEP 7: Otohime single-boulder sanity (mass and T) ----
print("\nSTEP 7  Otohime single-boulder cross-check")
mO = mass(164.0)
print(f"  m(164 m sphere, rho=1600) = {mO:.3e} kg")
print(f"  T from Otohime alone = G m /(g d^2) = {G*mO/(g_surface*d_eff**2):.3f}")
print("  (matches Step 2: Otohime dominates the local mass)")
print("="*64)
print("VERIFICATION COMPLETE")
