# =============================================================================
# code_predicted_estimates.py
# Documents the BASIS of the predicted (non-visited) trap numbers plotted as
# stars in Figures 3 and 5. Unlike the three visited asteroids, these points are
# NOT computed from a boulder catalogue; they are order-of-magnitude predictions
# for the small-body regime, included only to indicate where such bodies fall.
#
#   Dimorphos (~150 m across, post-DART):  T ~ 2-3
#       Taken from Robin et al. (2024), who characterised the dense post-impact
#       boulder field. The figures plot the representative value T = 2.3. This is
#       a cited estimate, not recomputed here. Hera (2026) will allow a direct
#       catalogue computation of T for Dimorphos.
#
#   50 m NEA (representative sub-100 m body):  T ~ a few to ~10
#       Estimated from the single-boulder form, Equation (7a):
#           T = (rho_b/rho) * D^3 / (8 * R * d_eff^2)
#       For a small rubble pile, d_eff (the regolith-transport scale) becomes a
#       few metres rather than 50 m, and a single block can be a large fraction
#       of the body. With the round assumptions below the estimate lands at a few;
#       the figures plot T = 6.8 as a representative value in this range. This is
#       an extrapolation, NOT a catalogue computation, and is presented as such.
#       (For a still smaller ~30 m body such as 1998 KY26, the same scaling gives
#       T ~ 10-100, as stated in Section 4.5.3.)
#
# Bottom line: the visited-body values (Table 1, Figs 1-2) are computed; the
# predicted points (Figs 3, 5 stars) are estimates with the bases above.
# =============================================================================
import numpy as np

G = 6.674e-11

def T_single_boulder(rho_b_over_rho, D, R, d_eff):
    """Equation (7a): dimensionless single-dominant-boulder trap number."""
    return rho_b_over_rho * D**3 / (8.0 * R * d_eff**2)

def main():
    print("Predicted small-body trap numbers (illustrative, not catalogue values)\n")

    print("Dimorphos (post-DART): T ~ 2-3  (Robin et al. 2024); figures plot 2.3.\n")

    # Representative 50 m NEA: radius 25 m, one ~20 m dominant block,
    # fine-regolith interaction scale d_eff ~ 3 m, boulder/bulk density ratio 1.4
    R = 25.0; D = 20.0; d_eff = 3.0; ratio = 1.4
    T = T_single_boulder(ratio, D, R, d_eff)
    print("50 m NEA example  (R=%.0f m, dominant block D=%.0f m, d_eff=%.0f m, "
          "rho_b/rho=%.1f):" % (R, D, d_eff, ratio))
    print("    Equation (7a) -> T ~ %.1f  (figures plot 6.8 as representative)\n" % T)

    # Smaller 30 m body (e.g. 1998 KY26): even higher T. The result is very
    # sensitive to the assumed regolith scale d_eff; a fine d_eff ~ 1.5 m gives:
    R = 15.0; D = 14.0; d_eff = 1.5; ratio = 1.4
    T2 = T_single_boulder(ratio, D, R, d_eff)
    print("30 m body example (R=%.0f m, dominant block D=%.0f m, d_eff=%.1f m, "
          "rho_b/rho=%.1f):" % (R, D, d_eff, ratio))
    print("    Equation (7a) -> T ~ %.0f  (in the stated 10-100 range; strongly "
          "d_eff-dependent)" % T2)

    print("\nThese are order-of-magnitude predictions for the sub-100 m regime, "
          "not\ncatalogue computations; they will be tested by Hera, OSIRIS-APEX, "
          "and Hayabusa2.")

if __name__ == '__main__':
    main()
