#!/usr/bin/env python3
"""
================================================================================
monte_carlo.py
Monte Carlo error analysis for the trap number 𝒯.

Runs N_REALIZATIONS independent realizations of the boulder population for each
asteroid. Named boulders are held fixed; the smaller population is resampled.
Reports 𝒯_max as mean +/- standard deviation.

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

Usage:
    python monte_carlo.py
================================================================================
"""

import numpy as np
from code_trap_number import ASTEROIDS, compute_trap_number_map

N_REALIZATIONS = 100


def monte_carlo(params, n=N_REALIZATIONS, d_eff=50.0):
    """Run n realizations and return arrays of T_max and surface fraction > 0.1."""
    lat_grid = np.linspace(-85, 85, 40)
    lon_grid = np.linspace(0, 357, 72)
    Tmaxes, fracs = [], []
    for seed in range(n):
        Tmap, _ = compute_trap_number_map(params, lat_grid, lon_grid,
                                           d_eff=d_eff, seed=seed)
        Tmaxes.append(Tmap.max())
        fracs.append(100 * np.mean(Tmap > 0.1))
    return np.array(Tmaxes), np.array(fracs)


if __name__ == '__main__':
    print(f"Monte Carlo: {N_REALIZATIONS} realizations per asteroid, d_eff = 50 m\n")
    print(f"{'Asteroid':<10}{'T_max':>16}{'range':>20}{'%>0.1':>14}")
    print("-" * 60)
    for name, params in ASTEROIDS.items():
        Tmaxes, fracs = monte_carlo(params)
        print(f"{name:<10}"
              f"{Tmaxes.mean():>8.3f} +/- {Tmaxes.std():<5.3f}"
              f"{f'[{Tmaxes.min():.3f}, {Tmaxes.max():.3f}]':>18}"
              f"{f'{fracs.mean():.1f} +/- {fracs.std():.1f}':>14}")

    # Expected output:
    #   Ryugu     0.44 +/- 0.05   [0.42, 0.61]    1.8 +/- 0.4   (synthetic; real ~= 0.5)
    #   Bennu     0.038 +/- 0.006 [0.036, 0.079]  0.0 +/- 0.0
    #   Itokawa   0.051 +/- 0.001 [0.051, 0.055]  0.0 +/- 0.0
