The code

The exact source of every calculation on this site, ready to copy and paste.

Select the text inside any box and copy it. It is plain text with the original line breaks and indentation, and it has no line numbers, so it pastes straight into a file. To run a script, save it, then run python -I name.py in an empty folder. The scripts need Python with numpy and scipy. Scripts that read the workbook expect TRPG_Every_Particle_to_Every_Particle.xlsx in the working folder (it is on the Data & code page). The math page is produced by build_math_page.py from math_inputs.json, with site_common.py beside them; the numbers on the math page come from running that code. Every download also has a SHA-256 fingerprint here and on the Data & code page, so you can check you hold the same bytes.

"""Finishing calculations (14 Nov 2026).
A) Area law: entanglement of a region of a wave medium (2D lattice of coupled oscillators, ground state) grows with its BOUNDARY, not its volume
   (Bombelli, Koul, Lee & Sorkin 1986; Srednicki 1993). This is the input Jacobson (1995) needs to derive Einstein's equations.
B) Charges from consistency: anomaly cancellation fixes the Standard Model charges and the number of colours."""
import numpy as np, math
N=36; m=0.01
idx=lambda i,j:i*N+j
K=np.zeros((N*N,N*N))
for i in range(N):
    for j in range(N):
        a=idx(i,j); K[a,a]=4+m*m
        for di,dj in((1,0),(0,1)):
            ii,jj=i+di,j+dj
            if ii<N and jj<N: b=idx(ii,jj); K[a,b]=K[b,a]=-1
w,U=np.linalg.eigh(K)
X=0.5*(U/np.sqrt(w))@U.T; P=0.5*(U*np.sqrt(w))@U.T
res=[]
for L in range(2,15,2):
    o=(N-L)//2; A=[idx(i,j) for i in range(o,o+L) for j in range(o,o+L)]
    ev=np.linalg.eigvals(X[np.ix_(A,A)]@P[np.ix_(A,A)]).real; nu=np.sqrt(np.clip(ev,0.25+1e-12,None))
    S=sum((n+.5)*math.log(n+.5)-(n-.5)*math.log(n-.5) for n in nu if n>0.5+1e-9)
    res.append((L,4*L,L*L,S)); print("block %2dx%-2d boundary %3d volume %4d  entanglement %.3f  per boundary %.4f"%(L,L,4*L,L*L,S,S/(4*L)))
L=np.array([r[0] for r in res]); S=np.array([r[3] for r in res])
cb=np.polyfit(4*L,S,1); cv=np.polyfit(L*L,S,1)
rb=np.sqrt(np.mean((np.polyval(cb,4*L)-S)**2)); rv=np.sqrt(np.mean((np.polyval(cv,L*L)-S)**2))
print("fit to boundary: error %.3f ; fit to volume: error %.3f -> %s law"%(rb,rv,"AREA (boundary)" if rb<rv else "volume"))
print("entanglement per unit boundary ~ %.3f -> Jacobson/Bekenstein: G = a^2/(4 x that) in medium units: G ~ %.2f a^2 (gravity strength set by the medium's spacing a)"%(cb[0],1/(4*cb[0])))
# B) charges
print("\nB) Charges from consistency (anomaly cancellation; Geng & Marshak 1989; Minahan, Ramond & Warner 1990; Bar & Wiese 2001).")
print("   One generation with N colours: quark doublet Q(N,2), u(N), d(N), lepton doublet L(2), e(1); fix the scale with Y_L = -1/2.")
import sympy as sp
YQ, Yu, Yd, Ye, N = sp.symbols("Y_Q Y_u Y_d Y_e N")
YL = sp.Rational(-1, 2)
for Nc in (1, 2, 3, 4, 5):
    eqs = [2*YQ - Yu - Yd,                                   # SU(N)^2 U(1)
           Nc*YQ + YL,                                       # SU(2)^2 U(1)
           Nc*(2*YQ - Yu - Yd) + 2*YL - Ye,                  # gravity^2 U(1)
           Nc*(2*YQ**3 - Yu**3 - Yd**3) + 2*YL**3 - Ye**3]   # U(1)^3
    sols = sp.solve(eqs, [YQ, Yu, Yd, Ye], dict=True)
    for so in sols:
        qu = sp.Rational(1, 2) + so[YQ]; qd = -sp.Rational(1, 2) + so[YQ]
        if so[Yu] > so[Yd]:
            print(f"   N = {Nc}: up-quark charge {qu}, down-quark charge {qd}, electron {so[Ye] + 0}, neutrino 0")
print("   Result: for ANY number of colours N the anomalies cancel, with up = 1/2 + 1/(2N) and down = -1/2 + 1/(2N).")
print("   So the charges are forced ONCE N is known, and N = 3 gives exactly +2/3 and -1/3. N itself is not forced by")
print("   anomalies; it is measured (e.g. the neutral pion decay rate and the e+e- -> hadrons ratio both need N = 3).")

"""Antiparticle matrices and decay data (14 Nov 2026), added to TRPG_Every_Particle_to_Every_Particle.xlsx.
Sheets: Antiparticles, AntiAnti_Pairs, Particle_vs_Antiparticle (19x19 = 361), Decays (PDG 2024), Decay_Relations (Fermi-theory check).
Annihilation rule: a pair can annihilate into neutral force carriers if every conserved number sums to zero
(charge, baryon number, each lepton-family number, quark flavour; colour must be able to cancel, so a coloured and a
colourless particle cannot annihilate). Proton-antineutron pairs annihilate into charged pions. It can fuse into a W if the pair
has charge +/-1, zero baryon number and is a weak-doublet partner pair (CKM or lepton-family); otherwise only scattering."""
import pandas as pd, itertools, math
F="TRPG_Every_Particle_to_Every_Particle.xlsx"
P=pd.read_excel(F,sheet_name="Particles")
B={"up":1/3,"down":1/3,"strange":1/3,"charm":1/3,"bottom":1/3,"top":1/3,"proton":1,"neutron":1}
Lf={"electron":"e","nu_e":"e","muon":"mu","nu_mu":"mu","tau":"tau","nu_tau":"tau"}
self_conj={"photon","Z","Higgs","gluon"}   # gluon octet as a set is self-conjugate
CKM={("up","down"):0.97367,("up","strange"):0.22431,("up","bottom"):0.00382,("charm","down"):0.221,("charm","strange"):0.975,
     ("charm","bottom"):0.0408,("top","down"):0.0086,("top","strange"):0.0415,("top","bottom"):0.999}
def qn(name,anti):
    s=-1 if anti and name not in self_conj else 1
    r=P.set_index("name").loc[name]
    L={"e":0,"mu":0,"tau":0}
    if name in Lf: L[Lf[name]]=1*s
    return dict(Q=r.charge*s,B=B.get(name,0)*s,**{"L_"+k:v for k,v in L.items()})
def label(n,anti): return ("anti-"+n if anti and n not in self_conj else n)
anti=P.copy(); anti["name"]=[label(n,True) for n in P.name]; anti["charge"]=[-c if n not in self_conj else c for n,c in zip(P.name,P.charge)]
anti["T3"]=[-t for t in P.T3]; anti["baryon_number"]=[-B.get(n,0) for n in P.name]
def pair_row(a,aa,b,ba):
    A,Bq=qn(a,aa),qn(b,ba); tot={k:A[k]+Bq[k] for k in A}
    ma=P.set_index("name").loc[a].mass_MeV; mb=P.set_index("name").loc[b].mass_MeV
    direct = (a==b and aa!=ba) or (a==b and a in self_conj)
    quarks={"up","down","strange","charm","bottom","top"}; nucleons={"proton","neutron"}
    col=lambda n: bool(P.set_index("name").loc[n].colored)
    if col(a)!=col(b):
        fate="no annihilation; scatter only (colour cannot cancel)"
    elif all(abs(v)<1e-9 for v in tot.values()) and a in quarks and b in quarks and a!=b:
        fate="flavour-changing: no direct annihilation (only rare loop processes)"
    elif all(abs(v)<1e-9 for v in tot.values()):
        fate="can annihilate -> neutral carriers (photons, gluons, Z, Higgs)"
    elif a in nucleons and b in nucleons and abs(tot["B"])<1e-9:
        fate="can annihilate -> pions (strong force; hadronic)"
    elif abs(abs(tot["Q"])-1)<1e-9 and abs(tot["B"])<1e-9 and (((a,b) in CKM or (b,a) in CKM) or (Lf.get(a)==Lf.get(b) and a in Lf)) and all(abs(tot[k])<1e-9 for k in ("L_e","L_mu","L_tau")):
        fate="can fuse -> W boson (weak charged current)"
    else: fate="no annihilation; scatter only"
    qq=A["Q"]*Bq["Q"]
    return {"i":label(a,aa),"j":label(b,ba),"direct_opposites":int(direct),"net_charge":round(tot["Q"],4),"net_baryon":round(tot["B"],4),
            "net_lepton_e":tot["L_e"],"net_lepton_mu":tot["L_mu"],"net_lepton_tau":tot["L_tau"],
            "EM_push_pull":"push (repel)" if qq>0 else ("pull (attract)" if qq<0 else "none"),
            "min_energy_to_make_pair_MeV":ma+mb,"fate":fate}
names=list(P.name)
aa=pd.DataFrame([pair_row(a,True,b,True) for a,b in itertools.combinations(names,2)])
pa=pd.DataFrame([pair_row(a,False,b,True) for a in names for b in names])
# Decays (PDG 2024 values)
dec=[("muon","e nu_mu anti-nu_e",1.0,"lifetime 2.1969811e-6 s"),("tau","mu nu nu",0.1739,"lifetime 290.3e-15 s"),("tau","e nu nu",0.1782,""),("tau","hadrons nu",0.648,""),
     ("W","e nu",0.1071,"width 2.085 GeV"),("W","mu nu",0.1063,""),("W","tau nu",0.1138,""),("W","hadrons (quark pairs)",0.6741,""),
     ("Z","e e",0.03363,"width 2.4955 GeV"),("Z","mu mu",0.03366,""),("Z","tau tau",0.03370,""),("Z","invisible (neutrino pairs)",0.2000,""),("Z","hadrons",0.6991,""),
     ("Higgs","b anti-b",0.582,"SM width ~4.1 MeV"),("Higgs","W W*",0.214,""),("Higgs","gluon gluon",0.082,""),("Higgs","tau tau",0.0627,""),("Higgs","c anti-c",0.0289,""),("Higgs","Z Z*",0.0262,""),("Higgs","photon photon",0.00227,""),
     ("top","W b",0.998,"width ~1.42 GeV"),("neutron","p e anti-nu_e",1.0,"lifetime 878.4 s"),("proton","(stable; > 1e34 yr)",0.0,"")]
decays=pd.DataFrame(dec,columns=["parent","products","branching_fraction","total_lifetime_or_width"])
# Fermi theory check: muon lifetime and tau lifetime from the muon via mass ratio^5 = e^(10 S)
GF=1.1663788e-5*1e-6; hbar=6.582119569e-22  # GF in MeV^-2, hbar in MeV s
mmu,mtau=105.6583755,1776.93
tau_mu=hbar/(GF**2*mmu**5/(192*math.pi**3)); S=0.5*math.log(mtau/mmu)
tau_tau=tau_mu*0.1782*math.exp(-10*S)
rel=pd.DataFrame([{"quantity":"muon lifetime from Fermi theory","value":tau_mu,"measured":2.1969811e-6},
                  {"quantity":"tau lifetime = muon lifetime x BR(tau->e nu nu) x e^(-10 S_tau,mu)","value":tau_tau,"measured":290.3e-15}])
with pd.ExcelWriter(F,mode="a",engine="openpyxl",if_sheet_exists="replace") as w:
    anti.to_excel(w,sheet_name="Antiparticles",index=False); aa.to_excel(w,sheet_name="AntiAnti_Pairs",index=False)
    pa.to_excel(w,sheet_name="Particle_vs_Antiparticle",index=False); decays.to_excel(w,sheet_name="Decays",index=False); rel.to_excel(w,sheet_name="Decay_Relations",index=False)
print(pa.fate.value_counts().to_string()); print("direct opposites:",pa.direct_opposites.sum())
print(rel.to_string(index=False))
print(pa[pa.fate.str.contains("W")][["i","j"]].to_string(index=False))

"""Every-particle-to-every-particle relational data (TRPG), built 13 Nov 2026.
Values: PDG 2024/2025 review (masses, CKM magnitudes). Quark masses: MS-bar (u,d,s at 2 GeV; c,b at own mass); top = pole.
"""
import numpy as np, pandas as pd, itertools, math
P = [ # name, family, mass MeV, charge, color(1=colored), spin, weak isospin T3 (left-handed), generation
("electron","lepton",0.51099895,-1,0,0.5,-0.5,1),("muon","lepton",105.6583755,-1,0,0.5,-0.5,2),("tau","lepton",1776.93,-1,0,0.5,-0.5,3),
("nu_e","neutrino",0.0,0,0,0.5,0.5,1),("nu_mu","neutrino",0.0,0,0,0.5,0.5,2),("nu_tau","neutrino",0.0,0,0,0.5,0.5,3),
("up","quark",2.16,2/3,1,0.5,0.5,1),("down","quark",4.70,-1/3,1,0.5,-0.5,1),("strange","quark",93.5,-1/3,1,0.5,-0.5,2),
("charm","quark",1273.0,2/3,1,0.5,0.5,2),("bottom","quark",4183.0,-1/3,1,0.5,-0.5,3),("top","quark",172570.0,2/3,1,0.5,0.5,3),
("photon","gauge boson",0.0,0,0,1,0,0),("gluon","gauge boson",0.0,0,1,1,0,0),("W","gauge boson",80369.2,1,0,1,1,0),
("Z","gauge boson",91188.0,0,0,1,0,0),("Higgs","scalar",125200.0,0,0,0,-0.5,0),
("proton","baryon",938.27208816,1,0,0.5,0.5,0),("neutron","baryon",939.56542052,0,0,0.5,-0.5,0)]
cols=["name","family","mass_MeV","charge","colored","spin","T3","gen"]
df=pd.DataFrame(P,columns=cols)
CKM={("up","down"):0.97367,("up","strange"):0.22431,("up","bottom"):0.00382,("charm","down"):0.221,("charm","strange"):0.975,
     ("charm","bottom"):0.0408,("top","down"):0.0086,("top","strange"):0.0415,("top","bottom"):0.999}
lep={("electron","nu_e"),("muon","nu_mu"),("tau","nu_tau")}
v=246219.65  # Higgs vacuum value, MeV
def W_amp(a,b):
    k=(a,b) if (a,b) in CKM else (b,a)
    if k in CKM: return CKM[k]
    if (a,b) in lep or (b,a) in lep: return 1.0
    return 0.0
rows=[]
for i,j in itertools.combinations(range(len(df)),2):
    A,B=df.iloc[i],df.iloc[j]
    mi,mj=A.mass_MeV,B.mass_MeV
    r={"particle_i":A["name"],"particle_j":B["name"],"m_i":mi,"m_j":mj}
    if mi>0 and mj>0:
        r["mass_ratio"]=max(mi,mj)/min(mi,mj); r["M_ij"]=(mi+mj)/math.sqrt(mi*mj); r["Y_ij"]=math.log10(r["M_ij"]); r["S_ij"]=0.5*math.log(mi/mj)
    qq=A.charge*B.charge
    r["charge_product"]=round(qq,4)
    r["EM_push_pull"]="push (repel)" if qq>0 else ("pull (attract)" if qq<0 else "none")
    r["strong_link"]=int(A.colored and B.colored)
    r["W_transition_amplitude"]=W_amp(A["name"],B["name"])
    r["T3_product"]=A.T3*B.T3
    yi=math.sqrt(2)*mi/v if A.family not in("gauge boson",) else 0; yj=math.sqrt(2)*mj/v if B.family not in("gauge boson",) else 0
    r["Higgs_coupling_product"]=yi*yj
    r["same_generation"]=int(A.gen==B.gen and A.gen>0)
    rows.append(r)
pairs=pd.DataFrame(rows)
# Test: Cabibbo-type relations |V| vs sqrt(mass ratio) = exp(S)
tests=[]
for (u,d),V in CKM.items():
    if u=="up" and d=="down" or u=="charm" and d=="strange" or u=="top" and d=="bottom": continue
    m=dict(zip(df.name,df.mass_MeV))
    cands={"sqrt(m_lighter_down / m_heavier_down)":None}
    tests.append((u,d,V))
gst=[("V_us","sqrt(m_d/m_s)",CKM[("up","strange")],math.sqrt(4.70/93.5)),
     ("V_ub","sqrt(m_u/m_t)",CKM[("up","bottom")],math.sqrt(2.16/172570)),
     ("V_cb","sqrt(m_s/m_b)",CKM[("charm","bottom")],math.sqrt(93.5/4183)),
     ("V_cb","sqrt(m_c/m_t)",CKM[("charm","bottom")],math.sqrt(1273/172570)),
     ("V_td","sqrt(m_d/m_b)",CKM[("top","down")],math.sqrt(4.70/4183)),
     ("V_ts","sqrt(m_s/m_b)",CKM[("top","strange")],math.sqrt(93.5/4183))]
gt=pd.DataFrame(gst,columns=["CKM element","relational guess = exp(S_ij)","measured","predicted"])
gt["percent_off"]=(gt.predicted/gt.measured-1)*100
print(gt.round(4).to_string(index=False))
with pd.ExcelWriter("TRPG_Every_Particle_to_Every_Particle.xlsx") as w:
    df.to_excel(w,sheet_name="Particles",index=False); pairs.to_excel(w,sheet_name="All_Pairs",index=False); gt.to_excel(w,sheet_name="Mixing_vs_Mass_Test",index=False)
pairs.to_csv("TRPG_All_Pairs.csv",index=False)
print(len(pairs),"pairs")
# summary: how many pairs linked by each force
print("EM-linked:",(pairs.charge_product!=0).sum()," strong-linked:",pairs.strong_link.sum()," W-linked:",(pairs.W_transition_amplitude>0).sum())

# ---- Scale_Views (added 13 Nov 2026, Matthew and Melissa's idea: each particle sets its own scale view of the one universe) ----
# Constants: CODATA 2018; H0 = 67.4 km/s/Mpc (Planck 2018).
hbar=1.054571817e-34;c=2.99792458e8;G=6.67430e-11;H0=67.4e3/3.0857e22;MeVJ=1.602176634e-13
mPl=math.sqrt(hbar*c/G); RH=c/H0; Muni=c**3/(2*G*H0)
massive=df[df.mass_MeV>0].reset_index(drop=True)
view=pd.DataFrame({"reference_particle":massive.name})
for _,B in massive.iterrows(): view["%s_in_ref_units"%B["name"]]=massive.mass_MeV.map(lambda mr:B.mass_MeV/mr)
kg=massive.mass_MeV*MeVJ/c**2
view.insert(1,"Compton_length_fm",hbar/(kg*c)*1e15)
view.insert(2,"gravity_vs_quantum_1_over_alphaG",(mPl/kg)**2)
view.insert(3,"Hubble_radius_in_ref_lengths",RH/(hbar/(kg*c)))
view.insert(4,"universe_mass_in_ref_particles",Muni/kg)
with pd.ExcelWriter("TRPG_Every_Particle_to_Every_Particle.xlsx",mode="a",engine="openpyxl",if_sheet_exists="replace") as w:
    view.to_excel(w,sheet_name="Scale_Views",index=False)
print(view.iloc[:,:5].to_string())

"""TRPG and the cosmos (test T4), 14 Nov 2026. Planck 2018: H0=67.4, Omega_L=0.685, Omega_k=0.001+-0.002."""
import math
c=2.99792458e8; G=6.674e-11; hbarc_MeVfm=197.3269804; Mpc=3.0857e22; MeV_J=1.602176634e-13
H0=67.4e3/Mpc; rho_c=3*H0**2/(8*math.pi*G)              # kg/m^3
uL=0.685*rho_c*c**2                                      # J/m^3 dark energy density
uL_MeVfm3=uL/MeV_J*1e-45
B=66.0
print("dark energy density = %.3e MeV/fm^3 ; proton-scale baseline B = %.0f MeV/fm^3 ; ratio B/rho_L = %.1e"%(uL_MeVfm3,B,B/uL_MeVfm3))
# energy scales (fourth roots) in natural units: E = (u * (hbar c)^3)^(1/4)
EL=(uL_MeVfm3*hbarc_MeVfm**3)**0.25*1e9   # meV
EB=(B*hbarc_MeVfm**3)**0.25
print("baseline energy scale: proton layer %.0f MeV ; dark-energy layer %.2f meV"%(EB,EL))
print("compare M2 neutrino masses (meV): 0.36, 8.7, 50.1 ; sqrt(m1*m2) = %.2f meV ; m2/4 = %.2f"%(math.sqrt(0.36*8.66),8.66/4))
# Hubble-sphere compactness
R=c/H0; M=(4*math.pi/3)*R**3*rho_c
print("Hubble sphere: R = %.2e m, M = %.2e kg, compactness 2GM/(Rc^2) = %.4f (exactly 1 for a flat universe)"%(R,M,2*G*M/(R*c**2)))

"""M15: 'degrees of separation' for particles (14 Nov 2026).
1) Small-world estimate L ~ ln N / ln k (Watts & Strogatz 1998) for people and for the universe's particles.
2) Degrees of separation in Matthew and Melissa's own every-particle table (K5): links = any direct interaction."""
import math, pandas as pd, itertools
for name,N,k in [("people (8e9, ~150 close contacts; Dunbar)",8e9,150),("particles in visible universe (~1e80), 10 links each",1e80,10),("same, 1000 links each",1e80,1000)]:
    print("%-58s ~%.1f steps"%(name,math.log(N)/math.log(k)))
df=pd.read_excel("TRPG_Every_Particle_to_Every_Particle.xlsx",sheet_name="All_Pairs")
def graph(use_higgs):
    E={}
    for _,r in df.iterrows():
        linked = r.charge_product!=0 or r.strong_link==1 or r.W_transition_amplitude>0 or (use_higgs and r.Higgs_coupling_product>0)
        if linked: E.setdefault(r.particle_i,set()).add(r.particle_j); E.setdefault(r.particle_j,set()).add(r.particle_i)
    return E
names=sorted(set(df.particle_i)|set(df.particle_j))
for use_higgs in (False,True):
    E=graph(use_higgs); dist=[]; unreached=0
    for a in names:
        seen={a:0}; fr=[a]
        while fr:
            nx=[]
            for u in fr:
                for v in E.get(u,()):
                    if v not in seen: seen[v]=seen[u]+1; nx.append(v)
            fr=nx
        for b in names:
            if b>a: (dist.append(seen[b]) if b in seen else None)
            if b>a and b not in seen: unreached+=1
    print("K5 table, Higgs links %s: max steps %d, average %.2f, pairs never connected %d"%("on" if use_higgs else "off",max(dist),sum(dist)/len(dist),unreached))
# Add force-carrier couplings (photon<->charged, gluon<->quarks, W/Z<->all weak-charged, Higgs<->massive)
P=pd.read_excel("TRPG_Every_Particle_to_Every_Particle.xlsx",sheet_name="Particles").set_index("name")
E=graph(True)
def link(a,b): E.setdefault(a,set()).add(b); E.setdefault(b,set()).add(a)
for n,r in P.iterrows():
    if n in("photon","gluon","W","Z","Higgs"): continue
    if r.charge!=0: link("photon",n)
    if r.colored: link("gluon",n)
    if r.T3!=0 or r.charge!=0: link("Z",n); link("W",n)
    if r.mass_MeV>0: link("Higgs",n)
link("W","photon"); link("W","Z"); link("Higgs","W"); link("Higgs","Z")
dist=[]; un=0
for a in names:
    seen={a:0}; fr=[a]
    while fr:
        nx=[]
        for u in fr:
            for v in E.get(u,()):
                if v not in seen: seen[v]=seen[u]+1; nx.append(v)
        fr=nx
    for b in names:
        if b>a:
            if b in seen: dist.append(seen[b])
            else: un+=1
print("with force carriers: max steps %d, average %.2f, never connected %d (of 171)"%(max(dist),sum(dist)/len(dist),un))

"""M11 expansion from Matthew and Melissa's 'full bag / full glass' and M12 uncertainty as wave math. 14 Nov 2026.
Planck 2018 values. Friedmann (1922), Lemaitre (1927)."""
import numpy as np, math
from scipy.integrate import quad
H0=67.4; Om=0.315; Or=9.1e-5; OL=1-Om-Or   # full glass: parts sum to exactly 1
print("full glass: matter %.3f + dark energy %.3f + radiation %.5f = %.4f (Planck: curvature 0.001 +/- 0.002)"%(Om,OL,Or,Om+OL+Or))
H=lambda a: H0*math.sqrt(Om/a**3+Or/a**4+OL)
Mpc_km=3.0857e19; Gyr_s=3.156e16
age=quad(lambda a: 1/(a*H(a)),1e-10,1)[0]*Mpc_km/Gyr_s
print("age from the expansion equation: %.2f billion years (Planck 2018: 13.787 +/- 0.020)"%age)
for a in [1e-9,1e-4,0.5,1,2]:
    print("size %.0e of today: matter share %.3f, radiation %.3f, dark energy %.3f"%(a,(Om/a**3)/(H(a)/H0)**2,(Or/a**4)/(H(a)/H0)**2,OL/(H(a)/H0)**2))
# local conservation check: rho_m * a^3 constant (same amount of 'water'), rho_L constant per volume
print("matter amount in a stretching region (rho*a^3): constant -> same water; dark energy per volume constant -> its total grows: global total not conserved, only local (as in Einstein's theory)")
# M12: uncertainty as pure wave math (Fourier): Gaussian pulse gives dx*dk = 1/2 exactly
x=np.linspace(-50,50,2**16); s=2.0; psi=np.exp(-x**2/(4*s*s)); p=abs(psi)**2; p/=np.trapezoid(p,x)
dx=math.sqrt(np.trapezoid(x**2*p,x))
k=np.fft.fftshift(np.fft.fftfreq(len(x),x[1]-x[0]))*2*np.pi; P=abs(np.fft.fftshift(np.fft.fft(psi)))**2; P/=np.trapezoid(P,k)
dk=math.sqrt(np.trapezoid(k**2*P,k)); print("wave pulse: spread in place x spread in wavenumber = %.4f (Heisenberg minimum 0.5)"%(dx*dk))

"""Where the '2' in resistance n = 1 + 2GM/(r c^2) comes from (14 Nov 2026).
Light deflection at the Sun's edge = (1+gamma)/2 * 4GM/(R c^2); gamma = how much mass-energy stretches SPACE relative to how much it slows TIME (PPN; Eddington 1922, Will 1971)."""
import math
G,c,M,R=6.674e-11,2.998e8,1.989e30,6.957e8
full=4*G*M/(R*c*c)*206265
cases=[("scalar resistance (one number, couples to trace of energy; Nordstrom 1913)",None),
       ("time-only resistance (Newton/Soldner 1801)",0.0),
       ("time + space resistance, equal shares (tensor; Einstein 1915)",1.0)]
for name,g in cases:
    d=0.0 if g is None else (1+g)/2*full
    print("%-75s %.3f arcsec"%(name,d))
print("measured: 1.75 arcsec; Cassini (2003) gamma - 1 = (2.1 +/- 2.3) x 10^-5")
# split of the 2: time part and space part each contribute GM/(r c^2)
r=R; print("at the Sun's edge: time part %.2e + space part %.2e = %.2e = n - 1"%(G*M/(r*c*c),G*M/(r*c*c),2*G*M/(r*c*c)))

"""Properties 1 (long range) and 2 (acts on all energy equally) for TRPG's resistance field. 14 Nov 2026."""
hc_eVm=1.23984198e-6
for name,m_eV in [("containment field ripple P (M1, glueball)",1.72e9),("graviton upper limit (LIGO-Virgo-KAGRA GWTC-3, 2021)",1.27e-23)]:
    lam=hc_eVm/m_eV; print("%-55s mass %.2e eV -> reach ~ %.2e m"%(name,m_eV,lam))
print("reach of P = %.2f fm  -> P cannot be gravity (too short by > 10^30)"%(hc_eVm/1.72e9*1e15))
print("GW speed vs light (GW170817): equal to within ~3e-15 -> consistent with massless")
print("MICROSCOPE (2022): different materials fall alike to 1.5e-15 -> every TRPG layer must feel the same resistance to that level")

"""Tightening the gravity chain (14 Nov 2026).
A) Area law in Matthew and Melissa's crossing-wave medium: couplings vary as 1 + 0.5 cos(x) cos(y) (resistance from crossing waves).
B) Equal action on all layers: the area-law coefficient for fields of very different mass (different 'layers')."""
import numpy as np, math
def ent(N,m,crossing,L):
    idx=lambda i,j:i*N+j; K=np.zeros((N*N,N*N))
    for i in range(N):
        for j in range(N):
            a=idx(i,j)
            for di,dj in((1,0),(0,1)):
                ii,jj=i+di,j+dj
                if ii<N and jj<N:
                    xm,ym=(i+ii)/2,(j+jj)/2
                    k=1+0.5*math.cos(2*math.pi*xm/6)*math.cos(2*math.pi*ym/6) if crossing else 1.0
                    b=idx(ii,jj); K[a,b]=K[b,a]=-k; K[a,a]+=k; K[b,b]+=k
            K[a,a]+=m*m
    w,U=np.linalg.eigh(K); X=0.5*(U/np.sqrt(w))@U.T; P=0.5*(U*np.sqrt(w))@U.T
    out=[]
    for l in L:
        o=(N-l)//2; A=[idx(i,j) for i in range(o,o+l) for j in range(o,o+l)]
        nu=np.sqrt(np.clip(np.linalg.eigvals(X[np.ix_(A,A)]@P[np.ix_(A,A)]).real,0.25+1e-12,None))
        out.append(sum((n+.5)*math.log(n+.5)-(n-.5)*math.log(n-.5) for n in nu if n>0.5+1e-9))
    return np.array(out)
N=36; L=np.array([4,6,8,10,12,14])
S=ent(N,0.01,True,L)
eb=np.sqrt(np.mean((np.polyval(np.polyfit(4*L,S,1),4*L)-S)**2)); ev=np.sqrt(np.mean((np.polyval(np.polyfit(L*L,S,1),L*L)-S)**2))
print("A) crossing-wave medium: S =",np.round(S,3)," boundary-fit error %.3f vs volume-fit error %.3f (%.0fx better)"%(eb,ev,ev/eb))
print("B) area coefficient (entanglement per boundary unit) for different layer masses:")
for m in (0.01,0.1,0.3,1.0):
    S=ent(N,m,False,L); c=np.polyfit(4*L,S,1)[0]; print("   mass %.2f: coefficient %.4f"%(m,c))

"""TRPG mass-formula fix (13 Nov 2026).
1) Count Gemini's V51-V53 knobs. 2) Koide (1981) and Brannen (2006) as a zero/one-knob replacement in sqrt-mass space,
which is Matthew and Melissa's S_ij = 1/2 ln(mi/mj) space (sqrt(mi/mj) = e^S). 3) Random-mass test of how special the result is.
Masses: PDG 2024 (leptons pole; quarks MS-bar u,d,s @2GeV, c,b @own mass, t pole)."""
import numpy as np, math
me,mmu,mtau=0.51099895,105.6583755,1776.93
def Q(a,b,c): return (a+b+c)/(math.sqrt(a)+math.sqrt(b)+math.sqrt(c))**2
print("Gemini knobs: E0 (fit to proton), P0 (soliton), '/10', anchors 17 & 40, v=174.1 GeV, O(1) coeffs 0.218,0.415,0.439,1.71,2.66,0.9925 -> >=11 choices for 9 masses")
print("Koide Q leptons = %.7f (2/3 = %.7f)"%(Q(me,mmu,mtau),2/3))
# predict tau from e, mu with Q=2/3: solve for x=sqrt(mtau)
a,b=math.sqrt(me),math.sqrt(mmu); s=a+b
# (a^2+b^2+x^2) = 2/3 (s+x)^2 -> x^2/3 - 4/3 s x + (a^2+b^2) - 2/3 s^2 = 0
A_,B_,C_=1/3,-4/3*s,(a*a+b*b)-2/3*s*s
x=(-B_+math.sqrt(B_*B_-4*A_*C_))/(2*A_); print("tau predicted from e,mu: %.3f MeV vs %.2f measured (%.4f%%)"%(x*x,mtau,(x*x/mtau-1)*100))
# Brannen form sqrt(m_k)=M(1+sqrt2 cos(d+2pi k/3)); fit d from data, compare to 2/9
r=np.sqrt([me,mmu,mtau]); M=r.mean()
d=[math.acos((ri/M-1)/math.sqrt(2)) for ri in r]
print("Brannen scale M^2 = %.4f MeV; angle data:"%(M*M), [round(v,5) for v in d])
for dd in [2/9]:
    pred=[ (M*(1+math.sqrt(2)*math.cos(dd+2*math.pi*k/3)))**2 for k in range(3)]
    print("delta=2/9, one scale -> masses",[round(p,4) for p in sorted(pred)])
# Random test: log-uniform triples over 1e-1..1e4 MeV, chance |Q-2/3| <= |Q_lep-2/3|
rng=np.random.default_rng(1); N=2_000_000
m=10**rng.uniform(-1,4,(N,3)); q=m.sum(1)/np.sqrt(m).sum(1)**2
tol=abs(Q(me,mmu,mtau)-2/3); print("random chance of matching 2/3 this closely: %.2e"%np.mean(abs(q-2/3)<=tol))
quarks={"u":2.16,"d":4.70,"s":93.5,"c":1273,"b":4183,"t":172570}
for trip in [("u","c","t"),("d","s","b"),("s","c","b"),("c","b","t")]:
    print(trip,"Q=%.4f"%Q(*[quarks[k] for k in trip]))
# One-input version: fix delta=2/9, take scale only from the tau, predict e and mu
d=2/9; Mt=math.sqrt(mtau)/(1+math.sqrt(2)*math.cos(d))
pm=(Mt*(1+math.sqrt(2)*math.cos(d+4*math.pi/3)))**2; pe=(Mt*(1+math.sqrt(2)*math.cos(d+2*math.pi/3)))**2
print("from tau alone: e=%.5f (%.3f%%), mu=%.4f (%.4f%%)"%(pe,(pe/me-1)*100,pm,(pm/mmu-1)*100))
print("Brannen scale x3 = %.2f MeV vs proton 938.27 (%.2f%%)"%(3*Mt*Mt,(3*Mt*Mt/938.272-1)*100))

"""TRPG master action, step M1: pin the containment ('persistence'/space-resistance) potential V(P) with no tuning.
V(P) = (lam/4) v^4 (P^2-1)^2 ; P=1 empty space (full resistance), P=0 inside a proton (resistance removed).
Inputs: B from neutron-star + tidal window (K1,K2) and the lightest glueball mass from lattice QCD
(Morningstar & Peardon 1999; Chen et al. 2006: 0++ ~ 1.71-1.73 GeV). Model form: Friedberg & Lee (1977)."""
import math
hbarc=197.3269804; MeV4_per_MeVfm3=hbarc**3
for B in (57.1,60.1,63.7,62,70):   # proton-set band (proton_sets_B) and the old neutron-star window
    for mg in (1710.,1730.):
        Bn=B*MeV4_per_MeVfm3; lam=mg**4/(16*Bn); v=math.sqrt(mg**2/(2*lam))
        wall=hbarc/mg
        print(f"B={B} MeV/fm3  glueball={mg:.0f} MeV -> lam={lam:7.0f}  v={v:5.1f} MeV  wall thickness ~{wall:.3f} fm")
# Thin wall (0.11 fm) << core radius (~0.95 fm bag, 0.7 fm rms) -> MIT bag limit is valid, so K1-K3 carry over unchanged.
# Compare with the QCD vacuum energy from the gluon condensate (Shifman, Vainshtein & Zakharov 1979):
G2=0.012 # GeV^4, <alpha_s G^2/pi>
eps=(9/32)*G2*1e12/MeV4_per_MeVfm3
print(f"SVZ vacuum energy ~ {eps:.0f} MeV/fm3, about {eps/66:.1f}x the bag B (a known, unresolved gap)")

"""TRPG master action, step M2: the mass couplings Y as a relational 3-generation matrix (13 Nov 2026).
sqrt-mass matrix  S = M0 * [[1, r w, r w*],[r w*, 1, r w],[r w, r w*, 1]],  w = e^{i delta}  (cyclic / circulant)
eigenvalues M0(1 + 2r cos(delta + 2 pi k/3)); Koide Q=2/3  <=>  r = 1/sqrt2.
Refs: Koide 1981/1983; Foot 1994; Brannen 2006. Neutrino data: NuFIT 6.0 (2024). Cosmology: DESI DR2 (2025)."""
import numpy as np, math
me,mmu,mtau=0.51099895,105.6583755,1776.93
r=1/math.sqrt(2); d=2/9
M0=(math.sqrt(me)+math.sqrt(mmu)+math.sqrt(mtau))/3
w=np.exp(1j*d); S=M0*np.array([[1,r*w,r*np.conj(w)],[r*np.conj(w),1,r*w],[r*w,r*np.conj(w),1]])
ev=np.sort(np.linalg.eigvalsh(S)); print("matrix -> masses",np.round(ev**2,4),"measured",[me,mmu,mtau])
# 'equal self vs relational' check: |common part|^2 vs |relational part|^2 of the sqrt-mass vector
v=np.sqrt([me,mmu,mtau]); c=v.mean()*np.ones(3); p=v-c
print("self (common) part^2 = %.3f, relational part^2 = %.3f  (Koide = equal)"%(c@c,p@p))
# Neutrino test: impose Q=2/3 with measured splittings (normal ordering), allow one negative sqrt (as for Brannen)
dm21,dm31=7.49e-5,2.513e-3   # eV^2, NuFIT 6.0 normal ordering
from scipy.optimize import brentq
def Q(m,s): rt=np.array(s)*np.sqrt(m); return m.sum()/rt.sum()**2
for s in [(1,1,1),(-1,1,1)]:
    f=lambda m1: Q(np.array([m1,math.sqrt(m1**2+dm21),math.sqrt(m1**2+dm31)]),s)-2/3
    grid=np.logspace(-6,0,4000); vals=[f(x) for x in grid]
    roots=[brentq(f,grid[i],grid[i+1]) for i in range(len(grid)-1) if np.sign(vals[i])!=np.sign(vals[i+1])]
    for m1 in roots:
        ms=[m1,math.sqrt(m1**2+dm21),math.sqrt(m1**2+dm31)]
        print("signs",s,"-> m1,m2,m3 = %.5f %.5f %.5f eV, sum = %.4f eV"%(*ms,sum(ms)))
    if not roots: print("signs",s,"-> no solution")
# Inverted ordering check: m3 lightest, m1 = sqrt(m3^2+|dm32|+...), use dm32=-2.484e-3
dm32=-2.484e-3
for s in [(1,1,1),(1,1,-1),(1,-1,1),(-1,1,1)]:
    def f(m3,s=s):
        m2=math.sqrt(m3**2-dm32); m1=math.sqrt(m2**2-dm21); return Q(np.array([m1,m2,m3]),s)-2/3
    grid=np.logspace(-6,0,4000); vals=[f(x) for x in grid]
    roots=[brentq(f,grid[i],grid[i+1]) for i in range(len(grid)-1) if np.sign(vals[i])!=np.sign(vals[i+1])]
    print("inverted signs",s,"roots (m3 eV):",[round(x,5) for x in roots])

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
CONV=1.3234e-6; MSUN=1.4766
B=54.1; a=0.00042
def p_gem(e):
    if e<=4*B: return max(1e-4,0.12*e*(e/(4*B))**0.65)
    d=e-4*B; return d/3+a*d*d/(1+1.45*a*d)
egrid=np.logspace(-1,np.log10(5000),40000)
def make_inv(pf):
    pg=np.array([pf(e) for e in egrid]); pg=np.maximum.accumulate(pg+egrid*1e-12)
    return lambda p: np.interp(p,pg,egrid)
def tov(pf,einv,ec,psurf=1e-6):
    pc=pf(ec)
    def rhs(r,y):
        p,m=y
        if p<=psurf: return [0,0]
        e=einv(p)*CONV; pg=p*CONV
        return [-(e+pg)*(m+4*np.pi*r**3*pg)/(r*(r-2*m))/CONV, 4*np.pi*r*r*e]
    ev=lambda r,y: y[0]-psurf; ev.terminal=True
    r0=1e-4; s=solve_ivp(rhs,[r0,60],[pc,4/3*np.pi*r0**3*ec*CONV],events=ev,rtol=1e-8,atol=1e-12,max_step=0.02)
    return s.t[-1], s.y[1][-1]/MSUN
inv=make_inv(p_gem)
print("Gemini's own EOS, solved properly:")
rows=[]
for ec in [280,350,430,520,650,780,920,1080,1240,1380,1500,1650,2000,2500,3000]:
    R,M=tov(p_gem,inv,ec); rows.append((ec,M,R))
    cs2=(p_gem(ec+0.5)-p_gem(ec-0.5))
    print(f"eps_c={ec:5d}  M={M:.3f} Msun  R={R:.2f} km  cs2={cs2:.3f}")
# low-density part: p = 0.12 e (e/216)^0.65 at e=216 -> p=26 MeV/fm3 (nuclear saturation pressure ~ 0)
print("p at nuclear density e=150:",p_gem(150), " (real nuclear matter ~ 0-3 MeV/fm3)")


# ---- SLy reference check and bag-star scan ----
import numpy as np
# (uses tov(), p_gem, B, a defined above)
# low-density check for Gemini EOS
inv=make_inv(p_gem)
print("Gemini EOS, low central densities:")
for ec in [20,40,80,120,160,200,240]:
    R,M=tov(p_gem,inv,ec); print(f"  eps_c={ec} M={M:.3f} R={R:.1f}")
# SLy piecewise polytrope (Read et al. 2009), cgs, p in dyn/cm2 via K rho^G * c^2
c2=8.98755e20
crust=[(6.80110e-9,1.58425,2.44034e7),(1.06186e-6,1.28733,3.78358e11),(5.32697e1,0.62223,2.62780e12),(3.99874e-8,1.35692,None)]
G1,G2,G3=3.005,2.988,2.851; rho1,rho2=10**14.7,10**15.0; p1=10**34.384/c2
K1=p1/rho1**G1
K4c,G4c=crust[3][:2]; rho0=(K4c/K1)**(1/(G1-G4c))
pieces=[]; lo=0
for K,G,hi in crust[:3]: pieces.append([K,G,lo,hi]); lo=hi
pieces.append([K4c,G4c,lo,rho0]); pieces.append([K1,G1,rho0,rho1])
K2=K1*rho1**(G1-G2); pieces.append([K2,G2,rho1,rho2])
K3=K2*rho2**(G2-G3); pieces.append([K3,G3,rho2,1e17])
# energy-density constants a_i for continuity
a=0; prev=None
for pc in pieces:
    K,G,lo,hi=pc
    if prev is None: pc.append(0.0)
    else:
        Kp,Gp,_,_,ap=prev; e_lo=(1+ap)*lo+Kp/(Gp-1)*lo**Gp
        pc.append(e_lo/lo-1-K/(G-1)*lo**(G-1))
    prev=pc
MEVFM3=1.7827e12  # g/cm3 per MeV/fm3
rr=np.logspace(3,16.5,60000); P=[];E=[]
for r in rr:
    for K,G,lo,hi,ai in pieces:
        if lo<=r<hi: P.append(K*r**G); E.append((1+ai)*r+K/(G-1)*r**G); break
P=np.array(P)/MEVFM3; E=np.array(E)/MEVFM3
p_sly=lambda e: np.interp(e,E,P); inv_s=lambda p: np.interp(p,P,E)
print("\nReal nuclear EOS (SLy):  p at e=150:",p_sly(150))
best=[]
for ec in np.linspace(200,3000,57):
    R,M=tov(p_sly,inv_s,ec,psurf=1e-12); best.append((ec,M,R))
best=np.array(best); i=best[:,1].argmax()
print(f"SLy max mass {best[i,1]:.3f} Msun at R={best[i,2]:.2f} km (eps_c={best[i,0]:.0f})")
j=np.argmin(abs(best[:i,1]-1.4)); print(f"SLy 1.4 Msun star: R={best[j,2]:.2f} km (M={best[j,1]:.3f})")
import numpy as np

hc=197.327
for Bv in [54.1,57.0,60.0,70.0]:
  for av in [0.0,0.00042]:
    pf=lambda e,Bv=Bv,av=av: max((e-4*Bv)/3 + av*(e-4*Bv)**2/(1+1.45*av*(e-4*Bv)),0.0)
    einv=lambda p,pf=pf,Bv=Bv: brentq(lambda e: pf(e)-p, 4*Bv*0.999999, 1e5)
    out=[]
    for ec in np.linspace(4*Bv*1.01,3500,60):
        R,M=tov(pf,einv,ec,psurf=1e-7); out.append((ec,M,R))
    out=np.array(out); i=out[:,1].argmax(); j=np.argmin(abs(out[:i,1]-1.4))
    print(f"B={Bv:5.1f} (B^1/4={(Bv*hc**3)**.25:.1f} MeV) wall={av}: Mmax={out[i,1]:.2f} R_Mmax={out[i,2]:.2f}  R1.4={out[j,2]:.2f}")

# K2: real tidal Love number k2 and Lambda (Hinderer 2008 method), self-bound surface correction
# Validation: SLy gives Lambda(1.4)=299 (published ~300).
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
CONV=1.3234e-6; MSUN=1.4766
def k2_of(C,y):
    t1=8*C**5/5*(1-2*C)**2*(2+2*C*(y-1)-y)
    t2=2*C*(6-3*y+3*C*(5*y-8))+4*C**3*(13-11*y+C*(3*y-2)+2*C*C*(1+y))+3*(1-2*C)**2*(2-y+2*C*(y-1))*np.log(1-2*C)
    return t1/t2
def star(pf,dpde,einv,ec,es_surf=0.0,psurf=1e-9):
    pc=pf(ec)
    def rhs(r,Y):
        p,m,y=Y
        e=einv(max(p,psurf))*CONV; P=p*CONV
        f=1-2*m/r; cs2=dpde(einv(max(p,psurf)))
        dp=-(e+P)*(m+4*np.pi*r**3*P)/(r*r*f)
        Q=4*np.pi/f*(5*e+9*P+(e+P)/cs2)-6/(r*r*f)-(2*(m+4*np.pi*r**3*P)/(r*r*f))**2
        F=(1-4*np.pi*r*r*(e-P))/f
        dy=-(y*y+y*F+r*r*Q)/r
        return [dp/CONV,4*np.pi*r*r*e,dy]
    ev=lambda r,Y:Y[0]-psurf; ev.terminal=True
    r0=1e-4
    s=solve_ivp(rhs,[r0,60],[pc,4/3*np.pi*r0**3*ec*CONV,2.0],events=ev,rtol=1e-9,atol=1e-14,max_step=0.01)
    R=s.t[-1]; M=s.y[1][-1]; y=s.y[2][-1]
    y-=4*np.pi*R**3*es_surf*CONV/M   # surface density jump correction (self-bound star)
    C=M/R; k2=k2_of(C,y); Lam=2/3*k2/C**5
    return M/MSUN,R,k2,Lam
def seq(name,B,av):
    pf=lambda e: max((e-4*B)/3+av*(e-4*B)**2/(1+1.45*av*(e-4*B)),0.0)
    h=1e-3
    dpde=lambda e: (pf(e+h)-pf(max(e-h,4*B)))/(e+h-max(e-h,4*B))
    einv=lambda p: brentq(lambda e: pf(e)-p,4*B,1e5) if p>0 else 4*B
    rows=[star(pf,dpde,einv,ec,es_surf=4*B) for ec in np.geomspace(4*B*1.01,2500,70)]
    rows=np.array(rows); i=rows[:,0].argmax(); r=rows[:i+1]
    L14=float(np.exp(np.interp(1.4,r[:,0],np.log(r[:,3])))); k14=np.interp(1.4,r[:,0],r[:,2])  # log(Lambda) interpolated: Lambda ~ C^-5
    print(f"{name}: Mmax={rows[i,0]:.2f}  k2(1.4)={k14:.3f}  Lambda(1.4)={L14:.0f}")
    return rows
# Main run: tidal test (GW170817 needs Lambda(1.4) below about 580; heaviest star at least 2.0 Msun)
if __name__ == "__main__":
    for B, av in [(54.1, 0.0), (62.0, 0.0), (62.0, 0.0002), (66.0, 0.0002), (70.0, 0.0002), (70.0, 0.0003)]:
        seq(f"B={B} wall={av}", B, av)

"""Neutron stars with nothing adjusted (14 Nov 2026).

Inward pull B comes from the proton (proton_sets_B_2026-11-14.py: B = 60.1, band 57-64 MeV/fm^3).
Outward push ('wall') comes from the known strong-force strength, not a knob:
  quark matter pressure  p(mu) = (3/(4 pi^2)) [a4 mu^4 - m_s^2 mu^2] - B,   a4 = 1 - 2 alpha_s/pi
  (Freedman & McLerran 1977; Farhi & Jaffe 1984; Fraga, Pisarski & Schaffner-Bielich 2001; Alford, Braby, Paris & Reddy 2005)
  alpha_s from two-loop running of the measured alpha_s(MZ) = 0.1180 down to scale Lbar = X mu, X = 2 to 4
  (the usual renormalisation-scale uncertainty, reported as a band, not tuned). m_s = 100 MeV; no pairing gap.
Tests: heaviest star >= 2.0 Msun (PSR J0740+6620, J0348+0432); Lambda(1.4) < ~580 (GW170817).
"""
import numpy as np, math, importlib.util, os
from scipy.interpolate import interp1d
here = os.path.dirname(os.path.abspath(__file__))
spec = importlib.util.spec_from_file_location("tid", os.path.join(here, "neutron_star_tidal_fix_2026-11-13.py"))
tid = importlib.util.module_from_spec(spec); spec.loader.exec_module(tid)
hc = 197.3269

def alpha_s(Q):  # Q in GeV, two-loop, flavour thresholds at mb, mc
    def run(a, Q0, Q1, nf, n=400):
        b0 = (33-2*nf)/(12*math.pi); b1 = (153-19*nf)/(24*math.pi**2)
        t0, t1 = math.log(Q0**2), math.log(Q1**2); h = (t1-t0)/n
        f = lambda a: -(b0*a*a + b1*a**3)
        for _ in range(n):
            k1 = f(a); k2 = f(a+h*k1/2); k3 = f(a+h*k2/2); k4 = f(a+h*k3); a += h*(k1+2*k2+2*k3+k4)/6
        return a
    a = run(0.1180, 91.1876, 4.18, 5); a = run(a, 4.18, 1.27, 4); return run(a, 1.27, Q, 3)

def eos(B, X, ms=100.0, gap=0.0):
    mu = np.linspace(560.0/X, 1200, 3000)   # quark chemical potential, MeV
    a4 = np.array([1-2*alpha_s(X*m/1000)/math.pi for m in mu])
    p = (3/(4*math.pi**2))*(a4*mu**4 - (ms**2 - 4*gap**2)*mu**2)/hc**3 - B
    n = np.gradient(p, mu)               # quark density, fm^-3
    e = mu*n - p
    ok = p > 0
    i0 = np.argmax(ok)
    # surface: where p crosses 0
    es = np.interp(0.0, p[i0-1:i0+1], e[i0-1:i0+1])
    P, E = np.concatenate([[0.0], p[ok]]), np.concatenate([[es], e[ok]])
    return P, E, es, a4[i0]

def sequence(B, X, gap=0.0):
    P, E, es, a4s = eos(B, X, gap=gap)
    pf = interp1d(E, P, bounds_error=False, fill_value=(0.0, P[-1]))
    ef = interp1d(P, E, bounds_error=False, fill_value=(es, E[-1]))
    pf_ = lambda e: float(max(pf(e), 0.0))
    einv = lambda p: float(ef(p)) if p > 0 else es
    h = 0.5
    dpde = lambda e: max((pf_(e+h)-pf_(max(e-h, es)))/(e+h-max(e-h, es)), 1e-4)
    rows = np.array([tid.star(pf_, dpde, einv, ec, es_surf=es) for ec in np.geomspace(es*1.01, min(E[-1]*0.95, 3000), 70)])
    i = rows[:, 0].argmax(); r = rows[:i+1]
    L14 = float(np.exp(np.interp(1.4, r[:, 0], np.log(r[:, 3])))); R14 = np.interp(1.4, r[:, 0], r[:, 1])  # log(Lambda) interpolated: Lambda ~ C^-5
    return rows[i, 0], L14, R14, a4s, es

if __name__ == "__main__":
    print("alpha_s check: alpha_s(1.0 GeV) = %.3f, alpha_s(1.5 GeV) = %.3f" % (alpha_s(1.0), alpha_s(1.5)))
    print(" B (proton)  scale X   a4 at surface   Mmax (Msun)   R1.4 (km)   Lambda1.4   2.0 test   tidal test")
    for B in (57.1, 60.1, 63.7):
        for X in (2.0, 3.0, 4.0):
            M, L, R, a4, es = sequence(B, X)
            print(f"  {B:5.1f}      {X:3.1f}      {a4:.3f}          {M:.2f}         {R:5.2f}       {L:5.0f}      {'PASS' if M>=2.0 else 'fail'}       {'PASS' if L<580 else 'fail'}")

    print()
    print("What is missing: quark pairing (colour superconductivity; Alford, Rajagopal & Wilczek 1998), i.e. quarks")
    print("entangling in pairs at high density. Theory estimates the pairing gap at about 50-150 MeV; it is not measured.")
    print(" B      X    gap (MeV)   Mmax   R1.4   Lambda1.4   2.0 test   tidal test")
    for B in (57.1, 60.1, 63.7):
        for gap in (50.0, 100.0, 150.0):
            M, L, R, a4, es = sequence(B, 3.0, gap)
            print(f" {B:5.1f}  3.0   {gap:5.0f}      {M:.2f}   {R:5.2f}   {L:5.0f}      {'PASS' if M>=2.0 else 'fail'}       {'PASS' if L<580 else 'fail'}")

"""Re-run of Gemini V50/V51 proton 'predictions' using Gemini's own formulas and typed constants (13 Nov 2026)."""
import numpy as np
MN,MD,mpi=1085-3*48.91,1085+3*48.91,139.57      # Gemini: centroid 1085 & C_ss 48.91 typed -> split = 6*48.91 = measured gap
p=np.sqrt((MD**2-(MN+mpi)**2)*(MD**2-(MN-mpi)**2))/(2*MD); EN=np.sqrt(MN**2+p**2)
G=lambda f:(f**2/(12*np.pi*mpi**2))*((EN+MN)/(2*MD))*p**3
print("Delta width with Gemini f=2.18: %.1f MeV (Gemini claimed 114.2; measured 117). f needed for 117: %.2f"%(G(2.18),2.18*np.sqrt(117/G(2.18))))
Z2=1/1.214
print("gA from stated formula: %.4f (Gemini then adds +0.0517 by hand -> 1.2781)"%(Z2*(1.5129-(2/27)*0.214*1.5129)))
print("mu_p, mu_n from Z2*bare: %.3f, %.3f (Gemini typed 2.368, -1.518, then +/-0.412)"%(Z2*2.5424,Z2*-1.6949))
print("Delta m_np: uses PDG m_d-m_u as input and multiplies EM energy by a chosen 2.85")
print("von Laue: wall tension rescaled by scale_b=-int_q/int_b, so residual = 0 by construction")

"""Set the containment depth B from the proton and the Delta (14 Nov 2026), with no adjustable numbers.

Matthew and Melissa's containment balance: the bag radius R is where inward pull (B) balances outward quark push.
Model pieces (all from established physics, credited):
  - MIT bag: Chodos, Jaffe, Johnson, Thorn, Weisskopf (1974); DeGrand, Jaffe, Johnson, Kiskis (1975)
    E = (4pi/3) B R^3 + (3 x0 - Z0)/R + colour-magnetic term   (x0 = 2.0428 for massless quarks)
  - centre-of-mass correction M = sqrt(E^2 - <P^2>): Donoghue & Johnson (1980)
  - pion cloud ('leakage'): cloudy bag model, Theberge, Thomas & Miller (1980); Chew & Low static couplings
      Sigma_N = -(3 f^2/(4 pi^2 m^2)) [ I_N + (32/25) I_Delta ],  Sigma_Delta = -(3 f^2/(4 pi^2 m^2)) [ (8/25) PV I_N + I_Delta ]
      form factor u(kR) = 3 j1(kR)/(kR); f^2/4pi = 0.0778 (measured pion-nucleon coupling)
Fits: unknowns B and alpha_c to the measured N (938.9) and Delta (1232) masses; each radius from its own balance.
Z0 (the old free zero-point term) is set to 0 in the TRPG run: its job is done by the CM correction and the cloud.
"""
import numpy as np, math
from scipy.integrate import quad
from scipy.optimize import minimize_scalar, fsolve

hc = 197.3269; x0 = 2.04279; MPI = 138.0; DM = 293.1
f2 = 0.0778*4*math.pi
MN_obs, MD_obs = 938.92, 1232.0
MAG = 0.17650*8*0.25    # massless colour-magnetic integral x pair sum (N: -, Delta: +)

def u(kR): return 1.0 if kR < 1e-6 else 3*(math.sin(kR)-kR*math.cos(kR))/kR**3
def I(R, dE, pv=False):
    Rk = R/hc
    if not pv:
        return quad(lambda k: k**4*u(k*Rk)**2/(math.sqrt(k*k+MPI*MPI)*(math.sqrt(k*k+MPI*MPI)+dE)), 0, 4000, limit=400)[0]
    # principal value at omega = DM (Delta -> N pi is open)
    k0 = math.sqrt(DM**2-MPI**2)
    # 1/(w - DM) = (w + DM)/(k^2 - k0^2); write integrand as h(k)/(k-k0)
    h = lambda k: k**4*u(k*Rk)**2/math.sqrt(k*k+MPI*MPI)*(math.sqrt(k*k+MPI*MPI)+DM)/(k+k0)
    return quad(h, 0, 4000, weight='cauchy', wvar=k0, limit=400)[0]
def sigma(R, which="N"):
    c = -3*f2/(4*math.pi**2*MPI**2)
    return c*(I(R, 0) + 32/25*I(R, DM))

def sigmaD(R):
    c = -3*f2/(4*math.pi**2*MPI**2)
    return c*(8/25*I(R, -DM, pv=True) + I(R, 0))

def mass(R, B, ac, Z0, which, cm=True, cloud=True):
    sgn = -1 if which == "N" else +1
    E = 4*math.pi/3*B*R**3 + hc*(3*x0 - Z0)/R + sgn*MAG*ac*hc/R
    M = math.sqrt(max(E*E - (3*(x0*hc/R)**2 if cm else 0), 1.0))
    if cloud:
        M += sigma(R, "N") if which == "N" else sigmaD(R)
    return M
def balanced(B, ac, Z0, which, cm=True, cloud=True):
    # Containment balance at the wall: quark push against B sets R (the cloud sits outside the wall, as leakage,
    # and is added afterwards; cloudy bag model practice).
    r = minimize_scalar(lambda R: mass(R, B, ac, Z0, which, cm=cm, cloud=False), bounds=(0.3, 2.5), method="bounded", options={"xatol":1e-5})
    return mass(r.x, B, ac, Z0, which, cm=cm, cloud=cloud), r.x
def fit(Z0, **kw):
    def eqs(v):
        B, ac = abs(v[0]), abs(v[1])
        return [balanced(B, ac, Z0, "N", **kw)[0]-MN_obs, balanced(B, ac, Z0, "D", **kw)[0]-MD_obs]
    B, ac = fsolve(eqs, [57.0, 1.5])
    B, ac = abs(B), abs(ac)
    mN, RN = balanced(B, ac, Z0, "N", **kw); mD, RD = balanced(B, ac, Z0, "D", **kw)
    return B, ac, RN, RD, mN, mD
def core_rms(R):
    j0 = lambda x: math.sin(x)/x
    j1 = lambda x: math.sin(x)/x**2 - math.cos(x)/x
    rho = lambda r: (j0(x0*r/R)**2 + j1(x0*r/R)**2)*r*r
    return math.sqrt(quad(lambda r: rho(r)*r*r, 1e-9, R)[0]/quad(rho, 1e-9, R)[0])

if __name__ == "__main__":
    print("1) Check against the 1975 fit (Z0 = 1.84, no CM correction, no cloud): expect B^1/4 near 145 MeV, B near 57")
    B, ac, RN, RD, mN, mD = fit(1.84, cm=False, cloud=False)
    print(f"   B = {B:.1f} MeV/fm^3 (B^1/4 = {(B*hc**3)**0.25:.0f} MeV), alpha_c = {ac:.2f}, R_N = {RN:.2f} fm")
    print("2) TRPG run: Z0 = 0 (no free zero-point term), CM correction on, pion cloud (leakage) on")
    B, ac, RN, RD, mN, mD = fit(0.0)
    print(f"   B = {B:.1f} MeV/fm^3 (B^1/4 = {(B*hc**3)**0.25:.0f} MeV), alpha_c = {ac:.2f}")
    print(f"   proton bag radius R = {RN:.3f} fm, quark-core rms = {core_rms(RN):.3f} fm, cloud energy = {sigma(RN,'N'):.0f} MeV")
    print(f"   check: M_N = {mN:.1f}, M_Delta = {mD:.1f}")
    print("3) Sensitivity: pion-nucleon coupling +/- 5%")
    for s in (0.95, 1.05):
        f2 = 0.0778*4*math.pi*s
        Bs = fit(0.0)[0]
        print(f"   coupling x{s}: B = {Bs:.1f} MeV/fm^3")

"""Run all six quark masses to one common energy (2-loop QCD), then retest the Koide balance angle. 14 Nov 2026.
Inputs PDG 2024: alpha_s(MZ)=0.1180; u,d,s MS-bar at 2 GeV; c(mc)=1.273; b(mb)=4.183; t MS-bar(mt) ~ 162.5 GeV (from pole 172.57).
Running: d a/dln mu^2 = -b0 a^2 - b1 a^3 ; d ln m/d ln mu^2 = -g0 a - g1 a^2 ; a = alpha_s/pi (Tarrach; Chetyrkin; Vermaseren)."""
import math
from scipy.integrate import solve_ivp
MZ=91.1876
def b(nf): return (11-2*nf/3)/4, (102-38*nf/3)/16
def g(nf): return 1.0, (202/3-20*nf/9)/16
def nf_at(mu): return 6 if mu>162.5 else 5 if mu>4.183 else 4 if mu>1.273 else 3
def run(a0,mu0,mu1,lnm=0.0):
    # piecewise across thresholds, continuous matching (adequate at 2 loops for this test)
    th=[1.273,4.183,162.5]; pts=sorted([mu0,mu1]+[t for t in th if min(mu0,mu1)<t<max(mu0,mu1)],reverse=(mu1<mu0))
    a,l=a0,lnm
    for x,y in zip(pts,pts[1:]):
        nf=nf_at(math.sqrt(x*y)); b0,b1=b(nf); g0,g1=g(nf)
        f=lambda t,s:[-b0*s[0]**2-b1*s[0]**3, -g0*s[0]-g1*s[0]**2]
        s=solve_ivp(f,[math.log(x*x),math.log(y*y)],[a,l],rtol=1e-10,atol=1e-12).y[:,-1]; a,l=s
    return a,l
aZ=0.1180/math.pi
def a_at(mu): return run(aZ,MZ,mu)[0]
def m_to(m0,mu0,mu1): return m0*math.exp(run(a_at(mu0),mu0,mu1)[1])
inp={"u":(2.16e-3,2.0),"d":(4.70e-3,2.0),"s":(0.0935,2.0),"c":(1.273,1.273),"b":(4.183,4.183),"t":(162.5,162.5)}
print("alpha_s(2 GeV)=%.4f  alpha_s(mt)=%.4f"%(a_at(2.0)*math.pi,a_at(162.5)*math.pi))
def Qang(v):
    Q=sum(v)/sum(math.sqrt(x) for x in v)**2; return Q, math.degrees(math.acos(math.sqrt(1/(3*Q))))
for mu in (2.0,MZ,162.5,1000.0):
    m={k:m_to(m0,mu0,mu) for k,(m0,mu0) in inp.items()}
    up=Qang([m['u'],m['c'],m['t']]); dn=Qang([m['d'],m['s'],m['b']])
    print("at %7.1f GeV: u %.5f d %.5f s %.4f c %.3f b %.3f t %.1f GeV | up-type Q=%.4f angle %.2f | down-type Q=%.4f angle %.2f"%(mu,m['u'],m['d'],m['s'],m['c'],m['b'],m['t'],up[0],up[1],dn[0],dn[1]))

"""Quark sectors in the same relational (circulant) form as the leptons (14 Nov 2026).
sqrt-mass eigenvalues M0 (1 + 2 r cos(delta + 2 pi k/3)) give Koide Q = (1 + 2 r^2)/3, so r^2 = (3Q - 1)/2.
Leptons: Q = 2/3 -> r^2 = 1/2. Quarks: Q at a common scale (QCD running leaves Q unchanged; quark_common_scale script).
Uncertainties: PDG 2024 mass errors, Monte Carlo. Look-elsewhere: how often would simple sixths match this well by chance?"""
import numpy as np, math
rng = np.random.default_rng(0); N = 200000
# PDG 2024 MS-bar values (GeV) and 1-sigma errors, run to 2 GeV-equivalent ratios are scale-free for Q, so use one consistent set:
# values at mu = MZ from the quark_common_scale script (same running for all, so relative errors carry over)
mZ = {"u": (0.00128, 0.00128*0.07/2.16*1.0*10), "d": (0.00278, 0.00278*0.07/4.70*10),
      "s": (0.0552, 0.0552*0.008/0.0935*10/1.0), "c": (0.658, 0.658*0.0046/1.273), "b": (2.889, 2.889*0.007/4.183), "t": (170.3, 170.3*0.7/162.5)}
# relative errors (PDG 2024): u +0.49/-0.26 MeV of 2.16 (~17%), d +0.48/-0.17 of 4.70 (~7%), s +0.8/-0.3 of 93.5 (~6%),
# c 4.6/1273 (0.4%), b 7/4183 (0.2%), t(MS-bar) ~0.7/162.5 (0.4%)
rel = {"u": 0.17, "d": 0.07, "s": 0.06, "c": 0.0036, "b": 0.0017, "t": 0.0043}
cen = {k: v[0] for k, v in mZ.items()}
def Q(a, b, c): return (a+b+c)/(np.sqrt(a)+np.sqrt(b)+np.sqrt(c))**2
def r2(q): return (3*q-1)/2
def delta(ms):
    s = np.sqrt(ms); M = s.mean(); return math.acos((s[-1]/M-1)/(2*math.sqrt(r2(Q(*ms)))))
samp = {k: cen[k]*(1+rel[k]*rng.standard_normal(N)) for k in cen}
for k in samp: samp[k] = np.abs(samp[k])
qu = Q(samp["u"], samp["c"], samp["t"]); qd = Q(samp["d"], samp["s"], samp["b"])
for name, q, target, label in [("up-type (u,c,t)", qu, 5/6, "5/6"), ("down-type (d,s,b)", qd, 2/3, "2/3")]:
    rr = r2(q)
    print(f"{name}: Q = {np.median(q):.4f} +/- {q.std():.4f};  r^2 = {np.median(rr):.4f} +/- {rr.std():.4f};  "
          f"simple value {label} = {target:.4f} is {(np.median(rr)-target)/rr.std():+.1f} sigma away")
print("leptons: r^2 = 1/2 (Koide). Tested idea: r^2 = 3/6, 4/6, 5/6 for leptons, down-type, up-type -> FAILS (both quark sectors are ~7 sigma off).")
print("twist angles: leptons 2/9 = %.4f; up-type %.4f; down-type %.4f" % (2/9, delta(np.array([cen['u'], cen['c'], cen['t']])), delta(np.array([cen['d'], cen['s'], cen['b']]))))
# Look-elsewhere: chance that two random values (uniform on the allowed range r^2 in [0.5, 1)) both land within
# their measured distance of SOME sixth or ninth or twelfth (simple fractions a researcher might try)
fr = sorted({a/b for b in (2, 3, 4, 6, 8, 9, 12) for a in range(1, b) if 0.5 <= a/b < 1})
du = min(abs(np.median(r2(qu))-f) for f in fr); dd = min(abs(np.median(r2(qd))-f) for f in fr)
x = rng.uniform(0.5, 1.0, (N, 2)); near = lambda v, d: np.min(np.abs(v[:, None]-np.array(fr)[None, :]), axis=1) <= d
print(f"look-elsewhere: {len(fr)} simple fractions tried; chance both random values land as close: {np.mean(near(x[:,0],du)&near(x[:,1],dd)):.3f}")

"""Quark masses from linked Koide triples ('waterfall', Rivero 2011; cf. Rodejohann & Zhang 2011). 14 Nov 2026.
Each overlapping triple (with one negative sqrt where needed) is set to Q=2/3; start from top and bottom only.
Caveat: quark masses depend on scheme/scale (MS-bar), so % comparisons are approximate."""
import math
from scipy.optimize import brentq
m={"u":2.16,"d":4.70,"s":93.5,"c":1273,"b":4183,"t":172570}
def Q(vals,signs): r=[s*math.sqrt(v) for v,s in zip(vals,signs)]; return sum(vals)/sum(r)**2
print("Q(c,b,t)=%.4f  Q(-s,c,b)=%.4f  Q(d,s,-c)? see chain"%(Q([m['c'],m['b'],m['t']],[1,1,1]),Q([m['s'],m['c'],m['b']],[-1,1,1])))
def solve_light(heavy2,heavy1,sign_light):
    f=lambda x: Q([x,heavy2,heavy1],[sign_light,1,1])-2/3
    xs=[10**(k/200) for k in range(-600,1200)]; out=[]
    for a,b in zip(xs,xs[1:]):
        if f(a)*f(b)<0: out.append(brentq(f,a,b))
    return out
c=solve_light(m['b'],m['t'],1); print("charm from (b,t):",[round(v,1) for v in c],"measured",m['c'])
cc=min(c,key=lambda v:abs(v-m['c']))
s=solve_light(cc,m['b'],-1); print("strange from (c,b), negative root:",[round(v,1) for v in s],"measured",m['s'])

"""Matthew and Melissa (14 Nov): the wave pattern is energy spreading over regions of different resistance. Three checks.
1) Resistance as refractive index n = 1 + 2GM/(r c^2) around the Sun -> light bending & Shapiro delay (Eddington; Gordon 1923 optical metric).
2) Periodic resistance -> allowed/forbidden bands (Bloch 1928; Kronig-Penney 1931), via 1D transfer matrices.
3) Random resistance -> waves trapped in place (Anderson 1958 localization)."""
import numpy as np, math
G=6.674e-11; c=2.998e8; M=1.989e30; R=6.957e8
# 1) ray through n(r): deflection = integral of grad_perp n along straight path (weak field)
b=R; z=np.linspace(-2e12,2e12,2_000_001); r=np.sqrt(b*b+z*z)
dndb=-2*G*M/(c*c)*b/r**3
defl=-np.trapezoid(dndb,z); print("1) light bending at Sun's edge from resistance n(r): %.3f arcsec (measured 1.75; Einstein 1.75; Newton-only 0.87)"%(defl*206265))
au=1.496e11; rE=au; rM=1.52*au
shap=(2*G*M/c**3)*math.log(4*rE*rM/b**2); print("   extra delay radar Earth-Mars grazing Sun (one way, Shapiro): %.0f microseconds (measured ~ 100-250 us round trip range)"%(shap*1e6))
# 2) periodic resistance: alternating layers of index n1, n2; transfer matrix trace gives bands |Tr|/2<=1 allowed
def tr(k,n1=1.0,n2=1.5,a=0.5,bb=0.5):
    def layer(n,d):
        q=k*n; return np.array([[math.cos(q*d),math.sin(q*d)/n],[-n*math.sin(q*d),math.cos(q*d)]])
    T=layer(n1,a)@layer(n2,bb); return abs(np.trace(T))/2
ks=np.linspace(0.01,12,4000); allowed=np.array([tr(k)<=1 for k in ks])
edges=np.where(np.diff(allowed.astype(int))!=0)[0]
gaps=[(round(ks[edges[i]],2),round(ks[edges[i+1]],2)) for i in range(0,len(edges)-1,2) if not allowed[edges[i]+1]]
print("2) periodic resistance -> forbidden frequency gaps (bands between them):",gaps[:5])
# 3) random resistance: transmission through N random layers decays exponentially (localization)
rng=np.random.default_rng(0)
def T_random(N,k=3.0,spread=0.3):
    T=np.eye(2)
    for _ in range(N):
        n=1+spread*rng.random(); d=0.5; q=k*n
        T=np.array([[math.cos(q*d),math.sin(q*d)/n],[-n*math.sin(q*d),math.cos(q*d)]])@T
    t=2/abs(T[0,0]+T[1,1]+1j*(T[0,1]*1-T[1,0]/1)); return t*t
for N in [10,100,1000,3000]:
    print("3) %4d random layers: transmitted fraction %.2e"%(N,np.mean([T_random(N) for _ in range(20)])))

"""Black-hole ringdown and echo test for TRPG's 'central particle' black hole (13 Nov 2026).
Kerr l=m=2 fundamental mode fits: Berti, Cardoso & Will (2006). Echo delay: Abedi, Dykaar & Afshordi (2017).
GW150914 remnant: source mass 62 Msun, spin 0.67, z = 0.09 (LIGO/Virgo 2016)."""
import math
Msun_s=4.925491e-6; Msun_m=1476.625
M=62*(1+0.09); a=0.67
Mw=1.5251-1.1568*(1-a)**0.1292; Q=0.7000+1.4187*(1-a)**(-0.4990)
w=Mw/(M*Msun_s); f=w/(2*math.pi); tau=2*Q/w
print("GR (Kerr) ringdown: f = %.0f Hz, tau = %.1f ms   | measured 251 +/- 8 Hz, 4.0 +/- 0.3 ms"%(f,tau*1e3))
lP=1.616255e-35
dt=4*M*Msun_s*(1+1/math.sqrt(1-a*a))*math.log(M*Msun_m/lP)
print("If a reflecting central-particle surface sat a Planck length above the horizon: echoes every %.3f s"%dt)
# Compactness ladder C = 2GM/(Rc^2)
for n,C in [("Sun",4.24e-6),("white dwarf (Sirius B)",2*1.02*Msun_m/5.8e6),("neutron star 1.4 Msun, R 11.4 km (K1)",2*1.4*Msun_m/11.4e3),
            ("Buchdahl limit, any static star",8/9),("black hole horizon",1.0)]:
    print("%-40s C = %.3g"%(n,C))

"""The 2/9 twist and 1/sqrt2 link in the lepton matrix (M2). 14 Nov 2026. PDG 2024 masses and errors."""
import math, numpy as np
def delta(me,mmu,mtau):
    r=np.sqrt([me,mmu,mtau]); M=r.mean(); return math.acos((r[2]/M-1)/math.sqrt(2))
d0=delta(0.51099895,105.6583755,1776.93)
rng=np.random.default_rng(0)
ds=[delta(0.51099895,105.6583755+rng.normal(0,2.3e-6),1776.93+rng.normal(0,0.09)) for _ in range(20000)]
print("twist from data: %.6f +/- %.6f rad ; 2/9 = %.6f ; difference = %.1f sigma"%(d0,np.std(ds),2/9,(d0-2/9)/np.std(ds)))
Q=2/3; print("2/9 = Q/3 = (Koide ratio)/(number of generations) = %.6f"%(Q/3))
print("link 1/sqrt2: Koide Q=2/3 <=> shared part^2 = relational part^2 (shown in M2); r=1/sqrt2 is that balance")

"""Matthew and Melissa's crossing-wave picture of space: z = cos x + cos y (and its time dependence). 14 Nov 2026.
Checks: critical points (max / saddle / min), Gaussian curvature sign map, Gauss-Bonnet total, and standing-wave time flip."""
import numpy as np
n=1201; x=np.linspace(0,2*np.pi,n); X,Y=np.meshgrid(x,x)
zx,zy=-np.sin(X),-np.sin(Y); zxx,zyy=-np.cos(X),-np.cos(Y); zxy=0
K=(zxx*zyy-zxy**2)/(1+zx**2+zy**2)**2
dA=np.sqrt(1+zx**2+zy**2)*(x[1]-x[0])**2
print("critical points per cell: (0,0) peak z=2, (pi,pi) trough z=-2, (pi,0) & (0,pi) saddles z=0 -> 1 peak, 1 trough, 2 saddles")
print("Morse check (peaks - saddles + troughs) = 1 - 2 + 1 = 0 = Euler number of a repeating grid")
pos=(K>0); print("area fraction curved like a sphere (K>0): %.3f ; like a saddle (K<0): %.3f"%((dA*pos).sum()/dA.sum(),(dA*(K<0)).sum()/dA.sum()))
print("total curvature over one cell (Gauss-Bonnet says 0): %.2e ; peak |K| = %.3f"%((K*dA)[:-1,:-1].sum(),abs(K).max()))
# standing waves: z = (cos x + cos y) cos(wt): every half period peaks and troughs swap; saddles stay saddles but flip orientation; at wt=pi/2 flat
for wt in [0,np.pi/4,np.pi/2,np.pi]:
    a=np.cos(wt); print("wt=%.2f  height at (0,0)=%+.2f  curvature at (0,0)=%+.3f  at saddle (pi,0)=%+.3f"%(wt,2*a,a*a,-a*a))

# build_math_page.py - writes build/math.html ("The math, step by step").
#
# NOTHING on that page is typed.  Every constant and every measured reference value is read from math_inputs.json (each with its
# source), every derived value is computed here (logic ported from downloads/code/*.py), and every PASS / FAIL / PARTIAL badge is
# decided by code from a tolerance that is stated on the page.  After writing, a provenance check scans the page and refuses
# (non-zero exit) if any number is not one that this script emitted or part of an equation's own formula.
#
#   python -I build_math_page.py <build_dir>              # build + provenance check
#   python -I build_math_page.py <build_dir> --plant      # also type a number into a copy and prove the check catches it
import sys
sys.dont_write_bytecode = True
import html, importlib.util, json, math, pathlib, re, zipfile
import xml.etree.ElementTree as ET
import numpy as np
from scipy.integrate import quad, solve_ivp
from scipy.optimize import brentq

HERE = pathlib.Path(__file__).resolve().parent
spec = importlib.util.spec_from_file_location("site_common", HERE / "site_common.py")
sc = importlib.util.module_from_spec(spec); spec.loader.exec_module(sc)
BUILD = pathlib.Path(sys.argv[1])
INP = json.loads((HERE / "math_inputs.json").read_text(encoding="utf-8"))
def _walk(o):
    if isinstance(o, str):
        yield o
    elif isinstance(o, dict):
        for v in o.values(): yield from _walk(v)
    elif isinstance(o, list):
        for v in o: yield from _walk(v)
LOADED_TEXT = " ".join(_walk(INP))
def V(k): return INP[k]["value"]
def CL(k): return INP["record_claims"]["value"][k]
STDM = V("std_ns_mass")
def SRC(k): return INP[k]["source"]

# ------------------------------------------------------------------------------------------------ static check of THIS file
import ast
def static_check():
    """No string literal that reaches the page may contain a typed number, and no N()/I()/fr() call may take a numeric literal,
    except: equation strings (the 5th argument of check()), the ANCH anchor table, and the module docstring."""
    src = pathlib.Path(__file__).read_text(encoding="utf-8"); tree = ast.parse(src); bad = []
    exempt = set()
    for node in ast.walk(tree):
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") == "check" and len(node.args) > 4 and isinstance(node.args[4], ast.List):
            for d in ast.walk(node.args[4]): exempt.add(id(d))
        if isinstance(node, ast.Assign) and any(getattr(tg, "id", "") in ("ANCH", "NUM", "SUP") for tg in node.targets):
            for d in ast.walk(node): exempt.add(id(d))
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") == "eq":
            for d in ast.walk(node): exempt.add(id(d))
    for node in ast.walk(tree):
        if id(node) in exempt: continue
        if isinstance(node, ast.Constant) and isinstance(node.value, str) and typed_numbers(node.value):
            if node.value.startswith("# build_math_page") : continue
            bad.append((node.lineno, node.value[:60]))
        if isinstance(node, ast.Call) and getattr(node.func, "id", "") in ("N", "I") and node.args and isinstance(node.args[0], ast.Constant) and isinstance(node.args[0].value, (int, float)):
            bad.append((node.lineno, "numeric literal passed to " + node.func.id))
    return bad
def typed_numbers(s):
    """decimals and numbers of three or more digits typed into a string, unless the loaded data file itself contains them."""
    loaded = set(NUM.findall(LOADED_TEXT))
    return [tk for tk in NUM_S.findall(s) if ("." in tk or len(tk) >= 3) and tk not in loaded]
NUM_S = re.compile(r"(?<![A-Za-z_\d.\-])\d+(?:\.\d+)?")

# ------------------------------------------------------------------------------------------------ provenance machinery
NUM = re.compile(r"(?<![A-Za-z_\d.])\d+(?:\.\d+)?")
EMITTED, FORMULA = set(NUM.findall(LOADED_TEXT)), set()
def emit(s):
    EMITTED.update(NUM.findall(s)); return s
def eq(s):
    FORMULA.update(NUM.findall(s)); return s
SUP = str.maketrans("0123456789-", "⁰¹²³⁴⁵⁶⁷⁸⁹⁻")
def N(x, sig=4):
    """format a computed or loaded number; every digit group goes into the emitted set."""
    x = float(x)
    if x == 0: return emit("0")
    sign = "−" if x < 0 else ""; a = abs(x)
    if 1e-3 <= a < 1e6:
        dec = max(0, sig - 1 - int(math.floor(math.log10(a))))
        s = f"{a:.{dec}f}"
        if "." in s: s = s.rstrip("0").rstrip(".")
    else:
        e = int(math.floor(math.log10(a))); m = a / 10 ** e
        if round(m, sig - 1) >= 10: m /= 10; e += 1
        ms = f"{m:.{sig - 1}f}"
        if "." in ms: ms = ms.rstrip("0").rstrip(".")
        s = ms + "×10" + str(e).translate(SUP)
    return emit(sign + s)
def I(x): return emit(str(int(x)))
def fr(a, b): return emit(f"{int(a)}/{int(b)}")

if "--static" in sys.argv:
    b = static_check()
    for ln, s in b: print(ln, ascii(s))
    print("static offenders:", len(b)); sys.exit(1 if b else 0)
STATUSES = ["PASS", "FAIL", "PARTIAL", "NOT TESTABLE YET", "INPUT", "INHERITED", "NOT RE-RUN"]
CLS = {"PASS": "pass", "FAIL": "fail", "PARTIAL": "partial", "NOT TESTABLE YET": "nty", "INPUT": "input", "INHERITED": "inherited", "NOT RE-RUN": "norerun"}
PF = V("partial_factor")
def j_rel(x, ref, tol):
    d = abs(x - ref) / abs(ref)
    return "PASS" if d <= tol else ("PARTIAL" if d <= PF * tol else "FAIL")
def j_sigma(x, ref, sig, n=1):
    d = abs(x - ref) / sig
    return "PASS" if d <= n else ("PARTIAL" if d <= PF * n else "FAIL")
def j_range(x, lo, hi): return "PASS" if lo <= x <= hi else "FAIL"
MISMATCH = []
def rec_line(label, computed, claimed, sig=4):
    """the record's own stated number beside this run's number; a >1% difference is flagged and counted."""
    d = abs(computed - claimed) / abs(claimed) if claimed else abs(computed)
    ok = d <= V("tol_percent_1") or round(computed, 3) == claimed
    if not ok: MISMATCH.append(label)
    return (f"The record states {N(claimed, sig)}; this run gives {N(computed, sig)}. " + ("They agree." if ok else "They DIFFER."))

# ------------------------------------------------------------------------------------------------ record anchors
R = "record.html#"
ANCH = {
 "A": "item-a1", "K1": "item-k1", "K2": "item-k2", "K3": "item-k3", "K4": "item-k4", "K5": "item-k5",
 "K6": "k6.-each-particles-own-scale-view-matthew-13-nov", "K7": "k7.-the-particle-mass-formula-fixed-for-the-leptons",
 "K8": "k8.-geminis-other-proton-results-checked-line-by-line", "K9": "k9.-black-hole-ringdown-and-the-central-particle",
 "K10": "k10.-antiparticle-matrices-and-decay-rates-matthews-saved-requests-workbook-updated",
 "M1": "m1.-the-containment-potential-vp-pinned-with-no-knobs",
 "M2": "m2.-the-mass-couplings-y-written-as-relations-between-generations",
 "M3": "m3.-the-layer-structure-g_ab-matthews-s-is-the-only-fair-distance",
 "M5": "m5.-space-as-crossing-waves-why-the-geometry-keeps-transitioning-matthew-original-idea-newly-written-down-14-nov",
 "M6": "m6.-quarks-a-hint-not-yet-a-result",
 "M7": "m7.-the-cause-of-the-waves-energy-spreading-through-uneven-resistance-matthew-14-nov",
 "M8": "m8.-the-cosmos-test-t4-first-pass", "M9": "m9.-the-factor-2-what-kind-of-resistance-gravity-must-be",
 "M10": "m10.-gravitys-other-two-properties-reaching-everywhere-and-acting-on-everything-equally",
 "M11": "m11.-the-expansion-equation-the-full-bag-and-the-full-glass-matthew-original-theory",
 "M12": "m12.-entanglement-and-uncertainty-inside-the-bag",
 "M14": "m14.-gravity-as-a-byproduct-of-entanglement-matthew-14-nov-gravity-is-accidental-as-a-result-of-entanglement",
 "M15": "m15.-quantum-physics-is-entanglement-the-universe-as-one-connected-web-matthew-from-his-original-meditations",
 "M16": "m16.-entangled-particles-are-touching-the-area-law-charges-from-consistency",
 "M17": "m17.-final-scorecard-14-nov-2026", "M18": "m18.-the-last-two-open-problems-honest-best-attempts",
 "P1": "p1.-the-master-equation-written-out-in-full",
 "Q3": "q3.-b-set-by-the-proton-with-nothing-adjusted-script-math_fixesproton_sets_b_2026-11-14.py",
 "Q4": "q4.-neutron-stars-with-nothing-adjusted-script-math_fixesneutron_stars_no_knobs_2026-11-14.py",
 "Q5": "q5.-quark-masses-one-idea-tested-and-it-fails-script-math_fixesquark_link_pattern_2026-11-14.py",
 "P2": "p2.-quark-masses-at-one-common-energy-script-math_fixesquark_common_scale_2026-11-14.py",
 "P4": "p4.-gravity-tightened-script-math_fixesgravity_tighten_2026-11-14.py",
 "P5": "p5.-the-29-twist-and-the-12-link-script-math_fixestwist_and_link_2026-11-14.py",
}

CHECKS = []
SECTIONS = {}
def check(cid, sec, title, rec, eqs, subst, computed, measured, tol, status, note="", pre=None):
    if pre: emit(pre)
    assert status in STATUSES, status
    CHECKS.append(dict(id=cid, sec=sec, title=title, rec=rec, eqs=[eq(e) for e in eqs], subst=subst, computed=computed,
                       measured=measured, tol=tol, status=status, note=note, pre=pre))

# ------------------------------------------------------------------------------------------------ computations
G, C = V("G_SI"), V("c_SI"); MSUN, RSUN = V("Msun_kg"), V("Rsun_m"); ASEC = V("arcsec_per_rad")
me, mmu, mtau = V("m_e"), V("m_mu"), V("m_tau")
Q23 = V("koide_target_num") / V("koide_target_den"); D29 = V("brannen_delta_num") / V("brannen_delta_den")
def Qk(a, b, c): return (a + b + c) / (math.sqrt(a) + math.sqrt(b) + math.sqrt(c)) ** 2

def bending():
    b = RSUN; z = np.linspace(-2e12, 2e12, 2_000_001); r = np.sqrt(b * b + z * z)
    dndb = -2 * G * MSUN / (C * C) * b / r ** 3
    return -np.trapezoid(dndb, z) * ASEC

def age_gyr():
    Om, Or = V("Omega_m"), V("Omega_r"); OL = 1 - Om - Or
    H = lambda a: V("H0") * math.sqrt(Om / a ** 3 + Or / a ** 4 + OL)
    return quad(lambda a: 1 / (a * H(a)), 1e-10, 1)[0] * V("Mpc_km") / V("Gyr_s"), OL, H

# --- TOV (units: p, e in MeV/fm^3, r and m in km) ---
CONV, MSUN_KM = V("tov_conv"), V("Msun_geo_km")
def tov(pf, einv, ec, psurf=1e-6):
    pc = pf(ec)
    def rhs(r, y):
        p, m = y
        if p <= psurf: return [0, 0]
        e = einv(p) * CONV; pg = p * CONV
        return [-(e + pg) * (m + 4 * np.pi * r ** 3 * pg) / (r * (r - 2 * m)) / CONV, 4 * np.pi * r * r * e]
    ev = lambda r, y: y[0] - psurf; ev.terminal = True
    r0 = 1e-4
    s = solve_ivp(rhs, [r0, 60], [pc, 4 / 3 * np.pi * r0 ** 3 * ec * CONV], events=ev, rtol=1e-8, atol=1e-12, max_step=0.02)
    return s.t[-1], s.y[1][-1] / MSUN_KM

def bag_eos(B, av):
    pf = lambda e: max((e - 4 * B) / 3 + av * (e - 4 * B) ** 2 / (1 + 1.45 * av * (e - 4 * B)), 0.0)
    def einv(p):
        if p <= 0: return 4 * B
        if av == 0: return 3 * p + 4 * B
        A = 1.45 * av / 3 + av; Bc = 1 / 3 - 1.45 * av * p
        return 4 * B + (-Bc + math.sqrt(Bc * Bc + 4 * A * p)) / (2 * A)
    return pf, einv

def bag_sequence(B, av):
    pf, einv = bag_eos(B, av)
    rows = np.array([(ec,) + tov(pf, einv, ec, psurf=1e-7) for ec in np.linspace(4 * B * 1.01, V("tov_scan_max"), 60)])
    i = rows[:, 2].argmax(); r = rows[:i + 1]
    return float(rows[i, 2]), float(np.interp(STDM, r[:, 2], r[:, 1])), float(rows[i, 1])

def gemini_sequence():
    B, a = V("bag_B_gemini"), V("bag_wall_gemini"); sc_, pw = V("gemini_p_scale"), V("gemini_p_power")
    def pf(e):
        if e <= 4 * B: return max(1e-4, sc_ * e * (e / (4 * B)) ** pw)
        d = e - 4 * B; return d / 3 + a * d * d / (1 + 1.45 * a * d)
    egrid = np.logspace(-1, np.log10(5000), 40000)
    pg = np.maximum.accumulate(np.array([pf(e) for e in egrid]) + egrid * 1e-12)
    einv = lambda p: np.interp(p, pg, egrid)
    rows = np.array([(ec,) + tov(pf, einv, ec) for ec in sorted(V("gemini_ns_ec_list_low") + V("gemini_ns_ec_list"))])
    i = rows[:, 2].argmax(); r = rows[:i + 1]
    return float(rows[i, 2]), float(np.interp(STDM, r[:, 2], r[:, 1])), pf(150.0)

def sly_tables():
    c2 = V("c2_cgs"); crust = V("sly_crust"); G1, G2, G3 = V("sly_G"); lr1, lr2 = V("sly_log_rho12")
    rho1, rho2 = 10 ** lr1, 10 ** lr2; p1 = 10 ** V("sly_log_p1") / c2; K1 = p1 / rho1 ** G1
    K4c, G4c = crust[3][:2]; rho0 = (K4c / K1) ** (1 / (G1 - G4c))
    pieces = []; lo = 0
    for K, Gm, hi in crust[:3]: pieces.append([K, Gm, lo, hi]); lo = hi
    pieces.append([K4c, G4c, lo, rho0]); pieces.append([K1, G1, rho0, rho1])
    K2 = K1 * rho1 ** (G1 - G2); pieces.append([K2, G2, rho1, rho2])
    K3 = K2 * rho2 ** (G2 - G3); pieces.append([K3, G3, rho2, 1e17])
    prev = None
    for pc in pieces:
        K, Gm, lo, hi = pc
        if prev is None: pc.append(0.0)
        else:
            Kp, Gp, _, _, ap = prev; e_lo = (1 + ap) * lo + Kp / (Gp - 1) * lo ** Gp
            pc.append(e_lo / lo - 1 - K / (Gm - 1) * lo ** (Gm - 1))
        prev = pc
    conv = V("mevfm3_gcm3"); rr = np.logspace(3, 16.5, 60000); P = []; E = []
    for r in rr:
        for K, Gm, lo, hi, ai in pieces:
            if lo <= r < hi: P.append(K * r ** Gm); E.append((1 + ai) * r + K / (Gm - 1) * r ** Gm); break
    return np.array(E) / conv, np.array(P) / conv   # note: arrays are built so that (E, P) rise together

def k2_of(Cc, y):
    t1 = 8 * Cc ** 5 / 5 * (1 - 2 * Cc) ** 2 * (2 + 2 * Cc * (y - 1) - y)
    t2 = (2 * Cc * (6 - 3 * y + 3 * Cc * (5 * y - 8)) + 4 * Cc ** 3 * (13 - 11 * y + Cc * (3 * y - 2) + 2 * Cc * Cc * (1 + y))
          + 3 * (1 - 2 * Cc) ** 2 * (2 - y + 2 * Cc * (y - 1)) * np.log(1 - 2 * Cc))
    return t1 / t2

def tidal_star(pf, einv, dpde, ec, es_surf, psurf=1e-9):
    pc = pf(ec)
    def rhs(r, Y):
        p, m, y = Y
        e_ = einv(max(p, psurf)); e = e_ * CONV; P = p * CONV
        f = 1 - 2 * m / r; cs2 = dpde(e_)
        dp = -(e + P) * (m + 4 * np.pi * r ** 3 * P) / (r * r * f)
        Q = 4 * np.pi / f * (5 * e + 9 * P + (e + P) / cs2) - 6 / (r * r * f) - (2 * (m + 4 * np.pi * r ** 3 * P) / (r * r * f)) ** 2
        F = (1 - 4 * np.pi * r * r * (e - P)) / f
        return [dp / CONV, 4 * np.pi * r * r * e, -(y * y + y * F + r * r * Q) / r]
    ev = lambda r, Y: Y[0] - psurf; ev.terminal = True
    r0 = 1e-4
    s = solve_ivp(rhs, [r0, 60], [pc, 4 / 3 * np.pi * r0 ** 3 * ec * CONV, 2.0], events=ev, rtol=1e-9, atol=1e-14, max_step=0.01)
    R_ = s.t[-1]; M = s.y[1][-1]; y = s.y[2][-1]
    y -= 4 * np.pi * R_ ** 3 * es_surf * CONV / M
    Cc = M / R_; k2 = k2_of(Cc, y)
    return M / MSUN_KM, R_, k2, 2 / 3 * k2 / Cc ** 5

def bag_tidal(B, av, n=45):
    pf, einv = bag_eos(B, av); h = 1e-3
    dpde = lambda e: (pf(e + h) - pf(max(e - h, 4 * B))) / (e + h - max(e - h, 4 * B))
    rows = np.array([tidal_star(pf, einv, dpde, ec, 4 * B) for ec in np.linspace(4 * B * 1.05, 2500, n)])
    i = rows[:, 0].argmax(); r = rows[:i + 1]
    return float(rows[i, 0]), float(np.interp(STDM, r[:, 0], r[:, 3])), float(np.interp(STDM, r[:, 0], r[:, 2]))

# --- minimal xlsx reader (no pandas in the isolated interpreter) ---
def read_xlsx(path, sheet):
    ns = {"m": "http://schemas.openxmlformats.org/spreadsheetml/2006/main"}
    with zipfile.ZipFile(path) as z:
        ss = []
        if "xl/sharedStrings.xml" in z.namelist():
            for si in ET.fromstring(z.read("xl/sharedStrings.xml")).findall("m:si", ns):
                ss.append("".join(t.text or "" for t in si.iter("{%s}t" % ns["m"])))
        wb = ET.fromstring(z.read("xl/workbook.xml"))
        rels = {r.get("Id"): r.get("Target") for r in ET.fromstring(z.read("xl/_rels/workbook.xml.rels"))}
        rid = None
        for s in wb.find("m:sheets", ns):
            if s.get("name") == sheet: rid = s.get("{http://schemas.openxmlformats.org/officeDocument/2006/relationships}id")
        target = rels[rid].lstrip("/"); target = target if target.startswith("xl/") else "xl/" + target
        rows = []
        for row in ET.fromstring(z.read(target)).find("m:sheetData", ns):
            vals = {}
            for c in row:
                col = re.match(r"[A-Z]+", c.get("r")).group(0); t = c.get("t"); v = c.find("m:v", ns)
                if t == "inlineStr": val = "".join(x.text or "" for x in c.iter("{%s}t" % ns["m"]))
                elif v is None: val = None
                elif t == "s": val = ss[int(v.text)]
                elif t in ("str", "b", "e"): val = v.text
                else: val = float(v.text)
                vals[col] = val
            rows.append(vals)
    head = rows[0]; cols = sorted(head, key=lambda k: (len(k), k))
    return [{head[c]: r.get(c) for c in cols} for r in rows[1:]]

XLSX = BUILD / "downloads" / "TRPG_Every_Particle_to_Every_Particle.xlsx"
# ---- rev3: running the published scripts as published (python -I, separate working directory) and the no-knobs neutron-star port
import subprocess, shutil
from scipy.interpolate import interp1d
def arg_after(flag): return sys.argv[sys.argv.index(flag) + 1] if flag in sys.argv else None
OUTDIR = pathlib.Path(arg_after("--outputs")) if arg_after("--outputs") else None
PUBLISHED = ("proton_sets_B", "neutron_stars_no_knobs", "quark_link_pattern", "area_law_and_charges", "expansion_and_entanglement",
             "master_action_M1_containment_potential", "neutron_star_tidal_fix", "build_antiparticle_and_decays")
def clean_paths(s):
    """no local file-system paths on the public page: keep only the last component of any Windows path"""
    return re.sub(r"[A-Za-z]:[\\/][^\s\"']+", lambda m: "..." + "\\" + m.group(0).replace("/", "\\").split("\\")[-1], s)
def published(prefix):
    f_, rc_, o_, e_ = published_raw(prefix); return f_, rc_, clean_paths(o_), clean_paths(e_)
def published_raw(prefix):
    f = sorted((BUILD / "downloads" / "code").glob(prefix + "*.py"))[0]
    if OUTDIR is not None and (OUTDIR / (f.stem + ".txt")).exists():
        rc = None
        for ln in (OUTDIR / "_done.txt").read_text(encoding="utf-8").splitlines():
            if ln.startswith(f.stem + " rc="): rc = int(ln.split("rc=")[1].split()[0])
        return f.name, rc, (OUTDIR / (f.stem + ".txt")).read_text(encoding="utf-8", errors="replace"), (OUTDIR / (f.stem + ".err")).read_text(encoding="utf-8", errors="replace")
    run = HERE / V("script_rundir"); run.mkdir(exist_ok=True)
    shutil.copyfile(BUILD / "downloads" / "TRPG_Every_Particle_to_Every_Particle.xlsx", run / "TRPG_Every_Particle_to_Every_Particle.xlsx")
    pr = subprocess.run([sys.executable, "-I", str(f)], cwd=str(run), capture_output=True, text=True, timeout=7200)
    return f.name, pr.returncode, pr.stdout, pr.stderr
PUB = {p: published(p) for p in PUBLISHED}
for _p in PUB.values(): emit(_p[2]); emit(_p[3][-1500:])
def nums(line): return [float(x) for x in re.findall(r"-?\d+\.\d+|-?\d+", line)]

_nkf = sorted((BUILD / "downloads" / "code").glob("neutron_stars_no_knobs*.py"))[0]
_nkspec = importlib.util.spec_from_file_location("nk_published", _nkf); NKMOD = importlib.util.module_from_spec(_nkspec); _nkspec.loader.exec_module(NKMOD)
def nk_sequence(B, X, gap):
    """Matt's own sequence() from the published neutron_stars_no_knobs script, called directly (its main block is not run)."""
    M_, L_, R_, a4_, es_ = NKMOD.sequence(B, X, gap)
    return float(M_), float(L_), float(R_)


# ================================================================================================ SECTION A: home-page table
SEC_A = "The home-page results table and the neutrino prediction"
SECTIONS["A"] = SEC_A
bend = bending()
check("A1", "A", "Bending of starlight by the Sun", ("results.html", "Results page"),
      ["δ = ∫ ∇⊥n dz,  with n = 1 + 2GM/(rc²)"],
      f"M = {N(MSUN)} kg, R = {N(RSUN)} m, G = {N(G)}, c = {N(C)} → n − 1 = {N(2*G*MSUN/(RSUN*C*C))} at the Sun's edge",
      [f"Deflection = {N(bend)} arcsec (numerical integral along the ray)", rec_line("A1 bending", bend, CL("A1"))],
      f"{N(V('bending_measured'))} arcsec — {SRC('bending_measured')}", f"within {N(V('tol_percent_1')*100,2)}% of the measured value",
      j_rel(bend, V("bending_measured"), V("tol_percent_1")))
age, OLs, Hf = age_gyr()
check("A2", "A", "Age of the universe", (R + ANCH["M11"], "record M11"),
      ["t = ∫ da / (a H(a)),  H = H₀ √(Ω_m a⁻³ + Ω_r a⁻⁴ + Ω_Λ)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, Ω_m = {N(V('Omega_m'))}, Ω_r = {N(V('Omega_r'))}, Ω_Λ = 1 − Ω_m − Ω_r = {N(OLs,5)}",
      [f"Age = {N(age)} Gyr", rec_line("A2 age", age, CL("A2"))],
      f"{N(V('age_measured'))} ± {N(V('age_sigma'))} Gyr — {SRC('age_measured')}", "within 1 standard deviation of the measurement",
      j_sigma(age, V("age_measured"), V("age_sigma")))
bB = V("bag_B_gemini")
mm, r14, _ = bag_sequence(bB, 0.0)
# ---- rev3 A3: neutron stars with nothing adjusted (record Q4, home page)
Bp = V("bag_B_proton"); gaps = V("pairing_gap_range"); Xc = V("no_knobs_scale")
nk = {(Bv, gp): nk_sequence(Bv, Xc, gp) for Bv in Bp for gp in gaps}
nkM = [v[0] for v in nk.values()]; nkL = [v[1] for v in nk.values()]; nkR = [v[2] for v in nk.values()]
nk_sub = (f"B = {', '.join(N(b) for b in Bp)} MeV/fm\u00b3 (set by the proton, stage B9); pairing gap \u0394 = {N(gaps[0])} and {N(gaps[1])} MeV; scale factor X = {N(Xc)}; "
          f"strange quark mass {N(V('no_knobs_ms'))} MeV; \u03b1_s from {N(V('no_knobs_alpha_s_MZ'))} at M_Z by two-loop running")
nk_eq = ["p(\u03bc) = (3/4\u03c0\u00b2)[a\u2084\u03bc\u2074 \u2212 (m_s\u00b2 \u2212 4\u0394\u00b2)\u03bc\u00b2] \u2212 B,   a\u2084 = 1 \u2212 2\u03b1_s/\u03c0;   then the TOV equations"]
mt, ms_ = V("ns_mass_test"), V("ns_mass_sigma")
check("A3a", "A", "Heaviest neutron star (B from the proton, quark pairing)", (R + ANCH["Q4"], "record Q4"), nk_eq, nk_sub,
      [f"Maximum mass from {N(min(nkM),3)} to {N(max(nkM),3)} solar masses"] + [f"B = {N(k_[0])}, Δ = {N(k_[1])} MeV: M_max = {N(v_[0],3)}, Λ = {N(v_[1],3)}, R = {N(v_[2],3)} km" for k_, v_ in nk.items()] + [rec_line("A3a low", min(nkM), CL("Q4mlo"), 3), rec_line("A3a high", max(nkM), CL("Q4mhi"), 3)],
      f"at least {N(mt)} solar masses (the measured star is {N(V('ns_mass_measured'))} \u00b1 {N(ms_)}) \u2014 {SRC('ns_mass_test')}", "every case in the proton band at or above the test mass",
      "PASS" if min(nkM) >= mt else ("PARTIAL" if min(nkM) >= mt - ms_ else "FAIL"))
rlo, rhi = V("ns_radius_lo"), V("ns_radius_hi")
inside = [rlo <= r_ <= rhi for r_ in nkR]
check("A3b", "A", "Radius of a standard-mass neutron star (same cases)", (R + ANCH["Q4"], "record Q4"), ["same star sequences; R at M = 1.4 solar masses"], nk_sub,
      [f"R({N(STDM)}) from {N(min(nkR),3)} to {N(max(nkR),3)} km", rec_line("A3b low", min(nkR), CL("Q4Rlo"), 3), rec_line("A3b high", max(nkR), CL("Q4Rhi"), 3)],
      f"{N(rlo)} to {N(rhi)} km \u2014 {SRC('ns_radius_lo')}", "every case inside the measured range",
      "PASS" if all(inside) else ("PARTIAL" if any(inside) else "FAIL"),
      "Part of the predicted range sits below the lower edge of the quoted measurement." if not all(inside) else "")
check("A3c", "A", "Tidal deformability \u039b of a standard-mass neutron star (same cases)", (R + ANCH["Q4"], "record Q4"), ["\u039b = (2/3) k\u2082 / C\u2075, from the same sequences (log-interpolated in \u039b)"], nk_sub,
      [f"\u039b({N(STDM)}) from {N(min(nkL),3)} to {N(max(nkL),3)}", rec_line("A3c low", min(nkL), CL("Q4Llo"), 3), rec_line("A3c high", max(nkL), CL("Q4Lhi"), 3)],
      f"\u039b at or below {N(V('lambda_limit'))} \u2014 {SRC('lambda_limit')}", "every case at or below the limit", "PASS" if max(nkL) <= V("lambda_limit") else "FAIL")
rows_nk = [nums(ln) for ln in PUB["neutron_stars_no_knobs"][2].splitlines() if re.match(r"^\s+\d+\.\d\s+\d\.\d\s+\d\.\d{3}\s", ln)]
strong = [r_[3] for r_ in rows_nk[:9]] if len(rows_nk) >= 9 else []
check("A3d", "A", "The strong force alone, with no quark pairing", (R + ANCH["Q4"], "record Q4"), ["same equation of state with \u0394 = 0"],
      "from the published script's own printed table (nine cases: three values of B, three scale factors)",
      ([f"Maximum mass from {N(min(strong),3)} to {N(max(strong),3)} solar masses", rec_line("A3d low", min(strong), CL("Q4strong_lo"), 3), rec_line("A3d high", max(strong), CL("Q4strong_hi"), 3)] if strong else ["The published script's table could not be read."]),
      f"at least {N(mt)} solar masses \u2014 {SRC('ns_mass_test')}", "heaviest star at or above the test mass", ("PASS" if strong and min(strong) >= mt else "FAIL"),
      "Fails the two-solar-mass test, as the record says: this is why the pairing gap is needed.")
check("A3e", "A", "The predicted quark pairing gap", (R + ANCH["Q4"], "record Q4"), ["\u0394 such that both neutron-star tests pass"],
      f"gap range {N(gaps[0])} to {N(gaps[1])} MeV; theory estimates {N(V('pairing_theory_range')[0])} to {N(V('pairing_theory_range')[1])} MeV",
      ["No measurement of the pairing gap exists."], "not measured", "a measured pairing gap", "NOT TESTABLE YET",
      "A prediction that can fail: it has not been measured, as the record says.")
# leptons
Qlep = Qk(me, mmu, mtau)
a_, b_ = math.sqrt(me), math.sqrt(mmu); s_ = a_ + b_
A_, B_, C_ = 1 / 3, -4 / 3 * s_, (a_ * a_ + b_ * b_) - 2 / 3 * s_ * s_
x_ = (-B_ + math.sqrt(B_ * B_ - 4 * A_ * C_)) / (2 * A_); tau_pred = x_ * x_
Mt = math.sqrt(mtau) / (1 + math.sqrt(2) * math.cos(D29))
mu_pred = (Mt * (1 + math.sqrt(2) * math.cos(D29 + 4 * math.pi / 3))) ** 2
e_pred = (Mt * (1 + math.sqrt(2) * math.cos(D29 + 2 * math.pi / 3))) ** 2
tl = V("tol_lepton_mass")
check("A4a", "A", "Tau mass from the electron and muon (Koide)", (R + ANCH["K7"], "record K7"),
      ["(m_e + m_μ + m_τ) / (√m_e + √m_μ + √m_τ)² = 2/3, solved for m_τ"],
      f"m_e = {N(me,8)} MeV, m_μ = {N(mmu,10)} MeV; Q on the measured trio = {N(Qlep,7)}",
      [f"Predicted m_τ = {N(tau_pred,6)} MeV", rec_line("A4a tau", tau_pred, CL("A4a"), 6)],
      f"{N(mtau,6)} MeV — {SRC('m_tau')}", f"within {N(tl*100,2)}%", j_rel(tau_pred, mtau, tl))
check("A4b", "A", "Electron and muon masses from the tau alone", (R + ANCH["K7"], "record K7"),
      ["√m_k = M(1 + √2 cos(δ + 2πk/3)),  δ = 2/9,  M fixed by m_τ"],
      f"δ = {fr(V('brannen_delta_num'), V('brannen_delta_den'))} = {N(D29,6)}; M from m_τ = {N(mtau,6)} MeV",
      [f"Predicted m_e = {N(e_pred,6)} MeV ({N((e_pred/me-1)*100,2)}% off); predicted m_μ = {N(mu_pred,7)} MeV ({N((mu_pred/mmu-1)*100,2)}% off)",
       "Three masses from one input scale (plus the fixed twist)."],
      f"m_e = {N(me,8)} MeV, m_μ = {N(mmu,10)} MeV — {SRC('m_e')}", f"each within {N(tl*100,2)}%",
      "PASS" if j_rel(e_pred, me, tl) == "PASS" and j_rel(mu_pred, mmu, tl) == "PASS" else ("PARTIAL" if max(abs(e_pred/me-1), abs(mu_pred/mmu-1)) <= PF*tl else "FAIL"))
cab = math.sqrt(V("q_d") / V("q_s"))
check("A5", "A", "Cabibbo angle from quark masses (label: already known, Gatto, Sartori and Tonin 1968)", (R + ANCH["K4"], "record K4"),
      ["θ_C = e^S = √(m_d / m_s),  S = ½ ln(m_d / m_s)"],
      f"m_d = {N(V('q_d'))} MeV, m_s = {N(V('q_s'))} MeV ({SRC('q_d')})",
      [f"√(m_d/m_s) = {N(cab,4)}", rec_line("A5 cabibbo", cab, CL("A5"))],
      f"{N(V('ckm_Vus'))} — {SRC('ckm_Vus')}", f"within {N(V('tol_percent_1')*100,2)}% (the quark-mass uncertainty)", j_rel(cab, V("ckm_Vus"), V("tol_percent_1")),
      "The relation is Gatto, Sartori and Tonin (1968); it is known physics, not new here.")
# tau lifetime
GF, hb = V("GF_MeV2"), V("hbar_MeVs")
tau_mu_c = hb / (GF ** 2 * mmu ** 5 / (192 * math.pi ** 3)); S_tm = 0.5 * math.log(mtau / mmu)
tau_tau_c = tau_mu_c * V("BR_tau_e") * math.exp(-10 * S_tm) * 1e15
check("A6", "A", "Tau lifetime from the muon (label: already known, lepton universality)", (R + ANCH["K10"], "record K10"),
      ["τ_μ = 192π³ℏ / (G_F² m_μ⁵);   τ_τ = τ_μ · BR · e^(−10 S),  S = ½ ln(m_τ/m_μ)"],
      f"G_F = {N(GF)} MeV⁻², m_μ = {N(mmu,7)} MeV, m_τ = {N(mtau,6)} MeV, BR(τ→eνν) = {N(V('BR_tau_e'))} → S = {N(S_tm,5)}",
      [f"Muon lifetime from Fermi theory = {N(tau_mu_c*1,5)} s (measured {N(V('tau_mu_s'),7)} s: {N((tau_mu_c/V('tau_mu_s')-1)*100,2)}% off)",
       f"Tau lifetime = {N(tau_tau_c,4)} fs", rec_line("A6 tau lifetime", tau_tau_c, CL("A6"))],
      f"{N(V('tau_tau_fs'))} fs — {SRC('tau_tau_fs')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(tau_tau_c, V("tau_tau_fs"), V("tol_percent_1")),
      "e^(−10 S) is the same as (m_μ/m_τ)⁵, the Fermi-theory mass scaling, so the exponent is not a fitted number.")
# anomaly sums (used by A7, B5, M16)
Yv = {k: fr_ for k, fr_ in V("anomaly_Y").items()}
from fractions import Fraction as Fr
Y = {k: Fr(a, b) for k, (a, b) in Yv.items()}
cond = {"SU(3)²U(1)": 2 * Y["YQ"] - Y["Yu"] - Y["Yd"], "SU(2)²U(1)": 3 * Y["YQ"] + Y["YL"],
        "gravity²U(1)": 6 * Y["YQ"] - 3 * Y["Yu"] - 3 * Y["Yd"] + 2 * Y["YL"] - Y["Ye"],
        "U(1)³": 6 * Y["YQ"] ** 3 - 3 * Y["Yu"] ** 3 - 3 * Y["Yd"] ** 3 + 2 * Y["YL"] ** 3 - Y["Ye"] ** 3}
anom_ok = all(v == 0 for v in cond.values())
cols = {n: n * (Fr(2, 3) - Fr(1, 3)) + (0 - 1) for n in V("colours_tested")}
only3 = [n for n, v in cols.items() if v == 0]
anom_lines = [f"{k} = {N(float(v),3)}" for k, v in cond.items()]
col_line = "Sum of electric charges in one generation for N colours: " + "; ".join(f"N = {I(n)} → {N(float(v),3)}" for n, v in cols.items())
check("A7", "A", "Quark charges and three colours (anomaly cancellation)", (R + ANCH["M16"], "record M16"),
      ["2Y_Q − Y_u − Y_d = 0;  3Y_Q + Y_L = 0;  6Y_Q − 3Y_u − 3Y_d + 2Y_L − Y_e = 0;  6Y_Q³ − 3Y_u³ − 3Y_d³ + 2Y_L³ − Y_e³ = 0",
       "charge sum per generation = N (Q_u + Q_d) + Q_e + Q_ν = 0"],
      "hypercharges Y_Q = " + fr(*V("anomaly_Y")["YQ"]) + ", Y_u = " + fr(*V("anomaly_Y")["Yu"]) + ", Y_d = " + fr(*V("anomaly_Y")["Yd"]) + ", Y_L = " + fr(*V("anomaly_Y")["YL"]) + ", Y_e = " + fr(*V("anomaly_Y")["Ye"]),
      anom_lines + [col_line, f"Only N = {I(only3[0]) if len(only3)==1 else '?'} gives zero."],
      "all four anomaly sums must be exactly zero (Geng and Marshak 1989; Minahan, Ramond and Warner 1990)", "exact rational arithmetic: zero means zero",
      "PASS" if anom_ok and len(only3) == 1 else "FAIL")
# ringdown
Msrc, zz, aa = V("ringdown_M_source"), V("ringdown_z"), V("ringdown_spin"); k = V("kerr_fit"); Msun_s = V("Msun_s")
Mrem = Msrc * (1 + zz)
Mw = k[0] - k[1] * (1 - aa) ** k[2]; Qf = k[3] + k[4] * (1 - aa) ** (-k[5])
w_ = Mw / (Mrem * Msun_s); f_ = w_ / (2 * math.pi); tau_ = 2 * Qf / w_ * 1e3
dt_echo = 4 * Mrem * Msun_s * (1 + 1 / math.sqrt(1 - aa * aa)) * math.log(Mrem * V("Msun_m") / V("planck_length_m"))
sf, st = V("ringdown_f_sigma"), V("ringdown_tau_sigma")
check("A8a", "A", "Black-hole ringdown tone and fade (GW150914)", (R + ANCH["K9"], "record K9"),
      ["Mω = a₀ − a₁ (1−χ)^a₂;  Q = b₀ + b₁ (1−χ)^(−b₂);  f = ω/2π;  τ = 2Q/ω"],
      f"remnant mass = {N(Msrc)} × (1 + z) = {N(Mrem,4)} solar masses (z = {N(zz)}), spin χ = {N(aa)}; fit constants {', '.join(N(v,5) for v in k)} (Berti, Cardoso and Will 2006)",
      [f"f = {N(f_,3)} Hz, τ = {N(tau_,3)} ms", rec_line("A8 frequency", f_, CL("A8f"), 3)],
      f"f = {N(V('ringdown_f_measured'))} ± {N(sf)} Hz, τ = {N(V('ringdown_tau_measured'))} ± {N(st)} ms — {SRC('ringdown_f_measured')}",
      "each within 1 standard deviation",
      "PASS" if j_sigma(f_, V("ringdown_f_measured"), sf) == "PASS" and j_sigma(tau_, V("ringdown_tau_measured"), st) == "PASS" else "FAIL")
check("A8b", "A", "No echoes after the ringdown", (R + ANCH["K9"], "record K9"),
      ["Δt = 4M(1 + 1/√(1−χ²)) ln(M/ℓ_P)   (delay if a reflecting surface sat one Planck length outside the horizon)"],
      f"M = {N(Mrem,4)} solar masses, χ = {N(aa)}, ℓ_P = {N(V('planck_length_m'))} m",
      [f"Echo delay that would appear = {N(dt_echo,3)} s. TRPG predicts there is no such surface, so no echoes."],
      "searches by LIGO/Virgo and by Westerweck and others (2018) found no significant echoes (this run does not repeat that search; the result is as the record cites it)",
      "a negative prediction is consistent with a null result", "PASS",
      "Consistent with the null searches. This is a prediction of absence, so it can still fail if echoes are found later.")
# neutrinos
dm21, dm31, dm32 = V("dm21"), V("dm31"), V("dm32_inverted")
def Qs(m, s): rt = np.array(s) * np.sqrt(m); return m.sum() / rt.sum() ** 2
def roots_of(f):
    grid = np.logspace(-6, 0, 4000); vals = [f(x) for x in grid]
    return [brentq(f, grid[i], grid[i + 1]) for i in range(len(grid) - 1) if np.sign(vals[i]) != np.sign(vals[i + 1])]
nrm = []
for sg in [(1, 1, 1), (-1, 1, 1)]:
    f = lambda m1, sg=sg: Qs(np.array([m1, math.sqrt(m1 ** 2 + dm21), math.sqrt(m1 ** 2 + dm31)]), sg) - Q23
    for m1 in roots_of(f): nrm.append((sg, m1, math.sqrt(m1 ** 2 + dm21), math.sqrt(m1 ** 2 + dm31)))
inv = []
for sg in [(1, 1, 1), (1, 1, -1), (1, -1, 1), (-1, 1, 1)]:
    def f(m3, sg=sg):
        m2 = math.sqrt(m3 ** 2 - dm32); m1 = math.sqrt(m2 ** 2 - dm21); return Qs(np.array([m1, m2, m3]), sg) - Q23
    for m3 in roots_of(f):
        m2 = math.sqrt(m3 ** 2 - dm32); inv.append((sg, math.sqrt(m2 ** 2 - dm21), m2, m3))
pos = [r for r in nrm if r[0] == (-1, 1, 1)]
pred = pos[0] if pos else None
sum_pred = sum(pred[1:]) if pred else float("nan")
lim = V("nu_sum_limit")
nrm_lines = [f"Normal ordering, square-root signs ({', '.join(('−' if s<0 else '+') for s in sg)}): m₁, m₂, m₃ = {N(a*1e3,3)}, {N(b*1e3,3)}, {N(c*1e3,3)} meV, sum = {N((a+b+c),4)} eV" for sg, a, b, c in nrm]
inv_lines = [f"Inverted ordering, signs ({', '.join(('−' if s<0 else '+') for s in sg)}): sum = {N(a+b+c,3)} eV" for sg, a, b, c in inv]
check("A9a", "A", "Neutrino masses from the lepton rule (normal ordering)", (R + ANCH["M2"], "record M2"),
      ["Q = (m₁ + m₂ + m₃)/(±√m₁ ± √m₂ ± √m₃)² = 2/3,  m₂² = m₁² + Δm²₂₁,  m₃² = m₁² + Δm²₃₁"],
      f"Δm²₂₁ = {N(dm21)} eV², Δm²₃₁ = {N(dm31)} eV² ({SRC('dm21')})",
      (nrm_lines + [rec_line("A9 neutrino sum", sum_pred, CL("A9")), rec_line("A9 lightest mass (meV)", pred[1] * 1e3, V("nu_light_stated_meV"), 3)]) if pred else nrm_lines,
      f"not yet measured; the record's stated prediction is {', '.join(N(v) for v in V('nu_masses_stated'))} meV", "the solution with one negative square root must reproduce the stated masses within " + N(V('tol_percent_5')*100,2) + "%",
      ("PASS" if pred and all(abs(p*1e3/s-1) <= V("tol_percent_5") for p, s in zip(pred[1:], V("nu_masses_stated"))) else "FAIL"),
      "This stage checks the arithmetic of the prediction, not nature. Nature's answer is in the next two stages.")
check("A9b", "A", "Does the predicted neutrino mass sum respect the cosmology limit?", (R + ANCH["M2"], "record M2"),
      ["Σm_ν (normal-ordering solution) < Σm_ν limit;   Σm_ν (inverted-ordering solutions) < limit?"],
      f"limit = {N(lim)} eV ({SRC('nu_sum_limit')}); predicted sum = {N(sum_pred,4)} eV",
      inv_lines if inv else ["no inverted-ordering solution found"],
      f"{N(lim)} eV upper limit — {SRC('nu_sum_limit')}", "predicted sum below the limit, and every inverted solution above it (ruling inverted ordering out)",
      "PASS" if sum_pred < lim and inv and all(sum(r[1:]) > lim for r in inv) else "FAIL")
num_min = math.sqrt(dm21) + math.sqrt(dm31)
check("A9d", "A", "How close the predicted sum is to the smallest sum normal ordering allows", (R + ANCH["M2"], "record M2"),
      ["Σm_min = √Δm²₂₁ + √Δm²₃₁   (lightest neutrino massless)"], f"Δm²₂₁ = {N(dm21)} eV², Δm²₃₁ = {N(dm31)} eV²",
      [f"Smallest allowed sum = {N(num_min,4)} eV; predicted sum = {N(sum_pred,4)} eV; limit = {N(lim)} eV", rec_line("A9 minimum sum", num_min, CL("numin"), 4)],
      "the record: the sum is only just above the smallest allowed and just under the limit", "predicted sum between the smallest allowed sum and the limit", j_range(sum_pred, num_min, lim),
      "The distinctive TRPG number is the lightest mass, not the sum.")
check("A9c", "A", "Neutrino ordering: normal, as predicted", (R + ANCH["M2"], "record M2"),
      ["prediction: normal ordering (m₁ < m₂ < m₃); inverted ordering would falsify it"],
      "no data yet", ["Awaiting the JUNO reactor experiment (running since 2025) and the sky surveys."],
      "not yet measured", "the JUNO measurement of the ordering", "NOT TESTABLE YET",
      "This is the prediction that could prove TRPG wrong. If JUNO finds inverted ordering, or a sum clearly away from the predicted one, this piece of TRPG fails.")

# ================================================================================================ SECTION B: master equation terms
SECTIONS["B"] = "The master equation, term by term"
# area law on the plain lattice
def lattice(Nn, m, crossing=False, amp=0.5):
    idx = lambda i, j: i * Nn + j; K = np.zeros((Nn * Nn, Nn * Nn))
    for i in range(Nn):
        for j in range(Nn):
            a = idx(i, j)
            for di, dj in ((1, 0), (0, 1)):
                ii, jj = i + di, j + dj
                if ii < Nn and jj < Nn:
                    xm, ym = (i + ii) / 2, (j + jj) / 2
                    kk = 1 + amp * math.cos(2 * math.pi * xm / 6) * math.cos(2 * math.pi * ym / 6) if crossing else 1.0
                    b = idx(ii, jj); K[a, b] = K[b, a] = -kk; K[a, a] += kk; K[b, b] += kk
            K[a, a] += m * m
    w, U = np.linalg.eigh(K); X = 0.5 * (U / np.sqrt(w)) @ U.T; P = 0.5 * (U * np.sqrt(w)) @ U.T
    return idx, X, P
def lattice_area(Nn, m):
    idx = lambda i, j: i * Nn + j; K = np.zeros((Nn * Nn, Nn * Nn))
    for i in range(Nn):
        for j in range(Nn):
            a = idx(i, j); K[a, a] = 4 + m * m
            for di, dj in ((1, 0), (0, 1)):
                ii, jj = i + di, j + dj
                if ii < Nn and jj < Nn: b = idx(ii, jj); K[a, b] = K[b, a] = -1
    w, U = np.linalg.eigh(K); X = 0.5 * (U / np.sqrt(w)) @ U.T; P = 0.5 * (U * np.sqrt(w)) @ U.T
    return idx, X, P
def ent(Nn, m, crossing, L, amp=0.5, area_script=False):
    idx, X, P = lattice_area(Nn, m) if area_script else lattice(Nn, m, crossing, amp); out = []
    for l in L:
        o = (Nn - l) // 2; A = [idx(i, j) for i in range(o, o + l) for j in range(o, o + l)]
        nu = np.sqrt(np.clip(np.linalg.eigvals(X[np.ix_(A, A)] @ P[np.ix_(A, A)]).real, 0.25 + 1e-12, None))
        out.append(sum((n + .5) * math.log(n + .5) - (n - .5) * math.log(n - .5) for n in nu if n > 0.5 + 1e-9))
    return np.array(out)
Ng, mg = int(V("grid_N")), V("grid_mass"); Ls = np.array(V("area_blocks"))
S_area = ent(Ng, mg, False, Ls, area_script=True)
cb = np.polyfit(4 * Ls, S_area, 1); cv = np.polyfit(Ls * Ls, S_area, 1)
rb = float(np.sqrt(np.mean((np.polyval(cb, 4 * Ls) - S_area) ** 2))); rv = float(np.sqrt(np.mean((np.polyval(cv, Ls * Ls) - S_area) ** 2)))
Gcoef = 1 / (4 * cb[0])
check("B1a", "B", "Gravity term: the area law holds in the wave medium", (R + ANCH["M16"], "record M16"),
      ["S(region) = c × (boundary length)   rather than   c × (area inside)",
       "G = a² / (4c)   (Bekenstein–Jacobson relation)"],
      f"{I(Ng)}×{I(Ng)} lattice of coupled oscillators, mass {N(mg)}, square blocks of side {', '.join(I(l) for l in Ls)}",
      [f"Entanglement per unit boundary: {', '.join(N(s/(4*l),3) for s, l in zip(S_area, Ls))}",
       f"Boundary-fit error = {N(rb,3)}; volume-fit error = {N(rv,3)}; ratio = {N(rv/rb,3)}", f"c = {N(cb[0],3)}, so G = {N(Gcoef,3)} a²",
       rec_line("B1 G coefficient", Gcoef, CL("B1"), 3)],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993): the area law is the known result for this kind of system",
      f"the boundary fit must beat the volume fit by at least a factor {N(V('area_ratio_min'))}", "PASS" if rv / rb >= V("area_ratio_min") else "FAIL")
check("B1b", "B", "Gravity term: Newton's G in real units", (R + ANCH["P1"], "record P1"),
      ["G = a² × (the coefficient above), with a = the medium's smallest spacing"],
      "a is not fixed by any measurement; the record takes it to be the Planck length",
      ["Not computed: the coefficient is in units of a², and a is not derived."], "G is measured, but TRPG does not yet predict a", "a derived value of a", "NOT TESTABLE YET")
Hs = V("H0") * 1e3 / V("Mpc_m"); rho_c = 3 * Hs ** 2 / (8 * math.pi * G)
uL = V("Omega_L_stated") * rho_c * V("c_exact") ** 2; uL_MeVfm3 = uL / V("MeV_J") * 1e-45
EL = (uL_MeVfm3 * V("hbarc_MeVfm") ** 3) ** 0.25 * 1e9
check("B2a", "B", "Cosmic baseline resistance −ρ_Λ", (R + ANCH["M8"], "record M8"),
      ["ρ_c = 3H₀²/(8πG);   u_Λ = Ω_Λ ρ_c c²;   E = (u_Λ (ℏc)³)^(1/4)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, Ω_Λ = {N(V('Omega_L_stated'))} ({SRC('H0')})",
      [f"u_Λ = {N(uL_MeVfm3,3)} MeV/fm³; energy scale = {N(EL,3)} meV", rec_line("B2 density", uL_MeVfm3, V("de_density_stated"), 3), rec_line("B2 scale", EL, V("de_scale_meV_stated"), 3)],
      "taken from the measured dark-energy density: an input, not a TRPG prediction", "not applicable", "INPUT",
      "TRPG does not yet derive this value.")
nu_lo, nu_hi = min(V("nu_masses_stated")), max(V("nu_masses_stated"))
check("B2b", "B", "Dark-energy scale versus the neutrino masses", (R + ANCH["M8"], "record M8"),
      ["m_ν,min ≤ E_Λ ≤ m_ν,max"],
      f"E_Λ = {N(EL,3)} meV; predicted neutrino masses {N(nu_lo)} to {N(nu_hi)} meV",
      [f"E_Λ lies {'inside' if nu_lo <= EL <= nu_hi else 'outside'} the predicted range."], "known puzzle: Fardon, Nelson and Weiner (2004)", "inside the range",
      j_range(EL, nu_lo, nu_hi), "A coincidence of scale, not a derivation.")
hbc = V("hbarc_MeVfm"); B_list = V("bag_B_proton")
gl_list = [V("glueball_lo"), V("glueball_hi")]
rows_lv = []
for Bv in B_list:
    for mgl in gl_list:
        Bn = Bv * hbc ** 3; lam = mgl ** 4 / (16 * Bn); v_ = math.sqrt(mgl ** 2 / (2 * lam)); rows_lv.append((Bv, mgl, lam, v_, hbc / mgl))
lam_lo = min(r[2] for r in rows_lv); lam_hi = max(r[2] for r in rows_lv); v_lo = min(r[3] for r in rows_lv); v_hi = max(r[3] for r in rows_lv)
wall = rows_lv[0][4]; core_hi = V("core_radius_hi")
cen_rows = [r for r in rows_lv if r[0] == B_list[1]]; lam_c_lo = min(r[2] for r in cen_rows); lam_c_hi = max(r[2] for r in cen_rows); v_c = [r[3] for r in cen_rows]
check("B3a", "B", "Containment field P: stiffness and depth", (R + ANCH["M1"], "record M1"),
      ["V(P) = (λ/4) v⁴ (P² − 1)²;   depth V(0) − V(1) = B;   glueball mass m² = 2λv²;   so λ = m⁴/(16B),  v = m/√(2λ)"],
      "B = " + ", ".join(N(b) for b in B_list) + " MeV/fm³ (set by the proton, stage B9); glueball = " + ", ".join(N(m) for m in gl_list) + " MeV (Morningstar and Peardon 1999; Chen and others 2006)",
      [f"λ from {N(lam_lo,3)} to {N(lam_hi,3)}; v from {N(v_lo,3)} to {N(v_hi,3)} MeV", f"Wall thickness ≈ ℏc/m = {N(wall,3)} fm",
       f"At the central B = {N(B_list[1])}: λ from {N(lam_c_lo,3)} to {N(lam_c_hi,3)}; v from {N(min(v_c),3)} to {N(max(v_c),3)} MeV",
       rec_line("B3 lambda_low", lam_c_lo, CL("B3lo"), 3), rec_line("B3 lambda_high", lam_c_hi, CL("B3hi"), 3), rec_line("B3 v", sum(v_c) / len(v_c), CL("B3v"), 3)],
      "no measurement of λ or v exists; this stage checks that the wall is thin enough for the bag-model results to carry over",
      "wall thickness must be well below the proton core radius " + N(core_hi) + " fm", "PASS" if wall < 0.5 * core_hi else "FAIL",
      "The quartic shape (P² − 1)² is the simplest choice, not a derived one (record M1).")
eps = (V("gluon_vev_eps_prefactor_num") / V("gluon_vev_eps_prefactor_den")) * V("svz_G2") * 1e12 / hbc ** 3
ratio_svz = eps / B_list[1]
check("B3b", "B", "Containment depth versus the standard QCD vacuum energy", (R + ANCH["M1"], "record M1"),
      ["ε_SVZ = (9/32) ⟨α_s G²/π⟩;   compare B"],
      f"⟨α_s G²/π⟩ = {N(V('svz_G2'))} GeV⁴ ({SRC('svz_G2')}); B = {N(B_list[1])} MeV/fm³",
      [f"ε_SVZ = {N(eps,3)} MeV/fm³, which is {N(ratio_svz,3)} times B"],
      f"the standard QCD estimate is {N(eps,3)} MeV/fm³", f"B within a factor {N(V('svz_factor'))} of the QCD estimate",
      "PASS" if ratio_svz <= V("svz_factor") else "FAIL", "The record calls this a known, unresolved gap in mainstream physics. It is shown here as a FAIL because it is one.")
rng_ = np.random.default_rng(int(V("random_seed")))
mvals = rng_.uniform(0.1, 100, 6)
scale_f = V("scale_factors"); S_ok = True; cosh_ok = True
for lam_ in scale_f:
    for i_ in range(5):
        a1, a2 = mvals[i_], mvals[i_ + 1]
        S1 = 0.5 * math.log(a1 / a2); S2 = 0.5 * math.log(lam_ * a1 / (lam_ * a2))
        if abs(S1 - S2) > 1e-12: S_ok = False
        if abs((a1 + a2) / math.sqrt(a1 * a2) - 2 * math.cosh(S1)) > 1e-9: cosh_ok = False
check("B4", "B", "Scale layer σ: the distance S is the only one that gives no particle a privileged unit", (R + ANCH["M3"], "record M3"),
      ["S_ij = ½ ln(m_i/m_j);  S(λm_i, λm_j) = S(m_i, m_j);  M_ij = (m_i + m_j)/√(m_i m_j) = 2 cosh S_ij"],
      f"random positive masses; rescaling factors λ = {', '.join(N(v) for v in scale_f)}",
      [f"S unchanged by rescaling in every trial: {'yes' if S_ok else 'no'}", f"M_ij = 2 cosh S_ij in every trial: {'yes' if cosh_ok else 'no'}"],
      "exact identities", f"agreement to {N(V('tol_identity'))} relative", "PASS" if S_ok and cosh_ok else "FAIL",
      "This checks the algebra. That S is forced by the symmetry is the record's argument (M3).")
check("B5", "B", "Forces: charges and colours", (R + ANCH["M16"], "record M16"),
      ["−¼ F² for U(1) × SU(2) × SU(3);  charges forced by anomaly cancellation"],
      "see stage A7 for the sums", [f"Anomaly sums all zero: {'yes' if anom_ok else 'no'}.", "The force groups themselves are chosen, not derived: the record says so (T2 is PARTIAL)."],
      "anomaly cancellation (Geng and Marshak 1989; Minahan, Ramond and Warner 1990)", "exact zero", "PARTIAL",
      "Charges and colour count are forced; the groups U(1), SU(2), SU(3) are not.")
check("B6", "B", "Matter term ψ̄ iγDψ", (R + ANCH["P1"], "record P1"),
      ["ψ̄ iγ^μ D_μ ψ"], "Dirac (1928)", ["Taken unchanged from established physics; TRPG adds no number here."], "Dirac (1928)", "not applicable", "INHERITED")
# lepton matrix
r_ = 1 / math.sqrt(2); d_ = D29; M0 = (math.sqrt(me) + math.sqrt(mmu) + math.sqrt(mtau)) / 3
w_c = np.exp(1j * d_); Sm = M0 * np.array([[1, r_ * w_c, r_ * np.conj(w_c)], [r_ * np.conj(w_c), 1, r_ * w_c], [r_ * w_c, r_ * np.conj(w_c), 1]])
ev = np.sort(np.linalg.eigvalsh(Sm)) ** 2
vv = np.sqrt([me, mmu, mtau]); cc = vv.mean() * np.ones(3); pp = vv - cc
errs = [abs(ev[0] / me - 1), abs(ev[1] / mmu - 1), abs(ev[2] / mtau - 1)]
check("B7a", "B", "Mass couplings: the lepton matrix gives the three masses", (R + ANCH["M2"], "record M2"),
      ["S_ℓ = M₀ × circulant(1, e^{iδ}/√2, e^{−iδ}/√2);  masses = (eigenvalues of S_ℓ)²"],
      f"δ = {fr(V('brannen_delta_num'), V('brannen_delta_den'))} = {N(D29,6)}; M₀ = (√m_e + √m_μ + √m_τ)/3 = {N(M0,6)} MeV½",
      [f"Eigenvalues squared: {N(ev[0],6)}, {N(ev[1],7)}, {N(ev[2],6)} MeV"],
      f"{N(me,6)}, {N(mmu,7)}, {N(mtau,6)} MeV — {SRC('m_e')}", f"each within {N(V('tol_percent_1')*100,2)}%",
      "PASS" if all(e <= V("tol_percent_1") for e in errs) else "FAIL",
      "M₀ is computed from the measured masses here, so this stage shows the matrix form is consistent with the three masses, not that it predicts them from nothing.")
check("B7b", "B", "Self part equals relational part (Koide's balance)", (R + ANCH["M2"], "record M2"),
      ["|common part|² = |relational part|² of the square-root-mass vector"],
      f"common part² = {N(cc@cc,6)}, relational part² = {N(pp@pp,6)} MeV",
      [f"Ratio = {N((pp@pp)/(cc@cc),6)}; common part² versus the proton mass: {N((cc@cc/V('m_proton')-1)*100,3)}%"],
      "Koide (1981); Foot (1994)", f"ratio equal to 1 within {N(V('tol_lepton_mass')*100,2)}%", j_rel((pp@pp) / (cc@cc), 1.0, V("tol_lepton_mass")))
# ================================================================================================ SECTION C: chapter K
SECTIONS["C"] = "Chapter K: the fixes, one by one"
gm_mmax, gm_r14, gm_p150 = gemini_sequence()
check("K1a", "C", "Gemini's own neutron-star equation solved properly", (R + ANCH["K1"], "record K1"),
      ["dp/dr = −(ε + p)(m + 4πr³p)/(r(r − 2m)),  dm/dr = 4πr²ε  with Gemini's low-density formula"],
      f"B = {N(V('bag_B_gemini'))} MeV/fm³, low-density term p = {N(V('gemini_p_scale'))} ε (ε/4B)^{N(V('gemini_p_power'))}",
      [f"R({N(STDM)}) = {N(gm_r14,3)} km; maximum mass {N(gm_mmax,3)} solar masses", f"Pressure at ε = {N(V("nuclear_density_probe"))} MeV/fm³ is {N(gm_p150,3)} MeV/fm³ (real nuclear matter is about {N(V('nuclear_pressure_range')[0])} to {N(V('nuclear_pressure_range')[1])})", rec_line("K1a R1.4", gm_r14, CL("K1aR"), 3)],
      f"{N(V('ns_radius_lo'))} to {N(V('ns_radius_hi'))} km — {SRC('ns_radius_lo')}", "inside the measured range", j_range(gm_r14, V("ns_radius_lo"), V("ns_radius_hi")),
      "This is the version that was wrong. It is shown so the failure is on the page.")
Es, Ps = sly_tables()
sly_pf = lambda e: float(np.interp(e, Es, Ps)); sly_inv = lambda p: float(np.interp(p, Ps, Es))
sl = np.array([(ec,) + tov(sly_pf, sly_inv, ec, psurf=1e-12) for ec in np.linspace(V("sly_scan")[0], V("sly_scan")[1], 57)])
i_ = sl[:, 2].argmax(); jj = int(np.argmin(abs(sl[:i_, 2] - STDM)))
sly_m, sly_r = float(sl[i_, 2]), float(sl[jj, 1])
check("K1b", "C", "Solver check on a standard real equation of state (SLy)", (R + ANCH["K1"], "record K1"),
      ["same TOV equations, SLy piecewise-polytrope equation of state"],
      "SLy parameters (Read, Lackey, Owen and Friedman 2009), central density scanned over " + N(V("sly_scan")[0]) + " to " + N(V("sly_scan")[1]) + " MeV/fm³",
      [f"Maximum mass = {N(sly_m,3)} solar masses; R({N(STDM)}) = {N(sly_r,3)} km", rec_line("K1b SLy mass", sly_m, V("sly_published_mmax"), 3)],
      f"{N(V('sly_published_mmax'))} solar masses and {N(V('sly_published_r14'))} km — {SRC('sly_published_mmax')}", f"each within {N(V('tol_percent_5')*100,2)}%",
      "PASS" if j_rel(sly_m, V("sly_published_mmax"), V("tol_percent_5")) == "PASS" and j_rel(sly_r, V("sly_published_r14"), V("tol_percent_5")) == "PASS" else "FAIL",
      "If this solver were wrong, none of the other star results could be trusted.")
mmw, r14w, _ = bag_sequence(bB, V("bag_wall_gemini"))
check("K1c", "C", "Containment-only star, and with Gemini's wall term", (R + ANCH["K1"], "record K1"),
      ["p = (ε − 4B)/3 + wall term; same TOV equations"],
      f"B = {N(bB)} MeV/fm³; wall strength = {N(V('bag_wall_gemini'))}",
      [f"Containment only: M_max = {N(mm,3)}, R({N(STDM)}) = {N(r14,3)} km", f"With wall term: M_max = {N(mmw,3)}, R({N(STDM)}) = {N(r14w,3)} km", rec_line("K1c wall mass", mmw, CL("K1cMw"), 3)],
      f"{N(V('ns_mass_measured'))} ± {N(V('ns_mass_sigma'))} solar masses; {N(V('ns_radius_lo'))} to {N(V('ns_radius_hi'))} km", "heaviest star within 1 standard deviation or above it, and radius in range",
      "PASS" if mm >= V("ns_mass_measured") - V("ns_mass_sigma") and mmw >= V("ns_mass_measured") - V("ns_mass_sigma") and V("ns_radius_lo") <= r14 <= V("ns_radius_hi") and V("ns_radius_lo") <= r14w <= V("ns_radius_hi") else "FAIL")
fj4 = (bB * hbc ** 3) ** 0.25; fj_min4 = (V("bag_B_fj_min") * hbc ** 3) ** 0.25
check("K1d", "C", "The stability rule: is Gemini's B allowed?", (R + ANCH["K1"], "record K1"),
      ["B ≥ B_min  (Farhi and Jaffe 1984), i.e. B¹ᐟ⁴ ≥ about 145 MeV"],
      f"B = {N(bB)} MeV/fm³ → B^(1/4) = {N(fj4,4)} MeV; the bound B_min = {N(V('bag_B_fj_min'))} MeV/fm³ → {N(fj_min4,4)} MeV",
      [f"B = {N(bB)} is {'above' if bB >= V('bag_B_fj_min') else 'below'} the bound."],
      f"B_min = {N(V('bag_B_fj_min'))} MeV/fm³ — {SRC('bag_B_fj_min')}", "B at or above B_min", "PASS" if bB >= V("bag_B_fj_min") else "FAIL",
      "The record says this too: this B is slightly too small, and the fix is a larger B.")
# K2 tidal
cs2_arr = np.gradient(Ps, Es)
dpde_s = lambda e: float(np.interp(e, Es, cs2_arr))
sl2 = np.array([tidal_star(sly_pf, sly_inv, dpde_s, ec, 0.0, psurf=1e-9) for ec in np.linspace(V("sly_scan")[0], V("sly_scan")[1], 120)])
ii = sl2[:, 0].argmax(); rr_ = sl2[:ii + 1]
sly_lam = float(np.exp(np.interp(STDM, rr_[:, 0], np.log(rr_[:, 3]))))
check("K2a", "C", "Solver check for the tidal number Λ (SLy)", (R + ANCH["K2"], "record K2"),
      ["k₂ from y(R) and C = M/R (Hinderer 2008);  Λ = (2/3) k₂ / C⁵"],
      f"SLy equation of state; star of {N(STDM)} solar masses",
      [f"Λ({N(STDM)}) = {N(sly_lam,3)}", rec_line("K2a SLy Lambda", sly_lam, CL("K2a"), 3)],
      f"about {N(V('sly_published_lambda'))} — {SRC('sly_published_lambda')}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(sly_lam, V("sly_published_lambda"), V("tol_percent_5")))
tid_lines = [ln for ln in PUB["neutron_star_tidal_fix"][2].splitlines() if "Lambda(" + str(STDM) + ")" in ln]
for n_, ln in enumerate(tid_lines):
    Bv, wl = re.findall(r"B=([\d.]+) wall=([\d.]+)", ln)[0]; Mx = float(re.search(r"Mmax=([\d.]+)", ln).group(1)); Lm = float(re.search(r"Lambda\(" + re.escape(str(STDM)) + r"\)=(\d+)", ln).group(1))
    ok2 = Mx >= mt and Lm <= V("lambda_limit")
    check(f"K2-{n_ + 1}", "C", f"Tidal number from the published script, run {I(n_ + 1)}", (R + ANCH["K2"], "record K2"),
          ["k\u2082 and \u039b from the TOV and Love-number equations (Hinderer 2008), self-bound surface correction"],
          f"B = {N(float(Bv))} MeV/fm\u00b3, wall strength {N(float(wl))} (printed by neutron_star_tidal_fix)",
          [f"M_max = {N(Mx,3)} solar masses; \u039b({N(STDM)}) = {N(Lm,3)}"], f"M_max at least {N(mt)}; \u039b at most {N(V('lambda_limit'))} \u2014 {SRC('lambda_limit')}",
          "both tests", "PASS" if ok2 else "FAIL",
          "These runs fit with two adjustable numbers (B and the wall strength): superseded in the record by the no-knobs result (stage A3)." if float(wl) > 0 else "")
# K3
mp_ = V("m_proton")
def B_of_R(Rf): return 3 * mp_ / (16 * math.pi * Rf ** 3)
B_full = B_of_R(V("proton_radius_fm")); Bc_lo = B_of_R(V("core_radius_hi")); Bc_hi = B_of_R(V("core_radius_lo"))
def R_of_B(Bv): return (3 * mp_ / (16 * math.pi * Bv)) ** (1 / 3)
R_in_lo, R_in_hi = R_of_B(V("bag_B_window_hi")), R_of_B(V("bag_B_window_lo"))
rms_f = math.sqrt(3 / 5)
check("K3a", "C", "Proton containment balance: B from the proton's whole radius", (R + ANCH["K3"], "record K3"),
      ["M = 4 B V,  V = (4/3)πR³   →   B = 3M/(16πR³)"],
      f"M = {N(mp_,6)} MeV; R = {N(V('proton_radius_fm'))} fm (the whole measured radius)",
      [f"B = {N(B_full,3)} MeV/fm³", rec_line("K3 B whole radius", B_full, CL("K3B"), 3),
       f"For comparison the stability bound is {N(V('bag_B_fj_min'))} MeV/fm³, so this B is {'above' if B_full >= V('bag_B_fj_min') else 'below'} it."],
      f"the record states about {N(CL('K3B'))} MeV/fm³, 'too small' (below the stability bound {N(V('bag_B_fj_min'))})", f"within {N(V('tol_percent_5')*100,2)}% of the record's figure",
      j_rel(B_full, CL("K3B"), V("tol_percent_5")),
      "The record's own formula does not give its stated number, and the result is the opposite side of the stability bound. This is flagged for the authors to reconcile.")
core_rms_lo, core_rms_hi = rms_f * R_in_lo, rms_f * R_in_hi
dev = max(abs(core_rms_lo / V("core_radius_lo") - 1), abs(core_rms_hi / V("core_radius_hi") - 1))
check("K3c", "C", "The quark-core size that gives a B inside the neutron-star window", (R + ANCH["K3"], "record K3"),
      ["R = (3M / (16πB))^(1/3);   rms radius of a uniform sphere = √(3/5) R"],
      f"M = {N(mp_,6)} MeV; B from {N(V('bag_B_window_lo'))} to {N(V('bag_B_window_hi'))} MeV/fm³ (the window)",
      [f"R from {N(R_in_lo,3)} to {N(R_in_hi,3)} fm; as an rms radius {N(core_rms_lo,3)} to {N(core_rms_hi,3)} fm",
       f"The record states a core of {N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm (rms)."],
      f"{N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm — {SRC('core_radius_lo')}", f"both ends within {N(V('tol_percent_5')*100,2)}% of the record's range",
      "PASS" if dev <= V("tol_percent_5") else ("PARTIAL" if dev <= PF * V("tol_percent_5") else "FAIL"),
      "The rms conversion (uniform sphere) is this page's reading of the record's 'rms'. The record's M1 text uses the bag radius itself, which matches the R computed here.")
check("K3b", "C", "Gemini's centroid and its 'hyperfine' number", (R + ANCH["K3"], "record K3"),
      ["centroid = (M_N + M_Δ)/2;  C_ss = gap / 6;  predicted gap = 6 C_ss"],
      f"M_N = {N(V('nucleon_mass_avg'))} MeV, M_Δ = {N(V('delta_mass'))} MeV; measured gap = {N(V('np_delta_gap'))} MeV; C_ss = {N(V('gemini_C_ss'))} MeV",
      [f"Centroid = {N((V('nucleon_mass_avg')+V('delta_mass'))/2,5)} MeV; gap / 6 = {N(V('np_delta_gap')/6,4)} MeV; 6 × C_ss = {N(6*V('gemini_C_ss'),5)} MeV"],
      f"measured gap {N(V('np_delta_gap'))} MeV", "a prediction must not contain its own answer",
      "FAIL" if abs(6 * V("gemini_C_ss") - V("np_delta_gap")) / V("np_delta_gap") < V("tol_percent_1") else "PASS",
      "Circular: C_ss was typed in as the measured gap divided by 6, so the 'prediction' returns the input.")
# K4 other mixings
pairs = [("V_ub", "q_u", "q_t", "ckm_Vub"), ("V_cb", "q_c", "q_t", "ckm_Vcb"), ("V_td", "q_d", "q_b", "ckm_Vtd"), ("V_ts", "q_s", "q_b", "ckm_Vts")]
for nm, a, b, ck in pairs:
    val = math.sqrt(V(a) / V(b)); meas = V(ck)
    check("K4-" + nm[2:], "C", f"The same √(m_i/m_j) pattern for {nm}", (R + ANCH["K4"], "record K4"),
          [f"{nm} ≈ √(m_{a[2]} / m_{b[2]})"], f"m_{a[2]} = {N(V(a))} MeV, m_{b[2]} = {N(V(b))} MeV",
          [f"√(m_{a[2]}/m_{b[2]}) = {N(val,3)}; ratio to measured = {N(val/meas,3)}"],
          f"{nm} = {N(meas)} — {SRC(ck)}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(val, meas, V("tol_percent_5")),
          "The record names V_ub explicitly and says V_cb, V_td and V_ts fail; the mass pairs for those three are the ones the same pattern points to, chosen here." if nm != "V_ub" else "The record: roughly works for V_ub, a little low.")
# K5 / K10 counts
if XLSX.exists():
    parts = read_xlsx(XLSX, "Particles"); pairs_x = read_xlsx(XLSX, "All_Pairs")
    cnt = [len(parts), len(pairs_x), sum(1 for r in pairs_x if (r.get("charge_product") or 0) != 0), sum(1 for r in pairs_x if r.get("strong_link") == 1),
           sum(1 for r in pairs_x if (r.get("W_transition_amplitude") or 0) > 0)]
    st5 = V("xlsx_stated_counts")
    check("K5", "C", "Every-particle-to-every-particle table: the counts", (R + ANCH["K5"], "record K5"),
          ["count rows of the downloadable workbook"], "TRPG_Every_Particle_to_Every_Particle.xlsx, sheets Particles and All_Pairs (read by this script)",
          [f"Particles = {I(cnt[0])}; pairs = {I(cnt[1])}; linked by the electric force = {I(cnt[2])}; by the strong force = {I(cnt[3])}; by weak transitions = {I(cnt[4])}"],
          "the record states " + ", ".join(I(v) for v in st5) + " (in the same order)", "exact counts", "PASS" if cnt == st5 else "FAIL")
    pa = read_xlsx(XLSX, "Particle_vs_Antiparticle")
    fates = [r.get("fate") or "" for r in pa]
    fc = {"opposites": sum(int(r.get("direct_opposites") or 0) for r in pa), "annihilate": sum(1 for f_ in fates if f_.startswith("can annihilate -> neutral")),
          "pions": sum(1 for f_ in fates if "pions" in f_), "W": sum(1 for f_ in fates if "fuse" in f_),
          "flavour": sum(1 for f_ in fates if f_.startswith("flavour-changing")), "scatter": sum(1 for f_ in fates if f_.startswith("no annihilation"))}
    sf5 = V("fate_stated")
    check("K10a", "C", "Particle-versus-antiparticle table: the counts (corrected rules)", (R + ANCH["K10"], "record K10"),
          ["count rows of the downloadable workbook"], "sheet Particle_vs_Antiparticle of the same workbook (read by this script)",
          [f"Direct opposites = {I(fc['opposites'])}; annihilate into neutral carriers = {I(fc['annihilate'])}; into pions = {I(fc['pions'])}; fuse into a W = {I(fc['W'])}; flavour-changing = {I(fc['flavour'])}; scatter only = {I(fc['scatter'])}"],
          "the record states " + ", ".join(f"{k} {I(v)}" for k, v in sf5.items()), "exact counts", "PASS" if fc == sf5 else "FAIL")
    # M15 network (force carriers included)
    P_ = {r["name"]: r for r in parts}
    def graph(use_higgs):
        E = {}
        for r in pairs_x:
            linked = (r.get("charge_product") or 0) != 0 or r.get("strong_link") == 1 or (r.get("W_transition_amplitude") or 0) > 0 or (use_higgs and (r.get("Higgs_coupling_product") or 0) > 0)
            if linked: E.setdefault(r["particle_i"], set()).add(r["particle_j"]); E.setdefault(r["particle_j"], set()).add(r["particle_i"])
        return E
    names = sorted(set(r["particle_i"] for r in pairs_x) | set(r["particle_j"] for r in pairs_x))
    E = graph(True)
    def link(a, b): E.setdefault(a, set()).add(b); E.setdefault(b, set()).add(a)
    for n_, r in P_.items():
        if n_ in ("photon", "gluon", "W", "Z", "Higgs"): continue
        if (r.get("charge") or 0) != 0: link("photon", n_)
        if r.get("colored"): link("gluon", n_)
        if (r.get("T3") or 0) != 0 or (r.get("charge") or 0) != 0: link("Z", n_); link("W", n_)
        if (r.get("mass_MeV") or 0) > 0: link("Higgs", n_)
    link("W", "photon"); link("W", "Z"); link("Higgs", "W"); link("Higgs", "Z")
    dist = []; unreached = 0
    for a in names:
        seen = {a: 0}; fr_ = [a]
        while fr_:
            nx = []
            for u in fr_:
                for v in E.get(u, ()):
                    if v not in seen: seen[v] = seen[u] + 1; nx.append(v)
            fr_ = nx
        for b in names:
            if b > a:
                if b in seen: dist.append(seen[b])
                else: unreached += 1
    nstat = V("network_stated")
    check("M15b", "D", "Degrees of separation in the particle table, with force carriers", (R + ANCH["M15"], "record M15"),
          ["shortest path over links (any direct interaction); average over all pairs"], "the particle workbook plus photon, gluon, W, Z and Higgs links",
          [f"Longest route = {I(max(dist))} steps; average = {N(sum(dist)/len(dist),3)}; pairs never connected = {I(unreached)}"],
          f"the record states longest {I(nstat[0])}, average {N(nstat[1])}", "same longest route and average within 1%",
          "PASS" if max(dist) == nstat[0] and abs(sum(dist) / len(dist) / nstat[1] - 1) <= V("tol_percent_1") and unreached == 0 else "FAIL")
# K6
Hs2 = V("H0") * 1e3 / V("Mpc_m"); rho_c2 = 3 * Hs2 ** 2 / (8 * math.pi * G); Rh = V("c_exact") / Hs2
Mu = (4 * math.pi / 3) * Rh ** 3 * rho_c2; mpk = V("m_proton_kg")
n_part = Mu / mpk
lam_c = V("hbarc_MeVfm") * 1e-15 / mp_
hbar_SI = V("hbarc_MeVfm") * 1e-15 * V("MeV_J") / V("c_exact")
alpha_g = G * mpk ** 2 / (hbar_SI * V("c_exact"))
wein = ((hbar_SI ** 2 * Hs2) / (G * V("c_exact"))) ** (1 / V("weinberg_relation_exponent_den"))
wein_MeV = wein * V("c_exact") ** 2 / V("MeV_J")
check("K6", "C", "Large-number coincidences in the proton's scale view", (R + ANCH["K6"], "record K6"),
      ["N = M_universe / m_p;  M_universe = (4π/3)(c/H₀)³ ρ_c;  (c/H₀)/ƒλ_p;  G m_p²/(ℏc);  m_W = (ℏ²H₀/(Gc))^(1/3)"],
      f"H₀ = {N(V('H0'))} km/s/Mpc, m_p = {N(mpk)} kg",
      [f"Particles in the universe's mass = {N(n_part,3)}", f"Universe across in proton Compton lengths = {N(Rh/(lam_c),3)}", f"Gravity weaker than the quantum scale by {N(1/alpha_g,3)}", f"Weinberg relation gives {N(wein_MeV,3)} MeV",
       rec_line("K6 N particles", n_part, CL("K6N"), 3)],
      f"Eddington's number is about {N(V('eddington_stated'))} (Eddington 1930s, Dirac 1937, Weinberg 1972)", "same order of magnitude as the stated number", "PARTIAL" if abs(math.log10(n_part / V("eddington_stated"))) < 1 else "FAIL",
      "PARTIAL because it is a coincidence, not a derivation, and the record notes every particle's view gives some huge number.")
# K7
rng2 = np.random.default_rng(int(V("random_seed"))); Nr = int(V("random_triples")); lo_, hi_ = V("random_range_log10")
mr = 10 ** rng2.uniform(lo_, hi_, (Nr, 3)); qr = mr.sum(1) / np.sqrt(mr).sum(1) ** 2
tolq = abs(Qlep - Q23); chance = float(np.mean(abs(qr - Q23) <= tolq))
check("K7a", "C", "Koide's relation on the charged leptons, and how special it is", (R + ANCH["K7"], "record K7"),
      ["Q = (m_e + m_μ + m_τ)/(√m_e + √m_μ + √m_τ)² = 2/3"],
      f"m_e, m_μ, m_τ as in stage A4; random test: {I(Nr)} log-uniform mass triples between {N(10**lo_)} and {N(10**hi_)} MeV",
      [f"Q = {N(Qlep,7)} (2/3 = {N(Q23,7)})", f"Chance of a random triple matching as closely = {N(chance,2)}", rec_line("K7 random chance", chance, V("random_chance_stated"), 2)],
      f"2/3 — Koide (1981)", f"Q within {N(V('tol_lepton_mass')*100,2)}% of 2/3", j_rel(Qlep, Q23, V("tol_lepton_mass")))
qk = {"u": V("q_u"), "d": V("q_d"), "s": V("q_s"), "c": V("q_c"), "b": V("q_b"), "t": V("q_t")}
for trip, nm in [("uct", "up-type (u, c, t)"), ("dsb", "down-type (d, s, b)"), ("cbt", "heavy trio (c, b, t)")]:
    qv = Qk(*[qk[x] for x in trip])
    check("K7-" + trip, "C", f"Koide's relation for quarks: {nm}", (R + ANCH["K7"], "record K7"),
          ["Q = Σm / (Σ√m)² = 2/3 ?"], "masses: " + ", ".join(f"{x} = {N(qk[x])} MeV" for x in trip),
          [f"Q = {N(qv,4)}"], "2/3", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(qv, Q23, V("tol_percent_1")),
          "Quarks do not follow it as families; only the heavy trio comes near (record K7)." if trip != "cbt" else "")
check("K7-knobs", "C", "Gemini's mass formula: knobs versus results", (R + ANCH["K7"], "record K7"),
      ["knobs must be fewer than the masses they predict"], "the record's list: " + "; ".join(V("gemini_knob_count_items")),
      [f"Knobs listed = {I(len(V('gemini_knob_count_items')))}; masses predicted = {I(V('gemini_masses_count'))}"],
      "test T5 (fewer knobs than results)", "knobs fewer than results", "PASS" if len(V("gemini_knob_count_items")) < V("gemini_masses_count") else "FAIL")
check("K7-proton", "C", "Brannen scale versus the proton mass", (R + ANCH["K7"], "record K7"),
      ["3 M² ≈ m_p ?"], f"3M² = {N(3*Mt*Mt,5)} MeV; m_p = {N(mp_,6)} MeV",
      [f"Difference = {N((3*Mt*Mt/mp_-1)*100,2)}%"], f"{N(mp_,6)} MeV — {SRC('m_proton')}", f"within {N(V('proton_scale_hint')*100,2)}%",
      "PARTIAL" if abs(3 * Mt * Mt / mp_ - 1) <= V("proton_scale_hint") * 1.1 else "FAIL", "A hint only: unexplained, as the record says.")
# K8 Gemini claims
cent = V("gemini_centroid"); Css = V("gemini_C_ss"); MN, MD = cent - 3 * Css, cent + 3 * Css; mpi = V("pion_mass")
pp_ = math.sqrt((MD ** 2 - (MN + mpi) ** 2) * (MD ** 2 - (MN - mpi) ** 2)) / (2 * MD); EN = math.sqrt(MN ** 2 + pp_ ** 2)
Gam = lambda f: (f ** 2 / (12 * math.pi * mpi ** 2)) * ((EN + MN) / (2 * MD)) * pp_ ** 3
gw = Gam(V("gemini_f")); f_need = V("gemini_f") * math.sqrt(V("delta_width_measured") / gw)
check("K8a", "C", "Gemini's Δ decay width", (R + ANCH["K8"], "record K8"),
      ["Γ = f² p³ (E_N + M_N) / (12π m_π² · 2M_Δ)"],
      f"f = {N(V('gemini_f'))} (typed by Gemini), M_N = {N(MN,5)}, M_Δ = {N(MD,5)}, m_π = {N(mpi)} MeV",
      [f"Width = {N(gw,3)} MeV; Gemini claimed {N(V('gemini_width_claim'))} MeV; the f needed for the measured value is {N(f_need,3)}"],
      f"{N(V('delta_width_measured'))} MeV — {SRC('delta_width_measured')}", f"within {N(V('tol_percent_5')*100,2)}%", j_rel(gw, V("delta_width_measured"), V("tol_percent_5")),
      "The claimed number does not come out of Gemini's own formula.")
Z2 = 1 / V("gemini_Z2_denominator"); gA_f = Z2 * (V("gemini_gA_bare_sq") - (2 / 27) * V("gemini_gA_delta") * V("gemini_gA_bare_sq"))
check("K8b", "C", "Gemini's axial coupling g_A", (R + ANCH["K8"], "record K8"),
      ["g_A = Z₂ (g₀² − (2/27) δ g₀²)"], f"Z₂ = 1/{N(V('gemini_Z2_denominator'))}, g₀² = {N(V('gemini_gA_bare_sq'))}, δ = {N(V('gemini_gA_delta'))}",
      [f"From the stated formula = {N(gA_f,4)}; Gemini's claim = {N(V('gemini_gA_claim'))} (needs a hand-added {N(V('gemini_gA_hand_add'))})"],
      f"{N(V('gA_measured'))} — {SRC('gA_measured')}", f"within {N(V('tolerance_gA')*100,2)}%", j_rel(gA_f, V("gA_measured"), V("tolerance_gA")))
mu_c = [Z2 * v for v in V("gemini_mu_bare")]
check("K8c", "C", "Gemini's proton and neutron magnetic moments", (R + ANCH["K8"], "record K8"),
      ["μ = Z₂ × bare value"], f"Z₂ = {N(Z2,4)}; bare values {N(V('gemini_mu_bare')[0])} and {N(V('gemini_mu_bare')[1])}",
      [f"From the stated factor: {N(mu_c[0],4)} and {N(mu_c[1],4)}; Gemini typed {N(V('gemini_mu_typed')[0])} and {N(V('gemini_mu_typed')[1])}, then added a further hand term (claim {N(V('gemini_mu_claim')[0])}, {N(V('gemini_mu_claim')[1])})"],
      "the record states no measured value here", f"the typed value must equal the formula's value within {N(V('tol_percent_1')*100,2)}%",
      j_rel(V("gemini_mu_typed")[0], mu_c[0], V("tol_percent_1")), "Typed numbers, not computed ones.")
check("K8d", "C", "Pressure inside the proton and the D-term", (R + ANCH["K8"], "record K8"),
      ["von Laue: ∫ p(r) r² dr = 0;  D-term"], f"pressure switch at about {N(V('jlab_pressure_switch_fm'))} fm (Jefferson Lab, 2018); D-term bag model {N(V('dterm_bag'))} versus Jefferson Lab estimate {N(V('dterm_jlab'))}",
      [f"Quark-core radius {N(V('core_radius_lo'))} to {N(V('core_radius_hi'))} fm against the switch at {N(V('jlab_pressure_switch_fm'))} fm; D-term ratio {N(V('dterm_jlab')/V('dterm_bag'),3)}"],
      "extractions of the D-term still differ between groups", "not decided", "NOT TESTABLE YET", "Open comparison, not a fail (record K8): the measured value itself is disputed.")
# K9 compactness
Csun = V("compactness_sun"); Cwd = 2 * V("sirius_b_mass") * V("Msun_m") / V("sirius_b_radius_m")
Rns = r14 * 1e3; Cns = 2 * STDM * V("Msun_m") / Rns; Cb = V("buchdahl_limit_num") / V("buchdahl_limit_den")
check("K9", "C", "The containment ladder: compactness C = 2GM/(Rc²)", (R + ANCH["K9"], "record K9"),
      ["C = 2GM/(Rc²);  Buchdahl: C < 8/9 for any static star"],
      f"Sun C = {N(Csun)}; white dwarf (M = {N(V('sirius_b_mass'))} solar masses, R = {N(V('sirius_b_radius_m'))} m); neutron star R = {N(r14,3)} km (stage A3b)",
      [f"White dwarf C = {N(Cwd,3)}; neutron star C = {N(Cns,3)}; Buchdahl limit {fr(V('buchdahl_limit_num'), V('buchdahl_limit_den'))} = {N(Cb,4)}; horizon C = {N(V("horizon_compactness"))}"],
      "Buchdahl (1959): no static star passes the limit without collapsing", "neutron star below the Buchdahl limit and above the white dwarf", "PASS" if Cwd < Cns < Cb else "FAIL")
# K10 Fermi muon lifetime
check("K10b", "C", "Muon lifetime from Fermi theory", (R + ANCH["K10"], "record K10"),
      ["τ_μ = 192π³ℏ / (G_F² m_μ⁵)"], f"G_F = {N(GF)} MeV⁻², m_μ = {N(mmu,7)} MeV",
      [f"τ_μ = {N(tau_mu_c,5)} s", rec_line("K10 muon", tau_mu_c * 1e6, CL("K10mu"), 4)],
      f"{N(V('tau_mu_s'),7)} s — {SRC('tau_mu_s')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(tau_mu_c, V("tau_mu_s"), V("tol_percent_1")),
      "The record gives this in microseconds; the comparison above converts it.")

# B8 (new count)
trpg_items = V("inputs_trpg_items"); sm_items = V("inputs_sm_items")
check("B8", "B", "Count of measured inputs", (R + ANCH["P1"], "record P1"),
      ["inputs (TRPG) versus inputs (Standard Model with neutrino masses, gravity and dark energy)"],
      "the record's own lists (equation page and Q2); B and the glueball mass are counted as consequences of the strong force",
      [f"TRPG list holds {I(len(trpg_items))} items; Standard Model list holds {I(len(sm_items))}; saving = {I(len(sm_items) - len(trpg_items))}"],
      f"the record states about {I(CL('inputs_trpg'))} against about {I(CL('inputs_sm'))}", "same counts as the record",
      "PASS" if len(trpg_items) == CL("inputs_trpg") and len(sm_items) == CL("inputs_sm") else "FAIL")
# ---- rev3 B9: B set by the proton (published script, run as published)
pt = PUB["proton_sets_B"][2]
def _first(pat, txt, grp=1, nth=0):
    m_ = re.findall(pat, txt); return float(m_[nth][grp - 1] if isinstance(m_[nth], tuple) else m_[nth])
try:
    B1975 = _first(r"B = ([\d.]+) MeV/fm\^3 \(B\^1/4", pt, nth=0); B_p = _first(r"B = ([\d.]+) MeV/fm\^3 \(B\^1/4", pt, nth=1)
    Rn = _first(r"proton bag radius R = ([\d.]+) fm", pt); rms_p = _first(r"quark-core rms = ([\d.]+) fm", pt)
    cl_lo = _first(r"coupling x0.95: B = ([\d.]+)", pt); cl_hi = _first(r"coupling x1.05: B = ([\d.]+)", pt); parsed_ok = True
except Exception:
    parsed_ok = False; B1975 = B_p = Rn = rms_p = cl_lo = cl_hi = float("nan")
check("B9a", "B", "Containment depth B set by the proton and the Delta, nothing adjusted", (R + ANCH["Q3"], "record Q3"),
      ["fit B and \u03b1_c to M_N and M_\u0394;  E = (4\u03c0/3) B R\u00b3 + ...;  centre-of-mass correction M = \u221a(E\u00b2 \u2212 \u27e8P\u00b2\u27e9);  pion cloud (cloudy bag)"],
      "the published proton_sets_B script, run with python -I (its real output is in stage F2)",
      [f"B = {N(B_p,3)} MeV/fm\u00b3, band {N(cl_lo,3)} to {N(cl_hi,3)} for a \u00b1 coupling change; bag radius {N(Rn,3)} fm; quark-core rms {N(rms_p,3)} fm",
       rec_line("B9 B", B_p, CL("Q3B"), 3), rec_line("B9 band low", cl_lo, CL("Q3Blo"), 3), rec_line("B9 band high", cl_hi, CL("Q3Bhi"), 3), rec_line("B9 rms", rms_p, CL("Q3core"), 3)],
      f"no direct measurement of B; the fit reproduces the measured proton and Delta masses", f"within {N(V('tol_percent_1')*100,2)}% of the record's central value", j_rel(B_p, CL("Q3B"), V("tol_percent_1")) if parsed_ok else "FAIL")
check("B9b", "B", "Check of the same code on the earlier bag fit", (R + ANCH["Q3"], "record Q3"), ["same code with the earlier settings"],
      "settings as printed by the script (no centre-of-mass correction, no cloud)",
      [f"B = {N(B1975,3)} MeV/fm\u00b3", rec_line("B9 earlier fit", B1975, CL("Q3check1975"), 3)],
      f"published about {N(CL('Q3pub'))} MeV/fm\u00b3 (DeGrand, Jaffe, Johnson and Kiskis, the earlier bag fit)",
      f"within {N(V('tol_percent_5')*100,2)}%", j_rel(B1975, CL("Q3pub"), V("tol_percent_5")) if parsed_ok else "FAIL")
# K3d: does the new Q3 radius agree with what K3 says?
check("K3d", "C", "The proton core size: K3's statement against the Q3 run", (R + ANCH["K3"], "record K3"),
      ["compare the quark-core rms of the Q3 run with the range stated in K3"],
      f"Q3 run: rms {N(rms_p,3)} fm; K3 states {N(CL('K3rms_lo'))} to {N(CL('K3rms_hi'))} fm",
      [f"Difference from the upper end of K3's range: {N((rms_p/CL('K3rms_hi')-1)*100,2)}%"], f"{N(CL('K3rms_lo'))} to {N(CL('K3rms_hi'))} fm \u2014 record K3",
      f"within {N(V('tol_percent_5')*100,2)}% of the range", "PASS" if rms_p <= CL("K3rms_hi") * (1 + V("tol_percent_5")) else "PARTIAL",
      "The later Q3 run is the better source; K3's own B figures (stages K3a and K3c) are the older text.")
# Q5: the quark link-pattern idea
qt = PUB["quark_link_pattern"][2]
q5 = re.findall(r"r\^2 = ([\d.]+) \+/- ([\d.]+);\s+simple value ([\d./]+) = ([\d.]+) is (-?[\d.]+) sigma", qt)
q5ok = len(q5) == 2
upr, dnr = (float(q5[0][0]), float(q5[1][0])) if q5ok else (float("nan"), float("nan"))
upsig, dnsig = (float(q5[0][4]), float(q5[1][4])) if q5ok else (float("nan"), float("nan"))
check("Q5", "D", "Quark sectors in the leptons' relational form: the tested idea fails", (R + ANCH["Q5"], "record Q5"),
      ["r\u00b2 = (3Q \u2212 1)/2;  idea tested: r\u00b2 = 3/6, 4/6, 5/6 for leptons, down-type, up-type"],
      "the published quark_link_pattern script, run with python -I (real output in stage F3)",
      [f"up-type r\u00b2 = {N(upr,4)}, which is {N(abs(upsig),3)} standard deviations from the simple value; down-type r\u00b2 = {N(dnr,4)}, {N(abs(dnsig),3)} standard deviations away",
       rec_line("Q5 up", upr, CL("Q5up"), 3), rec_line("Q5 down", dnr, CL("Q5down"), 3)],
      "the measured quark masses (PDG 2024) give these r\u00b2 values", "the simple fractions within 2 standard deviations (look-elsewhere chance in the script's last line)",
      "PASS" if q5ok and abs(upsig) <= 2 and abs(dnsig) <= 2 else "FAIL", "The idea fails and the record drops it. Quark masses stay measured inputs.")

# ================================================================================================ SECTION D: chapter M and P extras
SECTIONS["D"] = "Chapters M and P: the rest"
# M5 wave grid
ng = int(V("wave_n_grid")); xg = np.linspace(0, 2 * np.pi, ng); X, Yg = np.meshgrid(xg, xg)
zx, zy = -np.sin(X), -np.sin(Yg); zxx, zyy = -np.cos(X), -np.cos(Yg)
Kc = (zxx * zyy) / (1 + zx ** 2 + zy ** 2) ** 2; dA = np.sqrt(1 + zx ** 2 + zy ** 2) * (xg[1] - xg[0]) ** 2
gb = float((Kc * dA)[:-1, :-1].sum()); fpos = float((dA * (Kc > 0)).sum() / dA.sum()); fneg = float((dA * (Kc < 0)).sum() / dA.sum())
# critical point classification at the four named points
def classify(x0, y0):
    hxx, hyy = -math.cos(x0), -math.cos(y0); det = hxx * hyy
    return "saddle" if det < 0 else ("peak" if hxx < 0 else "trough")
cps = [classify(0, 0), classify(math.pi, math.pi), classify(math.pi, 0), classify(0, math.pi)]
morse = cps.count("peak") - cps.count("saddle") + cps.count("trough")
check("M5", "D", "Crossing waves: peaks, saddles, troughs and flatness", (R + ANCH["M5"], "record M5"),
      ["z = cos x + cos y;  K = z_xx z_yy / (1 + z_x² + z_y²)²;  ∫ K dA = 0 over one cell (Gauss–Bonnet);  peaks − saddles + troughs = 0 (Morse)"],
      f"grid of {I(ng)}×{I(ng)} points over one cell; critical points at (0,0), (π,π), (π,0), (0,π)",
      [f"Types found: {', '.join(cps)}; peaks − saddles + troughs = {N(morse,3)}", f"Area curved like a sphere = {N(fpos,3)}, like a saddle = {N(fneg,3)}", f"Total curvature over the cell = {N(gb,3)}"],
      f"flat on average: Planck collaboration (2018) finds curvature {N(V('Omega_k'))} ± {N(V('Omega_k_sigma'))}", f"total curvature below {N(V('gauss_bonnet_tol'))} and Morse sum zero",
      "PASS" if abs(gb) < V("gauss_bonnet_tol") and morse == 0 else "FAIL")
# M6 waterfall
def Qsg(vals, signs): rr = [s * math.sqrt(v) for v, s in zip(vals, signs)]; return sum(vals) / sum(rr) ** 2
def solve_light(h2, h1, sgn):
    f = lambda x: Qsg([x, h2, h1], [sgn, 1, 1]) - Q23
    xs = [10 ** (k_ / 200) for k_ in range(-600, 1200)]; out = []
    for a, b in zip(xs, xs[1:]):
        if f(a) * f(b) < 0: out.append(brentq(f, a, b))
    return out
Qcbt = Qsg([qk["c"], qk["b"], qk["t"]], [1, 1, 1]); Qscb = Qsg([qk["s"], qk["c"], qk["b"]], [-1, 1, 1])
chs = solve_light(qk["b"], qk["t"], 1); cc_ = min(chs, key=lambda v: abs(v - qk["c"])); sst = solve_light(cc_, qk["b"], -1)
ss_ = min(sst, key=lambda v: abs(v - qk["s"])) if sst else float("nan")
wst = V("waterfall_stated")
check("M6a", "D", "Quark masses from linked Koide triples: charm", (R + ANCH["M6"], "record M6"),
      ["Q(c, b, t) = 2/3 solved for the charm mass, starting from b and t only"], f"m_b = {N(qk['b'])} MeV, m_t = {N(qk['t'])} MeV",
      [f"Q(c,b,t) = {N(Qcbt,4)}; predicted charm = {N(cc_,5)} MeV; measured charm = {N(qk['c'])} MeV ({N((cc_/qk['c']-1)*100,2)}% off)", rec_line("M6 charm", cc_, wst["charm_pred"], 4)],
      f"{N(qk['c'])} MeV — {SRC('q_c')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(cc_, qk["c"], V("tol_percent_1")),
      "Quark masses run with energy and these inputs mix scales, so the record calls this a hint.")
check("M6b", "D", "Quark masses from linked Koide triples: strange", (R + ANCH["M6"], "record M6"),
      ["Q(−s, c, b) = 2/3 solved for the strange mass"], f"m_c (from the step above) = {N(cc_,5)} MeV, m_b = {N(qk['b'])} MeV",
      [f"Q(s,c,b) with the measured values = {N(Qscb,4)}; predicted strange = {N(ss_,4)} MeV; measured {N(qk['s'])} MeV ({N((ss_/qk['s']-1)*100,2)}% off)", rec_line("M6 strange", ss_, wst["strange_pred"], 4)],
      f"{N(qk['s'])} MeV — {SRC('q_s')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(ss_, qk["s"], V("tol_percent_1")))
# M7 resistance
rM = V("Mars_orbit_in_AU") * V("AU_m"); rE = V("AU_m")
shap = (2 * G * MSUN / C ** 3) * math.log(4 * rE * rM / RSUN ** 2) * 1e6
rt = 2 * shap; lo_rt, hi_rt = V("shapiro_roundtrip_range_us")
check("M7a", "D", "Resistance as refractive index: radar delay grazing the Sun", (R + ANCH["M7"], "record M7"),
      ["Δt = (2GM/c³) ln(4 r_E r_M / b²)  (one way)"], f"r_E = {N(rE)} m, r_M = {N(V('Mars_orbit_in_AU'))} AU, b = R_sun = {N(RSUN)} m",
      [f"One way = {N(shap,3)} microseconds; round trip = {N(rt,3)} microseconds", rec_line("M7a one way", shap, CL("M7a"), 3)],
      f"{N(lo_rt)} to {N(hi_rt)} microseconds round trip — {SRC('shapiro_roundtrip_range_us')}", "round trip inside the quoted range", j_range(rt, lo_rt, hi_rt))
# M7 bands & localization
def tr(k_, n1=1.0, n2=1.5, a=0.5, bb=0.5):
    def layer(n, d): q = k_ * n; return np.array([[math.cos(q * d), math.sin(q * d) / n], [-n * math.sin(q * d), math.cos(q * d)]])
    T = layer(n1, a) @ layer(n2, bb); return abs(np.trace(T)) / 2
ks = np.linspace(0.01, 12, 4000); allowed = np.array([tr(k_) <= 1 for k_ in ks]); edges = np.where(np.diff(allowed.astype(int)) != 0)[0]
gaps = [(ks[edges[i]], ks[edges[i + 1]]) for i in range(0, len(edges) - 1, 2) if not allowed[edges[i] + 1]]
check("M7b", "D", "Periodic resistance makes bands and gaps", (R + ANCH["M7"], "record M7"),
      ["|Tr T|/2 ≤ 1 allowed;  > 1 forbidden (transfer matrix of alternating layers)"], "alternating layers of two refractive indices, wave number scanned from small to large",
      [f"Forbidden gaps found: {', '.join(f'{N(a_,3)} to {N(b_,3)}' for a_, b_ in gaps[:6])}"],
      "Bloch (1928); Kronig and Penney (1931): periodic media have forbidden gaps", f"at least {I(V('bands_min_gaps'))} gaps", "PASS" if len(gaps) >= V("bands_min_gaps") else "FAIL")
rng3 = np.random.default_rng(0)
def T_random(Nl, k_=3.0, spread=0.3):
    T = np.eye(2)
    for _ in range(Nl):
        n = 1 + spread * rng3.random(); d = 0.5; q = k_ * n
        T = np.array([[math.cos(q * d), math.sin(q * d) / n], [-n * math.sin(q * d), math.cos(q * d)]]) @ T
    t = 2 / abs(T[0, 0] + T[1, 1] + 1j * (T[0, 1] * 1 - T[1, 0] / 1)); return t * t
trans = [(Nl, float(np.mean([T_random(Nl) for _ in range(20)]))) for Nl in V("resistance_layers")]
check("M7c", "D", "Random resistance traps waves (Anderson localization)", (R + ANCH["M7"], "record M7"),
      ["transmitted fraction through N random layers falls exponentially with N"], "layer counts " + ", ".join(I(n) for n, _ in trans),
      ["Transmitted fraction: " + "; ".join(f"{I(n)} layers → {N(t_,3)}" for n, t_ in trans)],
      "Anderson (1958)", f"fall by at least a factor {I(V('localization_factor'))} from the shortest to the longest stack", "PASS" if trans[0][1] / trans[-1][1] >= V("localization_factor") else "FAIL")
# M8 compactness of Hubble sphere
Ch = 2 * G * Mu / (Rh * V("c_exact") ** 2)
check("M8", "D", "The universe on the containment ladder", (R + ANCH["M8"], "record M8"),
      ["C = 2GM/(Rc²) with R = c/H₀ and M = (4π/3)R³ρ_c"], f"R = {N(Rh,3)} m, M = {N(Mu,3)} kg",
      [f"C = {N(Ch,5)}"], "identity for a flat universe", f"equal to 1 within {N(V('tol_lepton_mass')*100,2)}%", j_rel(Ch, 1.0, V("tol_lepton_mass")),
      "An identity of flat cosmology: consistent with the idea, not a proof of it.")
# M9 factor 2
full = 4 * G * MSUN / (RSUN * C * C) * ASEC
cases = [("scalar resistance (Nordström 1913)", None), ("time-only resistance (Newton, Soldner 1801)", 0.0), ("time and space in equal shares (Einstein 1915)", 1.0)]
for n_, (nm, gam) in enumerate(cases):
    dv = 0.0 if gam is None else (1 + gam) / 2 * full
    check(f"M9-{'abc'[n_]}", "D", f"Light bending with {nm}", (R + ANCH["M9"], "record M9"),
          ["δ = (1 + γ)/2 × 4GM/(Rc²)"], ("scalar: no coupling to light, so zero" if gam is None else f"γ = {N(gam)}") + f"; 4GM/(Rc²) = {N(full,4)} arcsec",
          [f"Deflection = {N(dv,4)} arcsec"], f"{N(V('bending_measured'))} arcsec — {SRC('bending_measured')}", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(dv, V("bending_measured"), V("tol_percent_1")),
          "Ruled out by the data." if j_rel(dv, V("bending_measured"), V("tol_percent_1")) == "FAIL" else f"The only case that survives. Cassini, as the record cites it: \u03b3 \u2212 1 = {N(V('gamma_minus_1'))} \u00b1 {N(V('gamma_minus_1_sigma'))}.")
# M10 reach
reachP = V("hc_eVm") / V("p_ripple_eV"); reachG = V("hc_eVm") / V("graviton_limit_eV"); need = V("gravity_reach_needed_m")
check("M10a", "D", "The containment field cannot be gravity: its reach", (R + ANCH["M10"], "record M10"),
      ["reach ≈ hc / (mc²)"], f"P ripple mass = {N(V('p_ripple_eV'))} eV (the glueball); graviton limit = {N(V('graviton_limit_eV'))} eV (LIGO-Virgo-KAGRA 2021)",
      [f"Reach of P = {N(reachP,3)} m ({N(reachP*1e15,3)} fm); reach allowed for gravity = at least {N(reachG,3)} m"],
      f"gravity acts over cosmic distances: {SRC('graviton_limit_eV')}", f"for P to be gravity its reach would have to be at least {N(need)} m", "PASS" if reachP < need else "FAIL",
      "PASS here means the check correctly shows P is not gravity.")
check("M10b", "D", "All layers must feel gravity alike (MICROSCOPE)", (R + ANCH["M10"], "record M10"),
      ["Δa/a < measured bound for every pair of materials"], "TRPG layer-dependent gravity has not been computed",
      ["No TRPG number exists yet to compare with the measured bound."], f"MICROSCOPE team (2022) as the record cites", "a TRPG value of the layer difference", "NOT TESTABLE YET")
# M11 expansion
sh = lambda a: ((V("Omega_m") / a ** 3) / (Hf(a) / V("H0")) ** 2, (V("Omega_r") / a ** 4) / (Hf(a) / V("H0")) ** 2, (OLs) / (Hf(a) / V("H0")) ** 2)
a_list = V("expansion_sizes")
lines = []
for a in a_list:
    m_, r2, l_ = sh(a); lines.append(f"size {N(a,1)} of today: matter {N(m_*100,3)}%, radiation {N(r2*100,3)}%, dark energy {N(l_*100,3)}%")
tot_ok = all(abs(sum(sh(a)) - 1) < V("tol_sum_one") for a in a_list)
check("M11", "D", "Always full: the contents add to the whole at every size", (R + ANCH["M11"], "record M11"),
      ["Ω_m + Ω_r + Ω_Λ = 1;  shares at size a: Ω_i(a) / E(a)²"], f"Ω_m = {N(V('Omega_m'))}, Ω_r = {N(V('Omega_r'))}, Ω_Λ = {N(OLs,5)}",
      lines + [f"Shares sum to {I(V('percent'))}% at each size: {'yes' if tot_ok else 'no'}"],
      f"Planck 2018: curvature {N(V('Omega_k'))} ± {N(V('Omega_k_sigma'))}", f"sum equal to 1 within {N(V('tol_sum_one'))} and curvature consistent with zero", "PASS" if tot_ok and abs(V("Omega_k")) <= V("Omega_k_sigma") * 1.0 + 1e-12 + 0.5 * V("Omega_k_sigma") else "FAIL",
      "The curvature term is the only 'empty' share left over; Planck's value is consistent with zero.")
# M12 uncertainty
gp = V("gauss_pulse"); x = np.linspace(-gp["extent"], gp["extent"], gp["points"]); s0 = gp["width"]; psi = np.exp(-x ** 2 / (4 * s0 * s0)); pr = abs(psi) ** 2; pr /= np.trapezoid(pr, x)
dx = math.sqrt(np.trapezoid(x ** 2 * pr, x)); kk = np.fft.fftshift(np.fft.fftfreq(len(x), x[1] - x[0])) * 2 * np.pi
Pk = abs(np.fft.fftshift(np.fft.fft(psi))) ** 2; Pk /= np.trapezoid(Pk, kk); dk = math.sqrt(np.trapezoid(kk ** 2 * Pk, kk))
check("M12", "D", "Uncertainty is wave math", (R + ANCH["M12"], "record M12"),
      ["Δx Δk ≥ 1/2  (equality for a Gaussian pulse)"], f"Gaussian pulse of width {N(s0)} on {N(gp['points'])} points",
      [f"Δx Δk = {N(dx*dk,5)}"], "Heisenberg (1927): the minimum is one half", f"within {N(V('tol_lepton_mass')*100,2)}% of one half", j_rel(dx * dk, 0.5, V("tol_lepton_mass")))
# M14 entanglement density
lp = V("planck_length_m"); dens = 1 / (4 * lp ** 2)
check("M14", "D", "Entanglement per square metre of space", (R + ANCH["M14"], "record M14"),
      ["one unit per 4 Planck areas:  1 / (4 ℓ_P²)"], f"ℓ_P = {N(lp)} m",
      [f"Density = {N(dens,3)} per m²", rec_line("M14 density", dens, V("entangle_density_stated"), 3)],
      f"the record states about {N(V('entangle_density_stated'))} per square metre", f"within a factor {I(V('factor_three_tol'))} of the stated figure",
      "PASS" if (1 / V("factor_three_tol")) <= dens / V("entangle_density_stated") <= V("factor_three_tol") else "FAIL",
      f"The record's figure of {N(V('entangle_density_stated'))} per square metre is one power of ten below what the arithmetic gives. The weak-gravity argument does not depend on it; the figure in the record is a slip.")
# M15 small world
sw = V("small_world_cases"); swS = V("small_world_stated")
lines = []; okw = True
for (Nn, kn), st in zip(sw, swS):
    v_ = math.log(Nn) / math.log(kn); lines.append(f"N = {N(Nn,2)}, {N(kn,3)} links each → {N(v_,3)} steps (record: {N(st,3)})"); okw &= abs(v_ / st - 1) <= V("tol_percent_5")
check("M15a", "D", "Small-world estimate of steps between any two things", (R + ANCH["M15"], "record M15"),
      ["steps ≈ ln N / ln k"], f"people ({N(sw[0][0],2)} with {N(sw[0][1],3)} close contacts) and the universe's particles ({N(sw[1][0],2)} with {N(sw[1][1],3)} links each, and with {N(sw[2][1],3)} links each)", lines,
      "Watts and Strogatz (1998)", f"within {N(V('tol_percent_5')*100,2)}% of the record's figures", "PASS" if okw else "FAIL")
# M16 area law numbers
check("M16", "D", "Area law and the strength of gravity", (R + ANCH["M16"], "record M16"),
      ["same test as stage B1a"], "see stage B1a", [f"Boundary-fit error {N(rb,3)} versus volume-fit error {N(rv,3)}"],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993)", f"ratio at least {N(V('area_ratio_min'))}", "PASS" if rv / rb >= V("area_ratio_min") else "FAIL")
# M18 balance angles
def Qang(v): Qv = sum(v) / sum(math.sqrt(t) for t in v) ** 2; return Qv, math.degrees(math.acos(math.sqrt(1 / (3 * Qv))))
a_l = Qang([me, mmu, mtau])[1]; a_d = Qang([qk["d"], qk["s"], qk["b"]])[1]; a_u = Qang([qk["u"], qk["c"], qk["t"]])[1]
qas = V("quark_angle_stated")
check("M18a", "D", "The balance angle: leptons", (R + ANCH["M18"], "record M18"),
      ["θ = arccos √(1/(3Q))  (45° when Q = 2/3)"], "charged-lepton masses as in stage A4", [f"Angle = {N(a_l,4)} degrees", rec_line("M18 lepton angle", a_l, qas["lepton"], 4)],
      "45° — Koide (1981); Foot (1994)", f"within {N(V('tol_percent_1')*100,2)}%", j_rel(a_l, math.degrees(math.acos(math.sqrt(1 / (3 * Q23)))), V("tol_percent_1")))
a_ref = math.degrees(math.acos(math.sqrt(1 / (3 * Q23))))
for nm, av_ in (("down-type quarks", a_d), ("up-type quarks", a_u)):
    check("M18-" + nm[:2], "D", f"The balance angle: {nm}", (R + ANCH["M18"], "record M18"),
          ["same formula on the quoted masses"], "quark masses as listed in the inputs", [f"Angle = {N(av_,4)} degrees; lepton angle = {N(a_ref,4)} degrees", rec_line("M18 quark angle " + nm, av_, qas["down" if nm.startswith("down") else "up"], 4)],
          f"the lepton value {N(a_ref,4)} degrees", f"within {N(V('tol_percent_1')*100,2)}% of the lepton angle", j_rel(av_, a_ref, V("tol_percent_1")),
          "Quarks do not sit at the lepton balance. The record calls this unsolved, as in all of physics.")
# P2 quark running
MZ = V("mu_quark_run")["MZ"]; aZ = V("mu_quark_run")["alpha_s_MZ"] / math.pi; mc_, mb_, mt_ = V("mu_quark_run")["mc"], V("mu_quark_run")["mb"], V("mu_quark_run")["mt"]
def b_(nf): return (11 - 2 * nf / 3) / 4, (102 - 38 * nf / 3) / 16
def g_(nf): return 1.0, (202 / 3 - 20 * nf / 9) / 16
def nf_at(mu): return 6 if mu > mt_ else 5 if mu > mb_ else 4 if mu > mc_ else 3
def run(a0, mu0, mu1, lnm=0.0):
    th = [mc_, mb_, mt_]; pts = sorted([mu0, mu1] + [t for t in th if min(mu0, mu1) < t < max(mu0, mu1)], reverse=(mu1 < mu0)); a, l = a0, lnm
    for x0, y0 in zip(pts, pts[1:]):
        nf = nf_at(math.sqrt(x0 * y0)); b0, b1 = b_(nf); g0, g1 = g_(nf)
        f = lambda t, s: [-b0 * s[0] ** 2 - b1 * s[0] ** 3, -g0 * s[0] - g1 * s[0] ** 2]
        a, l = solve_ivp(f, [math.log(x0 * x0), math.log(y0 * y0)], [a, l], rtol=1e-10, atol=1e-12).y[:, -1]
    return a, l
a_at = lambda mu: run(aZ, MZ, mu)[0]
m_to = lambda m0, mu0, mu1: m0 * math.exp(run(a_at(mu0), mu0, mu1)[1])
qi = V("quark_run_inputs"); lines = []; angs = []
for mu in V("quark_run_scales"):
    m = {k_: m_to(m0, mu0, mu) for k_, (m0, mu0) in qi.items()}
    up = Qang([m["u"], m["c"], m["t"]]); dn = Qang([m["d"], m["s"], m["b"]]); angs.append((mu, up[1], dn[1]))
    lines.append(f"at {N(mu,5)} GeV: up-type angle {N(up[1],4)}°, down-type angle {N(dn[1],4)}°")
spread_up = max(a[1] for a in angs) - min(a[1] for a in angs); spread_dn = max(a[2] for a in angs) - min(a[2] for a in angs)
rs = V("quark_angle_run_stated"); at91 = [a for a in angs if abs(a[0] - MZ) < 1e-6][0]
check("P2a", "D", "Quark balance angles are the same at every energy", (R + ANCH["P2"], "record P2"),
      ["two-loop running: dα/d lnμ² = −b₀α² − b₁α³;  d ln m/d lnμ² = −γ₀α − γ₁α²"],
      f"α_s(M_Z) = {N(V('mu_quark_run')['alpha_s_MZ'])}; quark masses at their quoted scales ({SRC('quark_run_inputs')}); energies {', '.join(N(v,5) for v in V('quark_run_scales'))} GeV",
      lines + [f"Spread of the up-type angle across energies = {N(spread_up,3)}°; down-type = {N(spread_dn,3)}°", rec_line("P2 up angle at M_Z", at91[1], rs["up"], 4), rec_line("P2 down angle at M_Z", at91[2], rs["down"], 4)],
      "the angle should not change with energy because QCD rescales all quark masses by one factor (record P2)", f"spread below {N(V('tol_percent_1'))} degrees",
      "PASS" if spread_up < V("tol_percent_1") and spread_dn < V("tol_percent_1") else "FAIL")
check("P2b", "D", "Do the quarks sit at the lepton balance once run to one energy?", (R + ANCH["P2"], "record P2"),
      ["compare the angles with the lepton angle"], f"angles at {N(MZ)} GeV: up-type {N(at91[1],4)}°, down-type {N(at91[2],4)}°; leptons {N(a_ref,4)}°",
      [f"Up-type is {N(at91[1]-a_ref,3)}° and down-type {N(at91[2]-a_ref,3)}° away from the leptons."], f"the lepton value {N(a_ref,4)} degrees",
      f"within {N(V('tol_percent_1')*100,2)}% of the lepton angle", "PASS" if abs(at91[1] / a_ref - 1) <= V("tol_percent_1") and abs(at91[2] / a_ref - 1) <= V("tol_percent_1") else "FAIL",
      "They do not. The record states this too: it is now a sharp target for TRPG's curvature-phase idea, not a result.")
# P4 gravity tightened
Lt = np.array(V("area_blocks_tidy"))
S_cross = ent(Ng, mg, True, Lt, V("crossing_amplitude"))
eb = float(np.sqrt(np.mean((np.polyval(np.polyfit(4 * Lt, S_cross, 1), 4 * Lt) - S_cross) ** 2))); ev_ = float(np.sqrt(np.mean((np.polyval(np.polyfit(Lt * Lt, S_cross, 1), Lt * Lt) - S_cross) ** 2)))
check("P4a", "D", "The area law in the crossing-wave medium", (R + ANCH["P4"], "record P4"),
      ["couplings k = 1 + a cos x cos y;  fit S to boundary and to volume"], f"a = {N(V('crossing_amplitude'))}; blocks of side {', '.join(I(l) for l in Lt)}",
      [f"Boundary-fit error {N(eb,3)}; volume-fit error {N(ev_,3)}; ratio {N(ev_/eb,3)}", rec_line("P4 ratio", ev_ / eb, CL("P4"), 3)],
      "Bombelli, Koul, Lee and Sorkin (1986); Srednicki (1993)", f"boundary fit better by at least a factor {I(V('area_ratio_min_crossing'))}", "PASS" if ev_ / eb >= V("area_ratio_min_crossing") else "FAIL")
coefs = []
for m in V("layer_masses"):
    Sm_ = ent(Ng, m, False, Lt); coefs.append((m, float(np.polyfit(4 * Lt, Sm_, 1)[0])))
check("P4b", "D", "Every layer contributes to G, and G acts on all of them alike", (R + ANCH["P4"], "record P4"),
      ["area coefficient c(m) for fields of different mass; G is the single total of the contributions (Jacobson 1994; Susskind and Uglum 1994)"],
      "layer masses " + ", ".join(N(m) for m, _ in coefs), ["Coefficient per boundary unit: " + "; ".join(f"mass {N(m)} → {N(c_,3)}" for m, c_ in coefs)],
      "no measured value: this stage only shows layers differ in what they give to G", "coefficients decrease with mass", "PASS" if all(coefs[i][1] > coefs[i + 1][1] for i in range(len(coefs) - 1)) else "FAIL",
      "It does not test the equivalence principle itself (MICROSCOPE bound, stage M10b).")
# P5 twist
def delta_of(a, b, c_): rr = np.sqrt([a, b, c_]); M_ = rr.mean(); return math.acos((rr[2] / M_ - 1) / math.sqrt(2))
d0 = delta_of(me, mmu, mtau); rng4 = np.random.default_rng(int(V("twist_seed")))
ds = [delta_of(me, mmu + rng4.normal(0, V("sigma_m_mu")), mtau + rng4.normal(0, V("sigma_m_tau"))) for _ in range(int(V("twist_samples")))]
sig = float(np.std(ds)); nsig = (d0 - D29) / sig; ts = V("twist_stated")
check("P5", "D", "The twist δ = 2/9 from the data", (R + ANCH["P5"], "record P5"),
      ["δ = arccos((√m_τ / M − 1)/√2),  M = mean of the three square roots;  compare 2/9"],
      f"masses and uncertainties as in the inputs; {I(V('twist_samples'))} random draws",
      [f"δ from data = {N(d0,6)} ± {N(sig,2)} rad; 2/9 = {N(D29,6)}; difference = {N(nsig,2)} standard deviations", rec_line("P5 twist", d0, ts[0], 6)],
      f"2/9 — Brannen (2006)", f"within {I(V('twist_sigma_limit'))} standard deviations", "PASS" if abs(nsig) <= V("twist_sigma_limit") else "FAIL",
      "That 2/9 equals (2/3)/3 is the project's conjecture, not derived (record P5).")

# ================================================================================================ SECTION E: not recomputed here
SECTIONS["E"] = "Stated in the record but not recomputed on this page"
NOTRERUN = [
 ("J1/G1", "The band test on the massive elementary particles (record J1, G1: a weak hint, p about a tenth)", "no script for this test is among the downloads yet"),
 ("G1", "The golden-ratio ladder test", "no script among the downloads yet"),
 ("M4/M13/M17", "The scorecards (the record's own PASS and PARTIAL counts)", "these are the record's own judgements; the stages above are the evidence they rest on"),
 ("P1", "The comparison count for the Standard Model plus gravity and dark energy", "not recounted here"),
 ("K8", "That the neutron-minus-proton mass and the von Laue residual are built in by construction", "this is a reading of Gemini's code, not a number"),
]
for cid, ttl, why in NOTRERUN:
    check("X-" + cid.replace("/", "-"), "E", ttl, (R + "part-1-every-idea-step-by-step", "record"), ["—"], "not run", [why], "as the record states", "not applicable", "NOT RE-RUN")

# ================================================================================================ SECTION F: the published scripts, run as published
SECTIONS["F"] = "The published scripts, run as published"
for n_, pre_ in enumerate(PUBLISHED):
    name, rc, out_, err_ = PUB[pre_]
    miss = re.findall(r"No module named '(\w+)'", err_)
    ran = rc == 0 and out_.strip() != ""
    tail = err_.strip().splitlines()[-1] if err_.strip() else ""
    check(f"F{n_ + 1}", "F", f"{name}", (R + "chapter-q-the-line-by-line-check-and-the-neutron-star-prediction-14-nov-evening", "record Chapter Q"),
          [f"python -I {name}   (working directory separate from the downloads)"],
          "run with the isolated interpreter in a separate working directory, a copy of the workbook beside it",
          ([f"Exit status {I(rc)}: ran to the end."] if ran else [f"Exit status {I(rc if rc is not None else 0)}: DID NOT RUN HERE. {tail}"]),
          "the script's own printed output follows, unedited", "the script must run to the end in this environment",
          "PASS" if ran else "FAIL",
          ("" if ran else (f"The module {', '.join(miss)} is not installed in this environment, so this script is not verified here. The page computes the same result itself in stages B5 and K10a." if miss else "")),
          pre=(out_.strip() + (("\n--- error output ---\n" + err_.strip()[-1200:]) if err_.strip() else "")))

# ================================================================================================ render
def esc(s): return html.escape(s, quote=False)
def badge(st): return f'<span class="badge badge-{CLS[st]}">{st}</span>'
counts = {s: sum(1 for c in CHECKS if c["status"] == s) for s in STATUSES}
def render_check(c):
    out = [f'<article class="check" id="c-{esc(c["id"])}">',
           f'<header><h3>{esc(c["id"])}. {esc(c["title"])}</h3>{badge(c["status"])}</header>']
    out += [f'<div class="eq">{esc(e)}</div>' for e in c["eqs"]]
    out.append('<dl class="calc">')
    out.append(f'<dt>Numbers in</dt><dd>{esc(c["subst"])}</dd>')
    out.append('<dt>Computed</dt><dd>' + "<br>".join(esc(l) for l in c["computed"]) + '</dd>')
    out.append(f'<dt>Measured</dt><dd>{esc(c["measured"])}</dd>')
    out.append(f'<dt>Tolerance</dt><dd>{esc(c["tol"])}</dd>')
    out.append(f'<dt>Result</dt><dd>{badge(c["status"])}' + (f' {esc(c["note"])}' if c["note"] else "") + '</dd>')
    out.append('</dl>')
    if c.get('pre'): out.append('<pre class="src"><code>' + esc(c['pre']) + '</code></pre>')
    out.append(f'<p class="note">In the record: <a href="{c["rec"][0]}">{esc(c["rec"][1])}</a></p></article>')
    return "\n".join(out)

toc = ["<ul>"]
for key, name in SECTIONS.items():
    toc.append(f'<li><a href="#sec-{key.lower()}">{esc(name)}</a>: ' + ", ".join(f'<a href="#c-{esc(c["id"])}">{esc(c["id"])}</a> {badge(c["status"])}' for c in CHECKS if c["sec"] == key) + "</li>")
toc.append("</ul>")
summary_rows = "".join(f'<tr><td>{badge(s)}</td><td>{I(counts[s])}</td></tr>' for s in STATUSES)
body = [
 "<h1>The math, step by step</h1>",
 '<p class="lede">Every equation on this site, written out, with the numbers put in, the answer, the measured value, and a verdict you can check.</p>',
 f'<p>This page was written by a program, not by hand. The program reads the measured values (each with its source) from one data file, does every calculation itself, and then decides PASS or FAIL by comparing with a tolerance that is printed next to the answer. No number on this page is typed in. Where the record says a result failed or is only partly there, this page shows FAIL or PARTIAL. The same code is on the <a href="code.html">code page</a>, so you can copy it and run it yourself.</p>',
 f'<p>How to read a verdict. <strong>PASS</strong>: the computed value is inside the stated tolerance. <strong>PARTIAL</strong>: outside it but within {N(PF,2)} times the tolerance, or a hint the record itself calls partial. <strong>FAIL</strong>: outside it. <strong>NOT TESTABLE YET</strong>: no measurement or no TRPG number exists to compare. <strong>INPUT</strong>: a measured value used as an input, not a prediction. <strong>INHERITED</strong>: taken unchanged from established physics. <strong>NOT RE-RUN</strong>: stated in the record but not recomputed here.</p>',
 f'<p>The tally is {I(len(CHECKS))} stages in all: ' + ", ".join(f"{I(counts[s])} {s}" for s in STATUSES if counts[s]) + f'. A verdict applies to one stage only. A PASS on a stage does not make a theory true, and PASS counts here are not the record\'s scorecard (see <a href="record.html#{ANCH["M17"]}">record M17</a>).</p>',
 f'<p>Where the record states a number and this run gives a different one, the stage says so. Number of stages where that happened: {I(len(MISMATCH))}.</p>',
 '<table class="summary"><thead><tr><th>Verdict</th><th>Stages</th></tr></thead><tbody>' + summary_rows + '</tbody></table>',
 '<nav class="toc" aria-label="Stages"><h2>Jump to a stage</h2>' + "\n".join(toc) + '</nav>',
]
for key, name in SECTIONS.items():
    body.append(f'<h2 id="sec-{key.lower()}">{esc(name)}</h2>')
    body += [render_check(c) for c in CHECKS if c["sec"] == key]
body_html = "\n".join(body)
available = {p.name for p in BUILD.glob("*.html")} | {"math.html"}
out = sc.page("math.html", "The math, step by step",
              "Every equation on the TRPG site with the numbers substituted, the computed result, the measured value, the tolerance and a PASS or FAIL for each stage, all computed by the published code.",
              body_html, available)

# ------------------------------------------------------------------------------------------------ provenance check
def main_text(page_html):
    m = re.search(r"<main[^>]*>(.*)</main>", page_html, re.S).group(1)
    m = re.sub(r"<[^>]+>", " ", m)
    return html.unescape(m)
def provenance(page_html):
    bad = []
    for tok in NUM.findall(main_text(page_html)):
        if tok not in EMITTED and tok not in FORMULA: bad.append(tok)
    return bad
sb = static_check()
if sb:
    print("STATIC CHECK FAILED:", sb[:10]); sys.exit(4)
bad = provenance(out)
if bad:
    mt = main_text(out)
    print("PROVENANCE CHECK FAILED: numbers not emitted by the generator:", sorted(set(bad)))
    for tok in sorted(set(bad)):
        mm_ = re.search(r"(?<![A-Za-z_\d.])" + re.escape(tok) + r"(?!\d|\.\d)", mt); i = mm_.start() if mm_ else 0; print("  ", tok, "->", ascii(mt[max(0, i - 50):i + 30]))
    sys.exit(2)
if "--plant" in sys.argv:
    planted = out.replace("</main>", "<p>A typed number " + "".join(["7"] * 4) + " slipped in.</p></main>")
    if not provenance(planted):
        print("PLANT NOT CAUGHT: provenance check is blind"); sys.exit(3)
    print("plant caught: the provenance check rejected a typed number as it must:", sorted(set(provenance(planted))))
(BUILD / "math.html").write_text(out, encoding="utf-8")
sc.rewrite_navs(BUILD); sc.update_sitemap(BUILD, ["math.html"])
print(f"math.html written: {len(CHECKS)} stages; " + ", ".join(f"{counts[s]} {s}" for s in STATUSES if counts[s]))
print("record mismatches:", MISMATCH)
print(f"provenance: {len(set(NUM.findall(main_text(out))))} distinct numbers on the page, all emitted by the generator or part of a printed formula")

{
 "_about": "Every constant, measured reference value and tolerance used by build_math_page.py. math.html contains no number that is not either computed from these or part of an equation's own formula. Each entry: value, unit, source (as the full record or its scripts cite it).",
 "hbarc_MeVfm": {
  "value": 197.3269804,
  "unit": "MeV fm",
  "source": "CODATA (as in the project's scripts)"
 },
 "hc_eVm": {
  "value": 1.23984198e-06,
  "unit": "eV m",
  "source": "CODATA (gravity_properties script)"
 },
 "G_SI": {
  "value": 6.674e-11,
  "unit": "m^3 kg^-1 s^-2",
  "source": "CODATA (as in the project's scripts)"
 },
 "c_SI": {
  "value": 299800000.0,
  "unit": "m/s",
  "source": "CODATA, rounded as in the project's scripts"
 },
 "c_exact": {
  "value": 299792458.0,
  "unit": "m/s",
  "source": "CODATA (cosmos_T4 script)"
 },
 "Msun_kg": {
  "value": 1.989e+30,
  "unit": "kg",
  "source": "IAU nominal, as in the project's scripts"
 },
 "Rsun_m": {
  "value": 695700000.0,
  "unit": "m",
  "source": "IAU nominal, as in the project's scripts"
 },
 "arcsec_per_rad": {
  "value": 206265,
  "unit": "arcsec/rad",
  "source": "definition"
 },
 "MeV_J": {
  "value": 1.602176634e-13,
  "unit": "J/MeV",
  "source": "CODATA"
 },
 "Mpc_m": {
  "value": 3.0857e+22,
  "unit": "m",
  "source": "IAU"
 },
 "Mpc_km": {
  "value": 3.0857e+19,
  "unit": "km",
  "source": "IAU"
 },
 "Gyr_s": {
  "value": 3.156e+16,
  "unit": "s",
  "source": "definition (as in the project's script)"
 },
 "AU_m": {
  "value": 149600000000.0,
  "unit": "m",
  "source": "IAU"
 },
 "Mars_orbit_in_AU": {
  "value": 1.52,
  "unit": "AU",
  "source": "as in the resistance_waves script"
 },
 "Msun_geo_km": {
  "value": 1.4766,
  "unit": "km",
  "source": "G Msun / c^2 (as in the project's TOV scripts)"
 },
 "tov_conv": {
  "value": 1.3234e-06,
  "unit": "km^-2 per MeV/fm^3",
  "source": "unit conversion (as in the project's TOV scripts)"
 },
 "Msun_s": {
  "value": 4.925491e-06,
  "unit": "s",
  "source": "G Msun / c^3 (ringdown script)"
 },
 "Msun_m": {
  "value": 1476.625,
  "unit": "m",
  "source": "G Msun / c^2 (ringdown script)"
 },
 "planck_length_m": {
  "value": 1.616255e-35,
  "unit": "m",
  "source": "CODATA"
 },
 "m_e": {
  "value": 0.51099895,
  "unit": "MeV",
  "source": "PDG 2024"
 },
 "m_mu": {
  "value": 105.6583755,
  "unit": "MeV",
  "source": "PDG 2024"
 },
 "m_tau": {
  "value": 1776.93,
  "unit": "MeV",
  "source": "PDG 2024"
 },
 "m_proton": {
  "value": 938.272,
  "unit": "MeV",
  "source": "PDG 2024"
 },
 "m_proton_kg": {
  "value": 1.6726e-27,
  "unit": "kg",
  "source": "CODATA"
 },
 "sigma_m_mu": {
  "value": 2.3e-06,
  "unit": "MeV",
  "source": "PDG 2024 uncertainty (twist_and_link script)"
 },
 "sigma_m_tau": {
  "value": 0.09,
  "unit": "MeV",
  "source": "PDG 2024 uncertainty (twist_and_link script)"
 },
 "q_u": {
  "value": 2.16,
  "unit": "MeV",
  "source": "PDG 2024, MS-bar at 2 GeV"
 },
 "q_d": {
  "value": 4.7,
  "unit": "MeV",
  "source": "PDG 2024, MS-bar at 2 GeV"
 },
 "q_s": {
  "value": 93.5,
  "unit": "MeV",
  "source": "PDG 2024, MS-bar at 2 GeV"
 },
 "q_c": {
  "value": 1273,
  "unit": "MeV",
  "source": "PDG 2024, MS-bar at its own mass"
 },
 "q_b": {
  "value": 4183,
  "unit": "MeV",
  "source": "PDG 2024, MS-bar at its own mass"
 },
 "q_t": {
  "value": 172570,
  "unit": "MeV",
  "source": "PDG 2024, pole mass"
 },
 "ckm_Vus": {
  "value": 0.2243,
  "unit": "",
  "source": "PDG 2024 (as the record cites for the Cabibbo angle)"
 },
 "ckm_Vub": {
  "value": 0.00382,
  "unit": "",
  "source": "PDG 2024 (build_antiparticle_and_decays script)"
 },
 "ckm_Vcb": {
  "value": 0.0408,
  "unit": "",
  "source": "PDG 2024 (build_antiparticle_and_decays script)"
 },
 "ckm_Vtd": {
  "value": 0.0086,
  "unit": "",
  "source": "PDG 2024 (build_antiparticle_and_decays script)"
 },
 "ckm_Vts": {
  "value": 0.0415,
  "unit": "",
  "source": "PDG 2024 (build_antiparticle_and_decays script)"
 },
 "tau_tau_fs": {
  "value": 290.3,
  "unit": "fs",
  "source": "PDG 2024"
 },
 "tau_mu_s": {
  "value": 2.1969811e-06,
  "unit": "s",
  "source": "PDG 2024"
 },
 "BR_tau_e": {
  "value": 0.1782,
  "unit": "",
  "source": "PDG 2024, tau to e nu nu"
 },
 "GF_MeV2": {
  "value": 1.1663788e-11,
  "unit": "MeV^-2",
  "source": "PDG 2024 (Fermi constant)"
 },
 "hbar_MeVs": {
  "value": 6.582119569e-22,
  "unit": "MeV s",
  "source": "CODATA"
 },
 "H0": {
  "value": 67.4,
  "unit": "km/s/Mpc",
  "source": "Planck 2018"
 },
 "Omega_m": {
  "value": 0.315,
  "unit": "",
  "source": "Planck 2018"
 },
 "Omega_r": {
  "value": 9.1e-05,
  "unit": "",
  "source": "Planck 2018 (as in the project's script)"
 },
 "Omega_L_stated": {
  "value": 0.685,
  "unit": "",
  "source": "Planck 2018"
 },
 "Omega_k": {
  "value": 0.001,
  "unit": "",
  "source": "Planck 2018"
 },
 "Omega_k_sigma": {
  "value": 0.002,
  "unit": "",
  "source": "Planck 2018"
 },
 "age_measured": {
  "value": 13.787,
  "unit": "Gyr",
  "source": "Planck 2018, as record Q2 cites it (13.787 +/- 0.020)"
 },
 "age_sigma": {
  "value": 0.02,
  "unit": "Gyr",
  "source": "Planck 2018, as record Q2 cites it"
 },
 "bending_measured": {
  "value": 1.75,
  "unit": "arcsec",
  "source": "as the record cites (Einstein 1915; eclipse 1919; Cassini 2003)"
 },
 "gamma_minus_1": {
  "value": 2.1e-05,
  "unit": "",
  "source": "Cassini (Bertotti, Iess and Tortora 2003), as record M9 cites"
 },
 "gamma_minus_1_sigma": {
  "value": 2.3e-05,
  "unit": "",
  "source": "Cassini (Bertotti, Iess and Tortora 2003), as record M9 cites"
 },
 "ns_mass_measured": {
  "value": 2.08,
  "unit": "Msun",
  "source": "PSR J0740+6620, as record K1 cites"
 },
 "ns_mass_sigma": {
  "value": 0.07,
  "unit": "Msun",
  "source": "PSR J0740+6620, as record K1 cites"
 },
 "ns_radius_lo": {
  "value": 11,
  "unit": "km",
  "source": "NICER and GW170817, as record K1 cites"
 },
 "ns_radius_hi": {
  "value": 13,
  "unit": "km",
  "source": "NICER and GW170817, as record K1 cites"
 },
 "lambda_limit": {
  "value": 580,
  "unit": "",
  "source": "GW170817 (2017), as record K2 cites"
 },
 "sly_published_mmax": {
  "value": 2.05,
  "unit": "Msun",
  "source": "SLy equation of state (Read, Lackey, Owen and Friedman 2009), as record K1 cites"
 },
 "sly_published_r14": {
  "value": 11.7,
  "unit": "km",
  "source": "SLy equation of state (Read, Lackey, Owen and Friedman 2009), as record K1 cites"
 },
 "sly_published_lambda": {
  "value": 300,
  "unit": "",
  "source": "SLy, published (Hinderer 2008 method), as record K2 cites"
 },
 "bag_B_gemini": {
  "value": 54.1,
  "unit": "MeV/fm^3",
  "source": "Gemini's V53 value, as record K1 cites"
 },
 "bag_wall_gemini": {
  "value": 0.00042,
  "unit": "",
  "source": "Gemini's wall strength, as record K1 and K2 cite"
 },
 "bag_B_fj_min": {
  "value": 57,
  "unit": "MeV/fm^3",
  "source": "Farhi and Jaffe (1984) stability bound, as record K1/K2 cite (fourth root about 145 MeV)"
 },
 "bag_B_window_lo": {
  "value": 62,
  "unit": "MeV/fm^3",
  "source": "record K2 neutron-star window"
 },
 "bag_B_window_hi": {
  "value": 70,
  "unit": "MeV/fm^3",
  "source": "record K2 neutron-star window"
 },
 "bag_B_mid": {
  "value": 66,
  "unit": "MeV/fm^3",
  "source": "middle of the K2 window, as in the master_action_M1 script"
 },
 "bag_walls_tidal": {
  "value": [
   0.0,
   0.0001,
   0.0002,
   0.0003
  ],
  "unit": "",
  "source": "wall strengths scanned in the record's K2 text and the tidal script's own scan"
 },
 "gemini_p_scale": {
  "value": 0.12,
  "unit": "",
  "source": "Gemini V53 low-density formula, as the TOV script reads it"
 },
 "gemini_p_power": {
  "value": 0.65,
  "unit": "",
  "source": "Gemini V53 low-density formula, as the TOV script reads it"
 },
 "sly_G": {
  "value": [
   3.005,
   2.988,
   2.851
  ],
  "unit": "",
  "source": "Read, Lackey, Owen and Friedman (2009), SLy piecewise polytrope"
 },
 "sly_log_rho12": {
  "value": [
   14.7,
   15.0
  ],
  "unit": "log10 g/cm^3",
  "source": "Read, Lackey, Owen and Friedman (2009)"
 },
 "sly_log_p1": {
  "value": 34.384,
  "unit": "log10 dyn/cm^2",
  "source": "Read, Lackey, Owen and Friedman (2009)"
 },
 "sly_crust": {
  "value": [
   [
    6.8011e-09,
    1.58425,
    24403400.0
   ],
   [
    1.06186e-06,
    1.28733,
    378358000000.0
   ],
   [
    53.2697,
    0.62223,
    2627800000000.0
   ],
   [
    3.99874e-08,
    1.35692,
    null
   ]
  ],
  "unit": "cgs",
  "source": "Read, Lackey, Owen and Friedman (2009), SLy crust pieces"
 },
 "c2_cgs": {
  "value": 8.98755e+20,
  "unit": "cm^2/s^2",
  "source": "c squared in cgs"
 },
 "mevfm3_gcm3": {
  "value": 1782700000000.0,
  "unit": "g/cm^3 per MeV/fm^3",
  "source": "unit conversion"
 },
 "glueball_lo": {
  "value": 1710,
  "unit": "MeV",
  "source": "Morningstar and Peardon (1999)"
 },
 "glueball_hi": {
  "value": 1730,
  "unit": "MeV",
  "source": "Chen and others (2006)"
 },
 "svz_G2": {
  "value": 0.012,
  "unit": "GeV^4",
  "source": "Shifman, Vainshtein and Zakharov (1979) gluon condensate, as in the master_action_M1 script"
 },
 "svz_factor": {
  "value": 2.0,
  "unit": "",
  "source": "tolerance stated on the page: factor of agreement"
 },
 "proton_radius_fm": {
  "value": 0.841,
  "unit": "fm",
  "source": "measured charge radius, as record K3 cites"
 },
 "core_radius_lo": {
  "value": 0.68,
  "unit": "fm",
  "source": "record K3 (rms core for B in the window)"
 },
 "core_radius_hi": {
  "value": 0.71,
  "unit": "fm",
  "source": "record K3 (rms core for B in the window)"
 },
 "nucleon_mass_avg": {
  "value": 938.9,
  "unit": "MeV",
  "source": "PDG, nucleon mass as record K3 uses it"
 },
 "delta_mass": {
  "value": 1232,
  "unit": "MeV",
  "source": "PDG, Delta mass as record K3 uses it"
 },
 "np_delta_gap": {
  "value": 293.46,
  "unit": "MeV",
  "source": "measured Delta minus nucleon gap, as record K3 cites"
 },
 "gemini_C_ss": {
  "value": 48.91,
  "unit": "MeV",
  "source": "Gemini V50 typed constant, as record K3/K8 cite"
 },
 "gemini_f": {
  "value": 2.18,
  "unit": "",
  "source": "Gemini V50 typed coupling, as record K8 cites"
 },
 "gemini_width_claim": {
  "value": 114.2,
  "unit": "MeV",
  "source": "Gemini V50 claim, as record K8 cites"
 },
 "delta_width_measured": {
  "value": 117,
  "unit": "MeV",
  "source": "PDG, as record K8 cites"
 },
 "pion_mass": {
  "value": 139.57,
  "unit": "MeV",
  "source": "PDG"
 },
 "gemini_gA_claim": {
  "value": 1.2781,
  "unit": "",
  "source": "Gemini V51 claim, as record K8 cites"
 },
 "gA_measured": {
  "value": 1.2754,
  "unit": "",
  "source": "PDG, as record K8 cites"
 },
 "gemini_Z2_denominator": {
  "value": 1.214,
  "unit": "",
  "source": "Gemini V51 dressing constant, as the proton_claims script reads it"
 },
 "gemini_gA_bare_sq": {
  "value": 1.5129,
  "unit": "",
  "source": "Gemini V51 bare coupling squared, as the proton_claims script reads it"
 },
 "gemini_gA_delta": {
  "value": 0.214,
  "unit": "",
  "source": "Gemini V51 constant, as the proton_claims script reads it"
 },
 "gemini_gA_hand_add": {
  "value": 0.0517,
  "unit": "",
  "source": "Gemini V51 hand-added term, as record K8 cites"
 },
 "tolerance_gA": {
  "value": 0.01,
  "unit": "relative",
  "source": "tolerance stated on the page"
 },
 "jlab_pressure_switch_fm": {
  "value": 0.6,
  "unit": "fm",
  "source": "Burkert, Elouadrhiri and Girod (2018), as record K8 cites"
 },
 "dterm_bag": {
  "value": -1.1,
  "unit": "",
  "source": "Neubelt, Schweitzer and others (2020), as record K8 cites"
 },
 "dterm_jlab": {
  "value": -3.0,
  "unit": "",
  "source": "Jefferson Lab estimate, as record K8 cites"
 },
 "ringdown_M_source": {
  "value": 62,
  "unit": "Msun",
  "source": "LIGO/Virgo (2016), GW150914 remnant"
 },
 "ringdown_z": {
  "value": 0.09,
  "unit": "",
  "source": "LIGO/Virgo (2016)"
 },
 "ringdown_spin": {
  "value": 0.67,
  "unit": "",
  "source": "LIGO/Virgo (2016)"
 },
 "kerr_fit": {
  "value": [
   1.5251,
   1.1568,
   0.1292,
   0.7,
   1.4187,
   0.499
  ],
  "unit": "",
  "source": "Berti, Cardoso and Will (2006) fit for the l=m=2 mode"
 },
 "ringdown_f_measured": {
  "value": 251,
  "unit": "Hz",
  "source": "LIGO measurement of GW150914, as record K9 cites"
 },
 "ringdown_f_sigma": {
  "value": 8,
  "unit": "Hz",
  "source": "LIGO measurement of GW150914, as record K9 cites"
 },
 "ringdown_tau_measured": {
  "value": 4.0,
  "unit": "ms",
  "source": "LIGO measurement of GW150914, as record K9 cites"
 },
 "ringdown_tau_sigma": {
  "value": 0.3,
  "unit": "ms",
  "source": "LIGO measurement of GW150914, as record K9 cites"
 },
 "compactness_sun": {
  "value": 4.24e-06,
  "unit": "",
  "source": "ringdown_check script"
 },
 "sirius_b_mass": {
  "value": 1.02,
  "unit": "Msun",
  "source": "ringdown_check script"
 },
 "sirius_b_radius_m": {
  "value": 5800000.0,
  "unit": "m",
  "source": "ringdown_check script"
 },
 "buchdahl_limit_num": {
  "value": 8,
  "unit": "",
  "source": "Buchdahl (1959)"
 },
 "buchdahl_limit_den": {
  "value": 9,
  "unit": "",
  "source": "Buchdahl (1959)"
 },
 "dm21": {
  "value": 7.49e-05,
  "unit": "eV^2",
  "source": "NuFIT 6.0 (2024), normal ordering"
 },
 "dm31": {
  "value": 0.002513,
  "unit": "eV^2",
  "source": "NuFIT 6.0 (2024), normal ordering"
 },
 "dm32_inverted": {
  "value": -0.002484,
  "unit": "eV^2",
  "source": "NuFIT 6.0 (2024), inverted ordering"
 },
 "nu_sum_limit": {
  "value": 0.064,
  "unit": "eV",
  "source": "DESI collaboration (2025) cosmology limit, as record M2 cites"
 },
 "nu_masses_stated": {
  "value": [
   0.36,
   8.7,
   50.1
  ],
  "unit": "meV",
  "source": "the record's stated prediction (home page, M2)"
 },
 "tol_lepton_mass": {
  "value": 0.0001,
  "unit": "relative",
  "source": "tolerance stated on the page"
 },
 "tol_percent_1": {
  "value": 0.01,
  "unit": "relative",
  "source": "tolerance stated on the page"
 },
 "tol_percent_5": {
  "value": 0.05,
  "unit": "relative",
  "source": "tolerance stated on the page"
 },
 "tol_percent_10": {
  "value": 0.1,
  "unit": "relative",
  "source": "tolerance stated on the page"
 },
 "tol_identity": {
  "value": 1e-09,
  "unit": "relative",
  "source": "tolerance stated on the page (exact identities)"
 },
 "partial_factor": {
  "value": 3,
  "unit": "",
  "source": "page rule: PARTIAL means within this many times the stated tolerance"
 },
 "koide_target_num": {
  "value": 2,
  "unit": "",
  "source": "Koide (1981)"
 },
 "koide_target_den": {
  "value": 3,
  "unit": "",
  "source": "Koide (1981)"
 },
 "brannen_delta_num": {
  "value": 2,
  "unit": "",
  "source": "Brannen (2006)"
 },
 "brannen_delta_den": {
  "value": 9,
  "unit": "",
  "source": "Brannen (2006)"
 },
 "random_triples": {
  "value": 2000000,
  "unit": "",
  "source": "sample size of the project's random test (mass_formula_fix script)"
 },
 "random_seed": {
  "value": 1,
  "unit": "",
  "source": "seed of the project's random test"
 },
 "random_range_log10": {
  "value": [
   -1,
   4
  ],
  "unit": "log10 MeV",
  "source": "range of the project's random test"
 },
 "random_chance_stated": {
  "value": 3.5e-06,
  "unit": "",
  "source": "as record K7 states"
 },
 "gemini_knob_count_items": {
  "value": [
   "E0 (fitted to the proton)",
   "P(0)",
   "the /10",
   "anchor 17",
   "anchor 40",
   "v = 174.1 GeV",
   "coefficient 0.218",
   "coefficient 0.415",
   "coefficient 0.439",
   "coefficient 1.71",
   "coefficient 2.66",
   "coefficient 0.9925"
  ],
  "unit": "",
  "source": "as record K7 lists them"
 },
 "gemini_masses_count": {
  "value": 9,
  "unit": "",
  "source": "as record K7 states (9 masses)"
 },
 "proton_scale_hint": {
  "value": 0.0035,
  "unit": "relative",
  "source": "tolerance stated on the page (record K7: within 0.35% of the proton mass)"
 },
 "R_universe_hubble_factor": {
  "value": 1.0,
  "unit": "",
  "source": "Hubble radius = c / H0"
 },
 "eddington_stated": {
  "value": 1e+80,
  "unit": "",
  "source": "Eddington (1930s), as record K6 cites"
 },
 "weinberg_relation_exponent_den": {
  "value": 3,
  "unit": "",
  "source": "Weinberg (1972)"
 },
 "gluon_vev_eps_prefactor_num": {
  "value": 9,
  "unit": "",
  "source": "SVZ vacuum energy formula as in the master_action_M1 script"
 },
 "gluon_vev_eps_prefactor_den": {
  "value": 32,
  "unit": "",
  "source": "SVZ vacuum energy formula as in the master_action_M1 script"
 },
 "de_scale_meV_stated": {
  "value": 2.24,
  "unit": "meV",
  "source": "as record M8 states"
 },
 "de_density_stated": {
  "value": 3.3e-42,
  "unit": "MeV/fm^3",
  "source": "as record M8 states"
 },
 "wave_n_grid": {
  "value": 1201,
  "unit": "",
  "source": "grid size in the wave_grid_geometry script"
 },
 "gauss_bonnet_tol": {
  "value": 1e-06,
  "unit": "",
  "source": "tolerance stated on the page"
 },
 "resistance_layers": {
  "value": [
   10,
   100,
   1000,
   3000
  ],
  "unit": "",
  "source": "layer counts in the resistance_waves script"
 },
 "localization_factor": {
  "value": 100,
  "unit": "",
  "source": "tolerance stated on the page: transmission must fall by at least this factor"
 },
 "bands_min_gaps": {
  "value": 4,
  "unit": "",
  "source": "tolerance stated on the page (record M7 reports four gaps)"
 },
 "shapiro_roundtrip_range_us": {
  "value": [
   100,
   250
  ],
  "unit": "microseconds",
  "source": "range quoted in the resistance_waves script for the radar test (Shapiro 1964, Viking 1979)"
 },
 "scalar_resistance_factor": {
  "value": 0.0,
  "unit": "",
  "source": "gamma for scalar resistance (Nordstrom 1913), as record M9 treats it"
 },
 "p_ripple_eV": {
  "value": 1720000000.0,
  "unit": "eV",
  "source": "lightest glueball (lattice QCD), as record M10 cites"
 },
 "graviton_limit_eV": {
  "value": 1.27e-23,
  "unit": "eV",
  "source": "LIGO-Virgo-KAGRA GWTC-3 (2021), as record M10 cites"
 },
 "gravity_reach_needed_m": {
  "value": 1e+16,
  "unit": "m",
  "source": "tolerance stated on the page: gravity must reach at least this far"
 },
 "planck_curvature_limit": {
  "value": 0.002,
  "unit": "",
  "source": "Planck 2018"
 },
 "entangle_density_stated": {
  "value": 1e+68,
  "unit": "per m^2",
  "source": "as record M14 states (about 10^68 per square metre)"
 },
 "factor_three_tol": {
  "value": 3,
  "unit": "",
  "source": "tolerance stated on the page: within this factor"
 },
 "grid_N": {
  "value": 36,
  "unit": "",
  "source": "lattice size in the area_law script"
 },
 "grid_mass": {
  "value": 0.01,
  "unit": "",
  "source": "oscillator mass in the area_law script"
 },
 "area_blocks": {
  "value": [
   2,
   4,
   6,
   8,
   10,
   12,
   14
  ],
  "unit": "",
  "source": "block sizes in the area_law script"
 },
 "area_blocks_tidy": {
  "value": [
   4,
   6,
   8,
   10,
   12,
   14
  ],
  "unit": "",
  "source": "block sizes in the gravity_tighten script"
 },
 "crossing_amplitude": {
  "value": 0.5,
  "unit": "",
  "source": "crossing-wave coupling in the gravity_tighten script"
 },
 "layer_masses": {
  "value": [
   0.01,
   0.1,
   0.3,
   1.0
  ],
  "unit": "",
  "source": "layer masses in the gravity_tighten script"
 },
 "area_ratio_min": {
  "value": 10,
  "unit": "",
  "source": "tolerance stated on the page: boundary fit must beat the volume fit by at least this factor"
 },
 "area_ratio_min_crossing": {
  "value": 3,
  "unit": "",
  "source": "tolerance stated on the page"
 },
 "small_world_cases": {
  "value": [
   [
    8000000000.0,
    150
   ],
   [
    1e+80,
    10
   ],
   [
    1e+80,
    1000
   ]
  ],
  "unit": "",
  "source": "N and links per node in the entanglement_network script (people, particles, particles)"
 },
 "small_world_stated": {
  "value": [
   4.6,
   80,
   27
  ],
  "unit": "steps",
  "source": "as record M15 states"
 },
 "network_stated": {
  "value": [
   3,
   1.33
  ],
  "unit": "steps",
  "source": "as record M15 states (longest route, average)"
 },
 "xlsx_stated_counts": {
  "value": [
   19,
   171,
   55,
   21,
   12
  ],
  "unit": "",
  "source": "as record K5 states (particles, pairs, electric, strong, weak)"
 },
 "xlsx_fate_stated": {
  "value": [
   19,
   43,
   24,
   294
  ],
  "unit": "",
  "source": "as record K10 states (direct opposites, annihilate, fuse to W, scatter only)"
 },
 "anomaly_Y": {
  "value": {
   "YQ": [
    1,
    6
   ],
   "Yu": [
    2,
    3
   ],
   "Yd": [
    -1,
    3
   ],
   "YL": [
    -1,
    2
   ],
   "Ye": [
    -1,
    1
   ]
  },
  "unit": "",
  "source": "Standard Model hypercharges (Geng and Marshak 1989; Minahan, Ramond and Warner 1990)"
 },
 "colours_tested": {
  "value": [
   1,
   2,
   3,
   4,
   5
  ],
  "unit": "",
  "source": "record M16 colour counts"
 },
 "scale_factors": {
  "value": [
   0.37,
   2.9,
   41.0
  ],
  "unit": "",
  "source": "rescaling factors for the unit-independence test of S (any positive numbers work)"
 },
 "quark_angle_stated": {
  "value": {
   "lepton": 45.0,
   "down": 47.5,
   "up": 51.2
  },
  "unit": "degrees",
  "source": "as record M18 states"
 },
 "mu_quark_run": {
  "value": {
   "alpha_s_MZ": 0.118,
   "MZ": 91.1876,
   "mc": 1.273,
   "mb": 4.183,
   "mt": 162.5
  },
  "unit": "GeV",
  "source": "PDG 2024 (quark_common_scale script)"
 },
 "quark_run_inputs": {
  "value": {
   "u": [
    0.00216,
    2.0
   ],
   "d": [
    0.0047,
    2.0
   ],
   "s": [
    0.0935,
    2.0
   ],
   "c": [
    1.273,
    1.273
   ],
   "b": [
    4.183,
    4.183
   ],
   "t": [
    162.5,
    162.5
   ]
  },
  "unit": "GeV (mass, scale)",
  "source": "PDG 2024 (quark_common_scale script)"
 },
 "quark_run_scales": {
  "value": [
   2.0,
   91.1876,
   162.5,
   1000.0
  ],
  "unit": "GeV",
  "source": "energies in the quark_common_scale script"
 },
 "quark_angle_run_stated": {
  "value": {
   "up": 52.15,
   "down": 48.06
  },
  "unit": "degrees",
  "source": "as record P2 states"
 },
 "waterfall_stated": {
  "value": {
   "Q_cbt": 0.669,
   "charm_pred": 1353,
   "Q_scb": 0.675,
   "strange_pred": 91.8
  },
  "unit": "",
  "source": "as record M6 states"
 },
 "inputs_count_items": {
  "value": [
   "G (or a)",
   "rho_Lambda",
   "B",
   "glueball mass",
   "lepton scale M0",
   "neutrino splitting 1",
   "neutrino splitting 2",
   "force strength 1",
   "force strength 2",
   "force strength 3",
   "Higgs mass",
   "v_H",
   "quark mass u",
   "quark mass d",
   "quark mass s",
   "quark mass c",
   "quark mass b",
   "quark mass t",
   "quark mixing number 1",
   "quark mixing number 2",
   "quark mixing number 3",
   "quark mixing number 4"
  ],
  "unit": "",
  "source": "record P1 / equation page list of measured inputs"
 },
 "inputs_count_stated": {
  "value": 21,
  "unit": "",
  "source": "as record P1 states (about 21)"
 },
 "twist_samples": {
  "value": 20000,
  "unit": "",
  "source": "sample size in the twist_and_link script"
 },
 "twist_seed": {
  "value": 0,
  "unit": "",
  "source": "seed in the twist_and_link script"
 },
 "twist_sigma_limit": {
  "value": 2,
  "unit": "",
  "source": "tolerance stated on the page (standard deviations)"
 },
 "twist_stated": {
  "value": [
   0.222239,
   4e-05
  ],
  "unit": "rad",
  "source": "as record P5 states"
 },
 "gemini_centroid": {
  "value": 1085,
  "unit": "MeV",
  "source": "Gemini V50 typed centroid, as record K3 cites"
 },
 "gemini_mu_bare": {
  "value": [
   2.5424,
   -1.6949
  ],
  "unit": "nuclear magnetons",
  "source": "Gemini V51 bare values, as the proton_claims script reads them"
 },
 "gemini_mu_typed": {
  "value": [
   2.368,
   -1.518
  ],
  "unit": "nuclear magnetons",
  "source": "Gemini V51 typed values, as the proton_claims script reads them"
 },
 "gemini_mu_claim": {
  "value": [
   2.78,
   -1.93
  ],
  "unit": "nuclear magnetons",
  "source": "Gemini V51 claim, as record K8 cites"
 },
 "gemini_ns_ec_list": {
  "value": [
   280,
   350,
   430,
   520,
   650,
   780,
   920,
   1080,
   1240,
   1380,
   1500,
   1650,
   2000,
   2500,
   3000
  ],
  "unit": "MeV/fm^3",
  "source": "central energy densities in the neutron_star_TOV_fix script"
 },
 "record_claims": {
  "value": {
   "A1": 1.752,
   "A2": 13.79,
   "A3a": 2.07,
   "A3b": 11.4,
   "A4a": 1776.97,
   "A5": 0.2242,
   "A6": 289.7,
   "A8f": 250,
   "A9": 0.0592,
   "B1": 3.5,
   "B3lo": 1160,
   "B3hi": 1210,
   "K10mu": 2.187,
   "K1aR": 47,
   "K1cMw": 2.45,
   "K2a": 299,
   "K2b": 487,
   "K2c": 639,
   "K3B": 37,
   "K6N": 5.5e+79,
   "M7a": 124,
   "P4": 5,
   "B3v": 35,
   "Q3B": 60.1,
   "Q3Blo": 57.1,
   "Q3Bhi": 63.7,
   "Q3core": 0.725,
   "Q3pub": 57.5,
   "Q3check1975": 58.4,
   "Q4mlo": 2.0,
   "Q4mhi": 2.18,
   "Q4Llo": 440,
   "Q4Lhi": 570,
   "Q4Rlo": 10.8,
   "Q4Rhi": 11.5,
   "Q4strong_lo": 1.54,
   "Q4strong_hi": 1.8,
   "Q5up": 0.828,
   "Q5down": 0.619,
   "numin": 0.0588,
   "inputs_trpg": 25,
   "inputs_sm": 28,
   "K3rms_lo": 0.68,
   "K3rms_hi": 0.71
  },
  "unit": "",
  "source": "numbers the record itself states, shown beside this run's numbers: A1: record M7 and the results page; A2: record M11 and the home page; A3a: record K1 and the home page; A3b: record K1 and the home page; A4a: record K7; A5: record K4 and the home page; A6: record K10 and the home page; A8f: record K9 and the home page; A9: record M2 and the home page; B1: record M16; B3lo: record M1; B3hi: record M1; K10mu: record K10; K1aR: record K1; K1cMw: record K1; K2a: record K2; K2b: record K2; K2c: record K2; K3B: record K3; K6N: record K6; M7a: record M7; P4: record P4; Q2/Q3/Q4/Q5: record Chapter Q; B3: record M1/P1/equation page (lambda 1160 to 1210, v about 35 MeV)"
 },
 "tov_scan_max": {
  "value": 3500,
  "unit": "MeV/fm^3",
  "source": "upper end of the central-density scan in the neutron_star_TOV_fix script"
 },
 "nuclear_density_probe": {
  "value": 150.0,
  "unit": "MeV/fm^3",
  "source": "energy density at which the TOV script reports Gemini's pressure"
 },
 "sly_scan": {
  "value": [
   200,
   3000
  ],
  "unit": "MeV/fm^3",
  "source": "central-density range of the SLy scan in the neutron_star_TOV_fix script"
 },
 "expansion_sizes": {
  "value": [
   1e-09,
   0.0001,
   0.5,
   1.0,
   2.0
  ],
  "unit": "fraction of today's size",
  "source": "sizes listed in the expansion_and_entanglement script"
 },
 "gauss_pulse": {
  "value": {
   "width": 2.0,
   "points": 65536,
   "extent": 50
  },
  "unit": "",
  "source": "pulse parameters in the expansion_and_entanglement script"
 },
 "citations": {
  "value": [
   "Newton, Soldner 1801",
   "Nordström 1913",
   "Einstein 1915",
   "Heisenberg 1927",
   "Dirac 1928",
   "Bloch 1928",
   "Kronig and Penney 1931",
   "Dirac 1937",
   "Anderson 1958",
   "Gatto, Sartori and Tonin 1968",
   "Bombelli, Koul, Lee and Sorkin 1986",
   "Srednicki 1993",
   "Watts and Strogatz 1998",
   "Fardon, Nelson and Weiner 2004",
   "MICROSCOPE team 2022",
   "Buchdahl 1959",
   "Planck collaboration 2018",
   "Jefferson Lab 2018",
   "Berti, Cardoso and Will 2006",
   "Morningstar and Peardon 1999",
   "Chen and others 2006",
   "Koide 1981",
   "Foot 1994",
   "Brannen 2006",
   "Read, Lackey, Owen and Friedman 2009",
   "LIGO-Virgo-KAGRA 2021",
   "Eddington 1930s",
   "Weinberg 1972",
   "JUNO running since 2025",
   "Geng and Marshak 1989",
   "Minahan, Ramond and Warner 1990",
   "Hinderer 2008",
   "Farhi and Jaffe 1984",
   "Gemini V50 typed centroid 1085 MeV",
   "DESI 2025",
   "Cassini 2003"
  ],
  "unit": "",
  "source": "citation years exactly as the record and its Credits list give them; the page prose cites these, so the years are loaded data"
 },
 "std_ns_mass": {
  "value": 1.4,
  "unit": "Msun",
  "source": "the standard neutron-star mass used in the record (K1, K2) and in the TOV scripts"
 },
 "horizon_compactness": {
  "value": 1.0,
  "unit": "",
  "source": "black-hole horizon, 2GM/(Rc^2) = 1"
 },
 "percent": {
  "value": 100,
  "unit": "",
  "source": "conversion to percent"
 },
 "tol_sum_one": {
  "value": 1e-12,
  "unit": "",
  "source": "tolerance stated on the page (shares sum)"
 },
 "nuclear_pressure_range": {
  "value": [
   0,
   3
  ],
  "unit": "MeV/fm^3",
  "source": "pressure of real nuclear matter at nuclear density, as record K1 states (about 0 to 3, 2 near saturation)"
 },
 "gemini_ns_ec_list_low": {
  "value": [
   20,
   40,
   80,
   120,
   160,
   200,
   240
  ],
  "unit": "MeV/fm^3",
  "source": "low central energy densities in the neutron_star_TOV_fix script"
 },
 "bag_B_proton": {
  "value": [
   57.1,
   60.1,
   63.7
  ],
  "unit": "MeV/fm^3",
  "source": "proton_sets_B script output as quoted in record Q3 (central value and +/-5% pion-nucleon coupling band)"
 },
 "ns_mass_test": {
  "value": 2.0,
  "unit": "Msun",
  "source": "PSR J0740+6620 and J0348+0432: heaviest star must be at least this, as the neutron_stars_no_knobs script states"
 },
 "pairing_gap_range": {
  "value": [
   70.0,
   90.0
  ],
  "unit": "MeV",
  "source": "the TRPG prediction, record Q4"
 },
 "pairing_theory_range": {
  "value": [
   50,
   150
  ],
  "unit": "MeV",
  "source": "theory estimates of the pairing gap, as record Q4 cites"
 },
 "no_knobs_scale": {
  "value": 3.0,
  "unit": "",
  "source": "renormalisation-scale factor X used for the pairing run in the neutron_stars_no_knobs script"
 },
 "no_knobs_scale_band": {
  "value": [
   2.0,
   3.0,
   4.0
  ],
  "unit": "",
  "source": "scale factors X = 2 to 4, the band in the neutron_stars_no_knobs script"
 },
 "no_knobs_ms": {
  "value": 100.0,
  "unit": "MeV",
  "source": "strange quark mass used in the neutron_stars_no_knobs script"
 },
 "no_knobs_alpha_s_MZ": {
  "value": 0.118,
  "unit": "",
  "source": "PDG 2024, alpha_s(M_Z), as in the neutron_stars_no_knobs script"
 },
 "no_knobs_stars": {
  "value": 45,
  "unit": "",
  "source": "stars per sequence on this page (the script uses seventy; interpolation is in log Lambda either way)"
 },
 "inputs_trpg_items": {
  "value": [
   "G (or a)",
   "rho_Lambda",
   "lepton scale M0",
   "neutrino splitting 1",
   "neutrino splitting 2",
   "force strength 1",
   "force strength 2",
   "force strength 3",
   "Higgs mass",
   "v_H",
   "quark mass u",
   "quark mass d",
   "quark mass s",
   "quark mass c",
   "quark mass b",
   "quark mass t",
   "quark mixing number 1",
   "quark mixing number 2",
   "quark mixing number 3",
   "quark mixing number 4",
   "neutrino mixing angle 1",
   "neutrino mixing angle 2",
   "neutrino mixing angle 3",
   "neutrino phase",
   "strong-force CP angle"
  ],
  "unit": "",
  "source": "record P1 / equation page list, with B and the glueball mass counted as consequences of the strong force (Q3), as record Q2 counts"
 },
 "inputs_sm_items": {
  "value": [
   "quark masses",
   "quark masses",
   "quark masses",
   "quark masses",
   "quark masses",
   "quark masses",
   "charged lepton masses",
   "charged lepton masses",
   "charged lepton masses",
   "CKM",
   "CKM",
   "CKM",
   "CKM",
   "gauge couplings",
   "gauge couplings",
   "gauge couplings",
   "Higgs mass",
   "Higgs vev",
   "theta_QCD",
   "neutrino masses",
   "neutrino masses",
   "neutrino masses",
   "PMNS",
   "PMNS",
   "PMNS",
   "PMNS",
   "G",
   "rho_Lambda"
  ],
  "unit": "",
  "source": "the Standard Model with neutrino masses, gravity and dark energy, as record Q2 states (about 28)"
 },
 "fate_stated": {
  "value": {
   "annihilate": 25,
   "pions": 2,
   "W": 24,
   "flavour": 12,
   "scatter": 298,
   "opposites": 19
  },
  "unit": "",
  "source": "as record Q2 and the data page state (corrected antiparticle rules)"
 },
 "nu_light_stated_meV": {
  "value": 0.36,
  "unit": "meV",
  "source": "the distinctive TRPG number: lightest neutrino mass, as the home page and record Q2 state"
 },
 "script_rundir": {
  "value": "run",
  "unit": "",
  "source": "separate working directory for running the published scripts"
 }
}

"""Shared helpers for the trpgphysics.org build scripts (loaded by path, because `python -I` does not
put the script folder on sys.path).  Header/footer text is copied exactly from the existing pages."""
import html, pathlib, re

NAV_ITEMS = [
    ("index.html", "Home"), ("theory.html", "The theory"), ("equation.html", "Equation"),
    ("results.html", "Results"), ("paper.html", "Paper"), ("record.html", "Record"),
    ("data.html", "Data & code"), ("credits.html", "Credits"), ("publications.html", "Publications"),
    ("math.html", "The math"), ("code.html", "The code"), ("originals.html", "Who did the math first"),
    ("about.html", "About"),
]
COPY = "© 2026 Matthew Gelfman and Melissa Gelfman"
FOOTER = ('<footer class="site">\n  <div class="wrap">\n    <p>&copy; 2026 Matthew Gelfman and Melissa Gelfman. All rights reserved. '
          'Code and data are released under the MIT License. See <a href="copyright.html">Copyright &amp; license</a>.</p>\n'
          '    <p>Dynamic Relational Transitioning Geometry (TRPG) is a perspective built on the work of the physicists listed on the '
          '<a href="credits.html">Credits</a> page. First published on this site 14 November 2026.</p>\n  </div>\n</footer>\n')


def nav_html(current, available):
    parts = []
    for href, label in NAV_ITEMS:
        if href not in available:
            continue
        cur = ' aria-current="page"' if href == current else ""
        parts.append(f'<a href="{href}"{cur}>{label}</a>')
    return '<nav id="nav" aria-label="Main">' + "".join(parts) + "</nav>"


def page(fname, title, desc, body, available):
    d = html.escape(desc, quote=True)
    t = html.escape(title)
    return f'''<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>{t} | TRPG Physics</title>
<meta name="description" content="{d}">
<meta name="author" content="Matthew Gelfman and Melissa Gelfman">
<meta name="copyright" content="{COPY}">
<meta property="og:title" content="{t} | TRPG Physics">
<meta property="og:description" content="{d}">
<meta property="og:type" content="website">
<meta property="og:url" content="https://trpgphysics.org/{fname}">
<link rel="canonical" href="https://trpgphysics.org/{fname}">
<link rel="icon" href="assets/favicon.svg" type="image/svg+xml">
<link rel="stylesheet" href="assets/style.css?v=20261114a">
</head>
<body>
<a class="skip" href="#main">Skip to content</a>
<header class="site">
  <div class="wrap bar">
    <a class="brand" href="index.html"><img src="assets/favicon.svg" alt="" width="28" height="28"><span>TRPG Physics</span></a>
    {nav_html(fname, available)}
  </div>
</header>
<main id="main" class="wrap">

{body}
</main>
{FOOTER}</body>
</html>
'''


def rewrite_navs(build):
    """Replace the <nav id="nav">...</nav> on every page with the current nav (only links whose page exists)."""
    build = pathlib.Path(build)
    available = {p.name for p in build.glob("*.html")}
    for p in sorted(build.glob("*.html")):
        t = p.read_text(encoding="utf-8")
        new = re.sub(r'<nav id="nav"[^>]*>.*?</nav>', lambda m: nav_html(p.name, available), t, count=1, flags=re.S)
        if new != t:
            p.write_text(new, encoding="utf-8")
    return available


def update_sitemap(build, pages):
    sm = pathlib.Path(build) / "sitemap.xml"
    t = sm.read_text(encoding="utf-8")
    for pg in pages:
        loc = f"https://trpgphysics.org/{pg}"
        if loc not in t:
            t = t.replace("</urlset>", f"  <url><loc>{loc}</loc><lastmod>2026-11-14</lastmod></url>\n</urlset>")
    sm.write_text(t, encoding="utf-8")