# Code — The Trap Number 𝒯

Companion code for **"The Trap Number 𝒯: A Dimensionless Criterion for
Boulder–Boulder Gravitational Dynamics on Near-Earth Rubble-Pile Asteroid
Surfaces"** (prepared for submission).
This file maps each script to the equation, figure, or number it reproduces.

## What the code computes

The trap number is computed in the form given by **Equation 5**:

    𝒯 = G · M_local / (g_surface · d_eff²)

where M_local is the total boulder mass within d_eff of a surface point,
g_surface = GM/R² is the surface gravity, and d_eff = 50 m. This is the downslope
force ratio 𝒯 = G·M_local/(g_surface·sinα·d_eff²) (Eqs 3–5) evaluated at the
**steepest-slope limit sinα = 1** (manuscript Section 2.2 / 3.1), a conservative
lower bound; the slope-aware value is larger by 1/sinα. The equivalent scaling is
**Equation 6**, 𝒯 = (M_local/M_ast)(R/d_eff)². The surface is sampled on a
60 (latitude) × 120 (longitude) = 7,200-point grid (Section 3.1). 𝒯 > 1 means the
boulder pull exceeds the downslope drive (a convergence zone); 𝒯 = 1 is a nominal order-unity scale, not a proven bifurcation (Section 3.5).

## Adopted parameters (Section 3.1)

| Body | M (kg) | R (m) | g_surface (m/s²) | boulder ρ (kg/m³) |
|------|--------|-------|------------------|-------------------|
| Ryugu | 4.5 × 10¹¹ | 450 | 1.48 × 10⁻⁴ | 1600 |
| Bennu | 7.329 × 10¹⁰ | 245 | 8.2 × 10⁻⁵ | 1600 |
| Itokawa | 3.51 × 10¹⁰ | 165 | 8.6 × 10⁻⁵ | 2500 |

G = 6.674 × 10⁻¹¹ m³ kg⁻¹ s⁻²; d_eff = 50 m (default).

## Script → result map

| Script | Reproduces | Expected output |
|--------|-----------|-----------------|
| `code_trap_number.py` | Eq 5; synthetic model (Sec 3.2) | Ryugu 𝒯_max ≈ 0.42, Bennu 0.036, Itokawa 0.051 |
| `code_fig2_trap_map.py` | **Figure 2** (3-panel real-catalogue trap map: field + zoom + 1/r² panel) | 𝒯 ≈ 0.5 (0.34–0.67) at Otohime; bright disk = d_eff footprint, not a trap |
| `code_sfd_slope.py` | Ryugu size-frequency slope | cumulative slope -2.5 to -2.9 over the complete range; adopted -2.665 |
| `code_predicted_estimates.py` | basis of the predicted points in Figs 3, 5 | Dimorphos T~2.3 (Robin 2024); 50 m NEA T~6 from Eq 7a (illustrative) |
| `code_real_ryugu_trap_number.py` | real catalogue (Sec 3.2–3.3, 3.6) | 𝒯_max ≈ 0.5 (0.34–0.67) at Otohime; d_eff 30/50/75/100 m → 1.9/0.67/0.31/0.19 |
| `code_monte_carlo.py` | Sec 3.4 Monte Carlo | Ryugu 0.44 ± 0.05, Bennu 0.038, Itokawa 0.051 |
| `code_hessian_stability.py` | Sec 3.5 stability | saddle; eigenvalues −9.5 × 10⁻⁵, +1.0 × 10⁻⁵ s⁻² |
| `code_deff_sensitivity.py` | Sec 3.6 | 𝒯_max vs d_eff |
| `code_generate_figures.py` | Figures 1–6 | six PNGs saved into this folder |
| `code_stability_diagram.py` | Figure 6 (schematic stability diagram) | illustrative bowl, saddle, and nominal 𝒯≈1 scale (not a proven bifurcation) |
| `run_trap_number_real_catalogue.py` | generic driver | 𝒯_max ≈ 0.67 at the Otohime location |
| `verify_from_scratch.py` | independent recheck | every line PASS, "VERIFICATION COMPLETE" |
| `hessian_and_plot.py` | independent Hessian + map | saddle confirmed; map PNG saved here |
| `TrapNumber_Derivations.wl` | symbolic derivation | Equations 1–8 (Mathematica) |

## How to run and check (macOS / Linux)

    python3 -m pip install -r requirements.txt
    python3 verify_from_scratch.py

`verify_from_scratch.py` reloads the raw catalogue and recomputes g_surface, the
headline 𝒯_max (0.670 long-axis / 0.343 shape-corrected), the d_eff sweep, the
1/d_eff² scaling, the Eq 5 = Eq 6 identity, and a dimensional check, printing
PASS/FAIL for each. The figure scripts save their PNGs into this folder and print
the path.

## Data provenance

`Ryugu_Boulder_Attributes.csv` (also provided as `ryugu_boulders.csv` with simple
column names) is the published Ryugu boulder catalogue of **Ray et al. (2026),
*Scientific Reports* 16, 14404, https://doi.org/10.1038/s41598-026-40370-4**,
reused under CC BY 4.0 (columns: FID, Longitude_X, Latitude_Y, Length_m,
Orientation, Aspect). The synthetic asteroid-parameter and catalogue CSVs are
generated to the parameters in the table above.
