#!/usr/bin/env python3
"""
verify.py — the check behind "The Number That Ends the Argument" (a combine portal).

The portal's claim is methodological: a folk/schoolbook explanation is settled not by
argument but by *sizing* the named mechanism from its own governing law and comparing
the magnitude to what the observation requires. Sizing returns one of three verdicts:

  (A) FALLS ORDERS SHORT   — the mechanism is real but 10^1..10^5 too weak.
  (B) WRONG QUANTITY/CAUSE — the felt effect is real but the named variable isn't it.
  (C) CONFIRMED            — sizing endorses the (correct) mechanism.

This file recomputes every number the portal shows live, from first principles, and
asserts each matches the figure stated in the member stratum it summarizes. The portal's
in-browser JS uses the SAME formulas; if they ever disagree, one of them is lying.

Run:  python3 verify.py        (prints a table; exits non-zero on any mismatch)

No third-party deps. Physical constants are SI; sources are named per block.
"""

import math

PASS, FAILS = [], []

def check(name, got, want, tol_rel=0.05, unit="", note=""):
    """Assert `got` is within tol_rel of `want`; record the result."""
    ok = (want == 0 and abs(got) < 1e-12) or abs(got - want) <= tol_rel * abs(want)
    (PASS if ok else FAILS).append(name)
    flag = "ok " if ok else "XX "
    print(f"  [{flag}] {name:<46} got {got:>12.4g} {unit:<8} (want ~{want:g} {unit}) {note}")
    return ok

def band(title):
    print(f"\n== {title} ==")

# ----------------------------------------------------------------------------
band("1. ICE / pressure-melting  —  the-ice-that-pressure-didnt-melt  (A: falls short)")
# Clausius-Clapeyron on the ice-water line:  dT/dP = T*Dv / L
# Dv = v_liquid - v_solid (per kg); ice is LESS dense, so Dv < 0 -> negative slope.
T_melt = 273.15                       # K
rho_water, rho_ice = 999.8, 916.7     # kg/m^3
Dv = 1.0/rho_water - 1.0/rho_ice      # m^3/kg  (negative)
L_fus = 333.55e3                      # J/kg, latent heat of fusion of ice
dTdP = T_melt * Dv / L_fus            # K/Pa
check("Clausius-Clapeyron slope", dTdP*1e6, -0.074, unit="K/MPa", tol_rel=0.05)

m_skater, g = 70.0, 9.81              # kg, m/s^2
A_blade = 1e-4                        # m^2  (1 cm^2, a generous single-blade contact patch)
P_blade = m_skater*g / A_blade        # Pa
check("skater blade pressure", P_blade/1e6, 6.9, unit="MPa")

dT_press = dTdP * P_blade             # K, melting-point depression delivered by pressure
check("melting-point depression", abs(dT_press), 0.51, unit="K", note="~half a degree")

rink_below = 5.0                      # K below freezing (a mild rink; hard rinks are 5-10)
miss_ice = abs(dT_press) / rink_below # dimensionless miss factor M = delivered/needed
check("miss factor M (vs -5 C rink)", miss_ice, 0.102, note=f"~{1/miss_ice:.0f}x short")

# ----------------------------------------------------------------------------
band("2. CORIOLIS / the draining sink  —  coriolis-drain  (A: falls short, extreme)")
# Coriolis parameter f = 2 Omega sin(phi); its acceleration on the flow is a_c = f*U.
# A visible drain vortex has centripetal acceleration a_vortex = U^2 / r.
Omega = 7.2921e-5                     # rad/s, Earth's sidereal rotation
phi = math.radians(42.0)             # latitude (Boston-ish)
f_cor = 2*Omega*math.sin(phi)        # 1/s
U_drain, r_drain = 0.3, 0.05         # m/s, m  (a lively sink vortex)
a_cor = f_cor * U_drain              # m/s^2
a_vortex = U_drain**2 / r_drain      # m/s^2
check("Coriolis acceleration a_c", a_cor, 3e-5, unit="m/s^2", tol_rel=0.35)
check("vortex centripetal accel", a_vortex, 1.8, unit="m/s^2")
miss_cor = a_cor / a_vortex          # M
check("miss factor M", miss_cor, 1.6e-5, tol_rel=0.4, note=f"~{1/miss_cor:.0f}x short")
# The timescale form the portal also shows: 1/f is HOURS, a sink drains in SECONDS.
t_cor = 1.0/f_cor
check("Coriolis time 1/f (hours)", t_cor/3600, 2.85, unit="hr", tol_rel=0.1,
      note="Shapiro 1962 had to settle water a full DAY to see it")

# ----------------------------------------------------------------------------
band("3. SALT / the boil  —  the-salt-that-barely-moves-the-boil  (A short + WRONG SIGN)")
# Boiling-point ELEVATION (colligative):  dTb = i * Kb * b.  Note: sign is +,
# i.e. salt makes water boil at a HIGHER temperature -> the folk 'boils faster' is backwards.
i_vanthoff, Kb = 2.0, 0.512          # NaCl dissociates ~2; ebullioscopic constant K.kg/mol
dose_g_per_L = 10.0                  # a heaped teaspoon per litre (a heavy hand)
M_NaCl = 58.44                       # g/mol
b_molal = (dose_g_per_L / M_NaCl) / 1.0   # mol per kg water (1 L ~ 1 kg)
dTb = i_vanthoff * Kb * b_molal      # K, POSITIVE = elevation
check("boiling-point elevation dTb", dTb, 0.175, unit="K", tol_rel=0.15,
      note="under 0.2 C; and it RAISES the boil, opposite of 'faster'")
assert dTb > 0, "salt raises the boiling point; sign must be positive"

# ----------------------------------------------------------------------------
band("4. METAL vs WOOD / effusivity  —  metal-colder-than-wood  (B: wrong quantity)")
# Two semi-infinite bodies in contact settle to  T_c = (e1 T1 + e2 T2)/(e1+e2),
# where effusivity e = sqrt(k*rho*c). Every object here is at the SAME temperature (20 C):
# the felt 'cold' is a heat-loss RATE (effusivity), not a temperature difference.
def effusivity(k, rho, c): return math.sqrt(k*rho*c)
# Skin effusivity is quoted, not derived: skin's k,rho,c vary with perfusion/moisture,
# and the literature value used by the member stratum is ~1500 Ws^.5/m^2K. We adopt it.
e_skin  = 1500.0                         # human skin (literature; ~1000-1700 range)
e_steel = effusivity(50.0, 7850, 490)    # mild steel
e_oak   = effusivity(0.17, 700, 1700)    # oak, across grain
e_foam  = effusivity(0.033, 30, 1300)    # expanded polystyrene
check("effusivity steel", e_steel, 13500, unit="Ws^.5/m^2K", tol_rel=0.1)
check("effusivity oak",   e_oak,   450,   unit="Ws^.5/m^2K", tol_rel=0.2)

T_skin, T_obj = 33.0, 20.0           # fingertip, room objects — all objects identical at 20 C
def contact(e_obj): return (e_skin*T_skin + e_obj*T_obj)/(e_skin+e_obj)
check("contact temp on steel", contact(e_steel), 21.3, unit="C", tol_rel=0.05)
check("contact temp on oak",   contact(e_oak),   30.0, unit="C", tol_rel=0.05)
check("contact temp on foam",  contact(e_foam),  32.7, unit="C", tol_rel=0.05)
dT_object = T_obj - T_obj            # the folk 'cause' magnitude: the objects' temp difference
check("object temperature difference", dT_object, 0.0, unit="K",
      note="folk cause is exactly zero -> the miss isn't SIZE, it's KIND")

# ----------------------------------------------------------------------------
band("5. ROAST / carryover  —  the-roast-keeps-cooking  (C: CONFIRMED)")
# Sizing is not only a debunking tool. Thermal penetration depth delta ~ sqrt(alpha*t):
# if heat reaches the centre during the rest, the centre keeps cooking. It does.
alpha_meat = 1.2e-7                   # m^2/s, thermal diffusivity of muscle
for t_min in (10, 20):
    delta = math.sqrt(alpha_meat * t_min*60) * 100  # cm
    print(f"  [ -- ] penetration depth after {t_min:>2} min rest      delta = {delta:5.2f} cm")
delta_20 = math.sqrt(alpha_meat * 20*60) * 100
check("penetration depth @20min", delta_20, 1.2, unit="cm", tol_rel=0.1,
      note="~1 cm reaches a small roast's core -> carryover CONFIRMED (+3..14 C, per source)")
miss_roast = 1.0                     # lands where the observation is
check("miss factor M (confirmed ~1x)", miss_roast, 1.0, note="sizing endorses the mechanism")

# ----------------------------------------------------------------------------
band("SYNTHESIS  —  the three verdicts on one axis")
rows = [
    ("Ice / pressure-melting", "A short",  miss_ice,  f"{1/miss_ice:.0f}x short"),
    ("Coriolis / the sink",    "A short",  miss_cor,  f"{1/miss_cor:.0f}x short"),
    ("Salt / the boil",        "A short",  None,      "negligible + wrong sign"),
    ("Metal vs wood",          "B kind",   None,      "wrong quantity (dT = 0)"),
    ("Roast / carryover",      "C conf.",  miss_roast, "confirmed (~1x)"),
]
print("  myth                        verdict   log10(M)     reading")
for name, verdict, M, reading in rows:
    lm = f"{math.log10(M):+6.2f}" if M else "   -- "
    print(f"    {name:<26}{verdict:<9}{lm:>8}    {reading}")

# ----------------------------------------------------------------------------
print("\n" + "="*70)
print(f"  {len(PASS)} checks passed, {len(FAILS)} failed.")
if FAILS:
    print("  FAILED:", ", ".join(FAILS))
    raise SystemExit(1)
print("  All numbers reproduce, and match the figures in each member stratum.")
print("="*70)
