#!/usr/bin/env python3
"""
particle-bench — G's particle-scale predictions vs the measured particle world.

What this is: Parts I-III of G's *0 Predictions* (Appendix B), reproduced
EXACTLY as published — three closed-form, parameter-free predictions of
particle-scale quantities, confronted with the values that real precision
experiments produced (CODATA 2018/2022).

What this is NOT: a simulation of particle collisions, and it has no
affiliation with CERN or any accelerator programme. It is arithmetic vs
measurement — nothing more is claimed.

The three claims under test (formulas verbatim from the published verify.py):

  Part I   m_p/m_e     = 21^2*4 + 21*3 + 3^2 + a*21*(1 - 1/(84*pi))
                         + a^2*21*16/1836
  Part II  G           = a^21 * (1 + 1/pi) * hbar*c / m_e^2
  Part III (m_n-m_p)/m_e = 3*(1 - 1/(2*pi)) + a*(1 + 1/(2*pi))

where a is the fine-structure constant — the single measured dimensionless
input. Part II additionally uses CODATA m_e (and hbar, c — exact by SI
definition); the boundary is stated so a sceptic doesn't have to find it.

Gate (self-verifying; exits non-zero if any residual exceeds G's stated
tolerance — outside tolerance = the Part is dead):

  Part I   |residual| <= 0.017 ppb
  Part II  |residual| <= 1.0 %
  Part III |residual| <= 5.0 ppm

Source: G (Studio G), *0 Predictions*, Appendix B. DOI 10.5281/zenodo.19208226.
Licence of the source formulas: CC BY 4.0. Python stdlib only, < 1 s, offline.
"""

import math
import sys

# ---------------------------------------------------------------------------
# Measured inputs — exactly as in G's published script. Vintages are mixed and
# labelled per value: this alpha is the CODATA 2018 value (2022 revised the
# last digits to ...5643); the Part I comparison ratio is CODATA 2022; G's
# central value is identical in 2018 and 2022. Values stay verbatim as G
# published them — labels state which table each number came from.
# ---------------------------------------------------------------------------
ALPHA = 7.2973525693e-3    # fine-structure constant (dimensionless, CODATA 2018)
HBAR = 1.054571817e-34     # reduced Planck constant (J*s)
C = 2.99792458e8           # speed of light (m/s) — exact by SI definition
M_E = 9.1093837015e-31     # electron mass (kg)

# G's stated tolerances — the kill switches for each Part.
TOL_I_PPB = 0.017          # Part I:   proton/electron mass ratio
TOL_II_PCT = 1.0           # Part II:  gravitational constant
TOL_III_PPM = 5.0          # Part III: neutron-proton mass difference


def part_i():
    """Proton-to-electron mass ratio. Measured: CODATA 2022."""
    static = 21**2 * 4 + 21 * 3 + 3**2                    # = 1836
    first_order = ALPHA * 21 * (1 - 1 / (84 * math.pi))
    second_order = ALPHA**2 * 21 * 16 / 1836
    predicted = static + first_order + second_order
    measured = 1836.152673426
    residual_ppb = (predicted - measured) / measured * 1e9
    return predicted, measured, residual_ppb


def part_ii():
    """Gravitational constant. Measured: CODATA 2018/2022 (identical central value)."""
    alpha_g = ALPHA**21 * (1 + 1 / math.pi)
    predicted = alpha_g * HBAR * C / M_E**2
    measured = 6.67430e-11
    residual_pct = (predicted - measured) / measured * 100
    return predicted, measured, residual_pct


def part_iii():
    """Neutron-proton mass difference, in electron masses. Measured: CODATA."""
    static = 3 * (1 - 1 / (2 * math.pi))
    dynamic = ALPHA * (1 + 1 / (2 * math.pi))
    predicted = static + dynamic
    measured = 2.5309883
    residual_ppm = (predicted - measured) / measured * 1e6
    return predicted, measured, residual_ppm


def main():
    fails = 0

    print("particle-bench — G's particle-scale predictions vs measurement")
    print("(no collisions simulated; no CERN affiliation; arithmetic vs CODATA)")
    print("=" * 68)

    # Part I ----------------------------------------------------------------
    p, m, r = part_i()
    ok = abs(r) <= TOL_I_PPB
    fails += not ok
    print("Part I   — proton-to-electron mass ratio  m_p/m_e")
    print("  formula  : 21^2*4 + 21*3 + 3^2 + a*21*(1-1/(84pi)) + a^2*21*16/1836")
    print(f"  predicted: {p:.9f}")
    print(f"  measured : {m:.9f}   (CODATA 2022)")
    print(f"  residual : {r:.3f} ppb   tolerance {TOL_I_PPB} ppb   "
          f"{'PASS' if ok else 'DEAD'}")
    print()

    # Part II ---------------------------------------------------------------
    p, m, r = part_ii()
    ok = abs(r) <= TOL_II_PCT
    fails += not ok
    print("Part II  — gravitational constant  G")
    print("  formula  : a^21 * (1 + 1/pi) * hbar*c / m_e^2")
    print(f"  predicted: {p:.4e} N*m^2/kg^2")
    print(f"  measured : {m:.4e}   (CODATA 2018/2022)")
    print(f"  residual : {r:.2f} %   tolerance {TOL_II_PCT} %   "
          f"{'PASS' if ok else 'DEAD'}")
    print()

    # Part III --------------------------------------------------------------
    p, m, r = part_iii()
    ok = abs(r) <= TOL_III_PPM
    fails += not ok
    print("Part III — neutron-proton mass difference  (m_n - m_p)/m_e")
    print("  formula  : 3*(1 - 1/(2pi)) + a*(1 + 1/(2pi))")
    print(f"  predicted: {p:.10f}")
    print(f"  measured : {m:.10f}   (CODATA)")
    print(f"  residual : {r:.2f} ppm   tolerance {TOL_III_PPM} ppm   "
          f"{'PASS' if ok else 'DEAD'}")
    print()

    print("=" * 68)
    if fails:
        print(f"FAIL — {fails} of 3 particle-scale predictions outside "
              f"G's stated tolerance.")
        return 1
    print("PASS — all 3 particle-scale predictions inside G's stated "
          "tolerances.")
    print("(Reproducing the arithmetic is internal consistency, not proof "
          "that nature agrees.)")
    return 0


if __name__ == "__main__":
    sys.exit(main())