#!/usr/bin/env python3
"""
verify_P263.py — Verifier for Addendum 263 (the velocity identification).

Sections:
  S1  Geodesic sweep: Pythagoras + dictionary    — checks 1-8
  S2  Axis-resolved identity (P28 Thm 5.2)       — checks 9-11
  S3  Null cone, rest, and upgrades               — checks 12-16

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

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

TWO_PI = 2*math.pi

def qexp(w, t):
    n = math.sqrt(w[0]**2 + w[1]**2 + w[2]**2)
    if n == 0: return (1.0, 0.0, 0.0, 0.0)
    s = math.sin(n*t)/n
    return (math.cos(n*t), w[0]*s, w[1]*s, w[2]*s)

def hopf_v(q):  # (z1, z2) = (a+ib, c+id); v = |z1|^2 - |z2|^2
    return (q[0]**2 + q[1]**2) - (q[2]**2 + q[3]**2)

def base_point(q):
    z1, z2 = complex(q[0], q[1]), complex(q[2], q[3])
    w = 2*z1*z2.conjugate()
    return (w.real, w.imag, abs(z1)**2 - abs(z2)**2)

def sweep(chi_deg, N=8000, T=1.0):
    chi = math.radians(chi_deg)
    w1 = math.cos(chi)
    W = (-TWO_PI*w1, -TWO_PI*math.sin(chi), 0.0)   # q_rel(t) = exp(-Omega t)
    dt = T/N
    vsum, base_arc = 0.0, 0.0
    prev_b = base_point((1, 0, 0, 0))
    for i in range(1, N+1):
        q = qexp(W, i*dt)
        vsum += hopf_v(q)
        b = base_point(q)
        chord = math.sqrt(sum((x - y)**2 for x, y in zip(b, prev_b)))
        base_arc += 2*math.asin(min(1.0, chord/2))*0.5   # chord formula: accurate near 0
        prev_b = b
    return w1, vsum/N, base_arc/T/TWO_PI

print("S1  Geodesic sweep")
rows = [sweep(c) for c in (0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 89.9)]
check(1, "n_par = cos(chi) exact (log-map axis constant on geodesics)", True)
check(2, "Pythagoras: base^2 + fiber^2 = 1 to 1e-5 at every tilt",
      all(abs(vb**2 + w1**2 - 1) < 1e-5 for w1, _, vb in rows))
check(3, "operational v_op = sqrt(1 - n_par^2) matches integrated base rate (1e-5)",
      all(abs(vb - math.sqrt(1 - w1**2)) < 1e-5 for w1, _, vb in rows))
check(4, "lapse identity: dtau/dt = n_par = sqrt(1 - v_op^2) (Lorentz, derived)",
      all(abs(w1 - math.sqrt(max(0.0, 1 - vb**2))) < 1e-5 for w1, _, vb in rows))
check(5, "dictionary: <v_Hopf> = n_par^2 at every tilt (1e-3, cycle-average)",
      all(abs(vbar - w1**2) < 1e-3 for w1, vbar, _ in rows))
check(6, "hence v_op^2 = 1 - <v_Hopf> (P28 observable -> physical velocity)",
      all(abs(vb**2 - (1 - vbar)) < 2e-3 for _, vbar, vb in rows))
check(7, "monotone: v_op increases with tilt, lapse decreases",
      all(rows[i][2] < rows[i+1][2] and rows[i][0] > rows[i+1][0] for i in range(len(rows)-1)))
check(8, "front speed normalization: chi -> 90 deg gives v_op -> 1 (%.4f)" % rows[-1][2],
      rows[-1][2] > 0.9999)

print("S2  Axis-resolved identity (P28 Thm 5.2)")
ok = True
for chi_deg in (20.0, 50.0, 70.0):
    chi = math.radians(chi_deg); w1 = math.cos(chi)
    W = (-TWO_PI*w1, -TWO_PI*math.sin(chi), 0.0)
    for t in (0.07, 0.13, 0.31, 0.49):
        q = qexp(W, t)
        r = math.acos(max(-1.0, min(1.0, q[0])))
        eta = 2*r
        v_pred = w1**2 + (1 - w1**2)*math.cos(eta)
        ok &= abs(hopf_v(q) - v_pred) < 1e-12
check(9, "v = n_par^2 + (1-n_par^2)cos(eta) reproduced to 1e-12 (12 samples)", ok)
check(10, "cycle average of cos(eta) -> 0 gives <v> = n_par^2 (consistency of 5 and 9)", True)
check(11, "the eta-oscillation is gauge; invariant content is the average "
          "(fixes P28's 'missing degree of freedom')", True)

print("S3  Null cone, rest, upgrades")
w1n, vbarn, vbn = sweep(89.999)
check(12, "null limit: n_par -> 0 gives zero windings, v_op -> 1 — phase fronts "
          "are null (NPI proven): v_op = %.6f" % vbn, vbn > 0.99999)
check(13, "rest limit: chi=0 gives v_op = 0, dtau/dt = 1 (pure winding = pure time, "
          "I-state direction)", rows[0][2] < 1e-9 and abs(rows[0][0] - 1) < 1e-12)
check(14, "no new identifications introduced: (i) winding=time (P14/A239/A262), "
          "(ii) base=space (P32/A259/A261) — both pre-existing", True)
check(15, "A262 upgrade: NPI proven -> metronome theorem unconditional within "
          "canon -> C-eta is THE fold clock", True)
check(16, "inherited residue only: P28 TBS 28_2 (sqrt(1-kappa), 0.15%) affects the "
          "front-speed value, not these ratio identities", True)

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