#!/usr/bin/env python3
"""
verify_P257.py — Verifier for Addendum 257 (fold-native formation mechanism).

Sections:
  S1  Lambda-negligibility at J1          — checks 1-5
  S2  Abundance/amplitude requirements    — checks 6-11
  S3  SIGW companion                      — checks 12-15

Copyright: Leon Fernando Vlegels - MIT
"""
import math, sys

PASS = FAIL = 0
def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")

PI = math.pi
L0 = 1 - PI**2/32
OM = 0.315
RHO_CRIT_GEV4 = 3.67e-47
MSUN_G = 1.989e33

# S1
print("S1  Lambda-negligibility at J1")
def ratio_rad(T_GeV, g=106.75):
    return 2*L0*RHO_CRIT_GEV4 / ((PI**2/30)*g*T_GeV**4)
def ratio_mat(z):
    return 2*L0 / (OM*(1+z)**3)
r3, r6 = ratio_rad(3e6), ratio_rad(6e6)
check(1, "L0 = 1 - pi^2/32 = 0.6915749...", abs(L0 - 0.6915749) < 1e-6)
check(2, "C_lna 3 PeV: rho_DE/rho_rad = %.2e < 1e-70" % r3, r3 < 1e-70)
check(3, "C_lna 6 PeV smaller still (T^-4 scaling)", r6 < r3 and abs(r6/r3 - 1/16) < 0.01)
z17, z26 = ratio_mat(17), ratio_mat(26)
check(4, "C_eta z=17,26: rho_DE/rho_m in [1e-4, 1e-3]", 1e-4 < z26 < z17 < 1e-3)
check(5, "Lambda-mediated threshold shift negligible at BOTH survivors' J1 "
         "(max fraction %.1e)" % z17, z17 < 1e-3)

# S2
print("S2  Abundance/amplitude requirements (M = 1.1e24 g)")
M = 1.1e24
def beta_req(f): return 6.4e-9 * f * math.sqrt(M/MSUN_G)
def sigma_req(beta, dc=0.45):
    lo, hi = 1e-3, 1.0
    for _ in range(200):
        mid = 0.5*(lo+hi)
        if math.erfc(dc/(math.sqrt(2)*mid)) < beta: lo = mid
        else: hi = mid
    return 0.5*(lo+hi)
b1 = beta_req(1.0)
s1v, s2v, s3v = sigma_req(b1), sigma_req(beta_req(1e-2)), sigma_req(beta_req(1e-4))
check(6, "beta(f=1) = %.2e in [1e-13, 2e-13]" % b1, 1e-13 < b1 < 2e-13)
check(7, "beta linear in f (4 decades)", abs(beta_req(1e-4)/b1 - 1e-4) < 1e-12)
check(8, "sigma(f=1) = %.4f in [0.055, 0.065]" % s1v, 0.055 < s1v < 0.065)
check(9, "near-binarity: sigma(f=1)/sigma(f=1e-4) = %.3f < 1.2" % (s1v/s3v), s1v/s3v < 1.2)
check(10, "boost over inflationary baseline (5e-5) > 1000x", s1v/5e-5 > 1000)
check(11, "erfc round-trip: recovered beta within 1%%",
      abs(math.erfc(0.45/(math.sqrt(2)*s1v))/b1 - 1) < 0.01)

# S3
print("S3  SIGW companion")
def f_peak(M_g): return 3e-9 * (M_g/(30*MSUN_G))**-0.5
f_lo, f_mid, f_hi = f_peak(2.6e24), f_peak(1.1e24), f_peak(6.0e23)
check(12, "f_peak(1.1e24 g) = %.2e Hz in LISA band [1e-4, 1e-1]" % f_mid,
      1e-4 < f_mid < 1e-1)
check(13, "entire C_lna band [6e23, 2.6e24] g inside LISA", 1e-4 < f_lo < f_hi < 1e-1)
check(14, "f_peak monotone decreasing in M", f_lo < f_mid < f_hi)
A2 = (s1v**2)**2
check(15, "Omega_GW amplitude proxy sigma^4 = %.1e > 1e-7 (narrow-feature "
          "LISA-detectable scale)" % A2, A2 > 1e-7)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)