#!/usr/bin/env python3 """SPARC rotation-curve test of the RS causal-response (ILG) gravity model. Zero per-galaxy parameters. Every global constant is either derived from phi = (1+sqrt(5))/2 or is a declared, named seam (listed below). The same masks and error model are applied to a simple-mu MOND baseline so the two models are compared on identical footing. This is a self-contained public port of the preregistered verifier `scripts/analysis/gravity_paper_v07_sparc_verify.py` (RS-derived mode, prereg tag prereg/sparc-2026-04-26) from the Recognition Physics `reality` repository. Ported to numpy-only; outputs verified identical to the original on 2026-07-08. Data: SPARC (Spitzer Photometry and Accurate Rotation Curves), Lelli, McGaugh & Schombert 2016, AJ 152, 157. Public data from http://astroweb.cwru.edu/SPARC/ . Download sparc_data.zip next to this script and unzip it (creates ./data/). Global parameter ledger (nothing else is tuned, and nothing per galaxy): DERIVED from phi: alpha = 1 - 1/phi (dynamical exponent) C_xi = 2 phi^-4 (gas-fraction coupling) A = 1 + alpha_lock/2 where alpha_lock = (1-1/phi)/2 p = 1 - alpha_lock/4 Upsilon_star = phi (mass-to-light policy) DECLARED SEAMS: tau0_seconds = 7.29968246159347e-15 (RS tick, SI bridge) N_tau = 142.4, N_r = 126.3 (galactic ladder rungs) ERROR MODEL (paper Sec. 4.3, fixed before this test, shared with MOND): sigma0 = 10 km/s, f_floor = 0.05, alpha_beam = 0.3, sigma_asym = 0.10 (dwarf) / 0.05 (spiral), k_turb = 0.07, p_turb = 1.3 MOND BASELINE: simple mu(x) = x/(1+x), a0 = 1.23e-10 m/s^2 Usage: python3 sparc_rotation_test.py [--data-dir data] """ import argparse import math import os import numpy as np trapz = getattr(np, "trapezoid", None) or np.trapz PHI = (1.0 + math.sqrt(5.0)) / 2.0 ALPHA = 1.0 - 1.0 / PHI ALPHA_LOCK = 0.5 * (1.0 - 1.0 / PHI) C_XI = 2.0 * PHI ** (-4.0) A_AMP = 1.0 + ALPHA_LOCK / 2.0 P_EXP = 1.0 - ALPHA_LOCK / 4.0 UPSILON_STAR = PHI # Declared seams (see docstring). tau0 from the reality repo's # data/calibration_tau0_seconds_gravity.json gravity-domain seam. TAU0_SECONDS = 7.29968246159347e-15 N_TAU = 142.4 N_R = 126.3 C_LIGHT = 299792458.0 KPC_TO_M = 3.086e19 KM_TO_M = 1.0e3 ELL0_M = C_LIGHT * TAU0_SECONDS TAU_STAR_S = TAU0_SECONDS * PHI ** N_TAU R0_M = ELL0_M * PHI ** N_R R0_KPC = R0_M / KPC_TO_M A0_SI = (2.0 * math.pi * R0_M) / (TAU_STAR_S ** 2) ACCEL_EXP = ALPHA * 0.5 # acceleration-space exponent is alpha/2 # Error model (fixed, shared by both models) SIGMA0 = 10.0 F_FLOOR = 0.05 ALPHA_BEAM = 0.3 SIGMA_ASYM_DWARF = 0.10 SIGMA_ASYM_SPIRAL = 0.05 K_TURB = 0.07 P_TURB = 1.3 A0_MOND = 1.23e-10 # ---------------------------------------------------------------- data loading def parse_catalog(mrt_path): """Parse the SPARC Table1 catalog (whitespace-separated token order).""" rows = [] seen = set() with open(mrt_path, "r", encoding="utf-8", errors="replace") as f: for line in f: toks = line.strip().split() if len(toks) < 18: continue try: q = int(toks[17]) except ValueError: continue if q not in (1, 2, 3): continue name = toks[0].strip() if name in seen: continue seen.add(name) def sf(s): s = s.strip() if not s: return float("nan") try: return float(s) except ValueError: return float("nan") rows.append({ "Galaxy": name, "L36": sf(toks[7]), # 1e9 Lsun at 3.6um "Rdisk_kpc": sf(toks[11]), "MHI": sf(toks[13]), # 1e9 Msun "Vflat_kms": sf(toks[15]), "Q": q, }) if not rows: raise RuntimeError(f"parsed 0 catalog rows from {mrt_path}") return rows def load_rotmod(path): """Columns: rad vobs verr vgas vdisk vbul sbdisk sbbul.""" arr = np.loadtxt(path, comments="#", ndmin=2) if arr.shape[1] < 8: pad = np.zeros((arr.shape[0], 8 - arr.shape[1])) arr = np.hstack([arr, pad]) ok = np.all(np.isfinite(arr), axis=1) arr = arr[ok] ok = (arr[:, 0] > 0) & (arr[:, 1] > 0) & (arr[:, 2] > 0) return arr[ok] # ---------------------------------------------------------------- model pieces def n_profile(r_kpc): return 1.0 + A_AMP * (1.0 - np.exp(-((r_kpc / R0_KPC) ** P_EXP))) def disk_weighted_mean_n(r_kpc, n_r, sigma_profile): w = np.maximum(sigma_profile, 0.0) * np.maximum(r_kpc, 0.0) denom = trapz(w, r_kpc) if not np.isfinite(denom) or denom <= 0: return 1.0 num = trapz(n_r * w, r_kpc) if not np.isfinite(num): return 1.0 return float(num / denom) def zeta_profile(r_kpc, Rd_kpc): if not np.isfinite(Rd_kpc) or Rd_kpc <= 0: z = np.ones_like(r_kpc) else: z = 1.0 - 0.2 * np.tanh(r_kpc / Rd_kpc) return np.clip(z, 0.8, 1.2) def gas_fraction(L36, MHI): if not np.isfinite(L36) or not np.isfinite(MHI): return float("nan") m_star = UPSILON_STAR * L36 m_gas = 1.33 * MHI denom = m_star + m_gas if denom <= 0: return float("nan") return float(m_gas / denom) def xi_factor(f_gas, thresholds): if not np.isfinite(f_gas): return 1.0 idx = int(np.searchsorted(thresholds, f_gas, side="right")) idx = max(0, min(4, idx)) u_b = idx / 4.0 return float(1.0 + C_XI * math.sqrt(u_b)) def baryon_velocity(rot): vgas, vdisk, vbul = rot[:, 3], rot[:, 4], rot[:, 5] return np.sqrt(vgas * vgas + UPSILON_STAR * (vdisk * vdisk) + UPSILON_STAR * (vbul * vbul)) def acceleration_si(v_kms, r_kpc): r_m = r_kpc * KPC_TO_M v_mps = v_kms * KM_TO_M with np.errstate(divide="ignore", invalid="ignore"): return np.where(r_m > 0, (v_mps ** 2) / r_m, np.nan) def weight_w(r_kpc, vbar, rot, Rd_kpc, xi): a_bar = np.maximum(acceleration_si(vbar, r_kpc), 1e-30) term_acc = (A0_SI / a_bar) ** ACCEL_EXP n_raw = n_profile(r_kpc) sigma_star = UPSILON_STAR * (rot[:, 6] + rot[:, 7]) # sbdisk + sbbul mean_n = disk_weighted_mean_n(r_kpc, n_raw, sigma_star) n_norm = n_raw / mean_n return xi * n_norm * term_acc * zeta_profile(r_kpc, Rd_kpc) def find_outer_cut(r_kpc, vobs, Rd_kpc): if len(r_kpc) < 5 or not np.isfinite(Rd_kpc) or Rd_kpc <= 0: return None dvdr = np.gradient(vobs, r_kpc) bad = np.where((r_kpc >= 2.0 * Rd_kpc) & (dvdr < -5.0))[0] if bad.size == 0: return None return float(r_kpc[int(bad[0])]) def sigma_eff(r_kpc, vobs, verr, Rd_kpc, sigma_asym_frac): b_kpc = 0.3 * Rd_kpc s_beam = ALPHA_BEAM * (b_kpc * vobs) / (r_kpc + b_kpc) s_turb = K_TURB * vobs * (1.0 - np.exp(-r_kpc / max(Rd_kpc, 1e-9))) ** P_TURB return np.sqrt(verr ** 2 + SIGMA0 ** 2 + (F_FLOOR * vobs) ** 2 + s_beam ** 2 + (sigma_asym_frac * vobs) ** 2 + s_turb ** 2) def galaxy_type(vmax): if not np.isfinite(vmax): return "unknown" if vmax < 80: return "dwarf" if vmax <= 200: return "spiral" return "massive" def mond_velocity(vbar, r_kpc): a_bar = np.maximum(acceleration_si(vbar, r_kpc), 1e-30) a = 0.5 * (a_bar + np.sqrt(a_bar * a_bar + 4.0 * A0_MOND * a_bar)) return np.sqrt(np.maximum(a, 0.0) * (r_kpc * KPC_TO_M)) / KM_TO_M # ---------------------------------------------------------------- main analysis CHART_GALAXY = "NGC3198" def main(): ap = argparse.ArgumentParser(description="SPARC test of RS causal-response gravity") ap.add_argument("--data-dir", default="data") args = ap.parse_args() mrt = os.path.join(args.data_dir, "SPARC_Lelli2016c.mrt.txt") rotdir = os.path.join(args.data_dir, "rotmod") if not os.path.isfile(mrt) or not os.path.isdir(rotdir): raise SystemExit( f"data not found under '{args.data_dir}'. Download sparc_data.zip " "next to this script and unzip it (creates ./data/).") catalog = parse_catalog(mrt) print("=" * 74) print("SPARC ROTATION-CURVE TEST: RS causal-response model, zero per-galaxy knobs") print("=" * 74) print(f"phi = {PHI:.12f}") print(f"derived: alpha = {ALPHA:.6f} C_xi = {C_XI:.6f} " f"A = {A_AMP:.6f} p = {P_EXP:.6f} Upsilon* = phi") print(f"seams: tau0 = {TAU0_SECONDS:.6e} s N_tau = {N_TAU} N_r = {N_R}") print(f" -> r0 = {R0_KPC:.3f} kpc a0 = {A0_SI:.4e} m/s^2") print(f"baseline: MOND simple mu, a0 = {A0_MOND:.2e} m/s^2 (same masks, same errors)") # gas-fraction quintiles from the Q=1 sample fq1 = np.array([gas_fraction(r["L36"], r["MHI"]) for r in catalog if r["Q"] == 1]) fq1 = fq1[np.isfinite(fq1)] thresholds = np.quantile(fq1, [0.2, 0.4, 0.6, 0.8]) per_gal = [] rar_residuals = [] vflat_rotmod = {} chart_lines = [] for row in catalog: name = row["Galaxy"] path = os.path.join(rotdir, f"{name}_rotmod.dat") if not os.path.isfile(path): continue rot = load_rotmod(path) if rot.shape[0] == 0: continue r_kpc, vobs, verr = rot[:, 0], rot[:, 1], rot[:, 2] vmax = float(np.max(vobs)) Rd = row["Rdisk_kpc"] if not np.isfinite(Rd) or Rd <= 0: vdisk = rot[:, 4] Rd = float(r_kpc[int(np.argmax(vdisk))] / 2.2) if np.any(vdisk > 0) else 2.0 xi = xi_factor(gas_fraction(row["L36"], row["MHI"]), thresholds) vbar = baryon_velocity(rot) w = weight_w(r_kpc, vbar, rot, Rd, xi) v_model = np.sqrt(np.maximum(w, 0.0)) * vbar mask = (r_kpc >= 0.3 * Rd) & ((vobs / verr) >= 3.0) r_cut = find_outer_cut(r_kpc, vobs, Rd) if r_cut is not None: mask &= r_kpc <= r_cut if mask.sum() < 1: continue gtype = galaxy_type(vmax) s_asym = SIGMA_ASYM_DWARF if gtype == "dwarf" else SIGMA_ASYM_SPIRAL sig = sigma_eff(r_kpc, vobs, verr, Rd, s_asym) chi2N = float(np.sum(((vobs - v_model)[mask] / sig[mask]) ** 2) / mask.sum()) v_mond = mond_velocity(vbar, r_kpc) chi2N_m = float(np.sum(((vobs - v_mond)[mask] / sig[mask]) ** 2) / mask.sum()) per_gal.append((name, row["Q"], gtype, int(mask.sum()), chi2N, chi2N_m)) if row["Q"] == 1: a_obs = acceleration_si(vobs, r_kpc)[mask] a_mod = acceleration_si(v_model, r_kpc)[mask] ok = np.isfinite(a_obs) & np.isfinite(a_mod) & (a_obs > 0) & (a_mod > 0) if np.any(ok): rar_residuals.extend((np.log10(a_obs[ok]) - np.log10(a_mod[ok])).tolist()) vv = vobs[mask] if vv.size: vflat_rotmod[name] = float(np.mean(vv[-3:])) if name == CHART_GALAXY: chart_lines.append(f"\nchart data: {name} (Q={row['Q']}, Rd = {Rd:.2f} kpc)") chart_lines.append(f"{'r_kpc':>7} {'vobs':>7} {'verr':>6} {'vbar':>7} " f"{'v_RS':>7} {'v_MOND':>7} {'masked':>6}") for i in range(len(r_kpc)): chart_lines.append( f"{r_kpc[i]:7.2f} {vobs[i]:7.1f} {verr[i]:6.1f} {vbar[i]:7.1f} " f"{v_model[i]:7.1f} {v_mond[i]:7.1f} {str(bool(mask[i])):>6}") # ---------------- aggregates over the Q=1 (highest-quality) subsample q1 = [g for g in per_gal if g[1] == 1] chi_rs = np.array([g[4] for g in q1]) chi_mond = np.array([g[5] for g in q1]) n_points = int(sum(g[3] for g in q1)) print(f"\nper-galaxy chi^2/N over the Q=1 subsample ({len(q1)} galaxies, " f"{n_points} masked points):") print(f"{'galaxy':<14}{'type':<9}{'N':>4} {'RS chi2/N':>11} {'MOND chi2/N':>12}") for name, _q, gtype, n, c_rs, c_m in q1: print(f"{name:<14}{gtype:<9}{n:>4} {c_rs:>11.3f} {c_m:>12.3f}") med_rs = float(np.median(chi_rs)) med_m = float(np.median(chi_mond)) mean_rs = float(np.mean(chi_rs)) f_good = float(np.mean(chi_rs < 2.0)) print("\n" + "=" * 74) print("AGGREGATE (Q=1)") print("=" * 74) print(f"median chi^2/N RS causal-response: {med_rs:.4f} MOND: {med_m:.4f}") print(f"mean chi^2/N RS causal-response: {mean_rs:.4f} " f"MOND: {float(np.mean(chi_mond)):.4f}") print(f"fraction of galaxies with RS chi^2/N < 2: {f_good:.4f}") print(f"galaxies with chi^2/N > 5: RS {int(np.sum(chi_rs > 5))} " f"MOND {int(np.sum(chi_mond > 5))}") rar = np.array(rar_residuals) print(f"\nRAR scatter (log10 a_obs - log10 a_model, {rar.size} points): " f"{float(np.std(rar)):.4f} dex") # BTFR from rotmod-defined vflat (mean of outermost 3 masked points) xs, ys = [], [] for row in catalog: if row["Q"] != 1 or row["Galaxy"] not in vflat_rotmod: continue mb = UPSILON_STAR * row["L36"] * 1e9 + 1.33 * row["MHI"] * 1e9 vf = vflat_rotmod[row["Galaxy"]] if np.isfinite(mb) and mb > 0 and np.isfinite(vf) and vf > 0: xs.append(math.log10(vf)) ys.append(math.log10(mb)) if len(xs) >= 10: slope, intercept = np.polyfit(np.array(xs), np.array(ys), 1) scatter = float(np.std(np.array(ys) - (slope * np.array(xs) + intercept))) print(f"BTFR ({len(xs)} galaxies): slope = {slope:.4f}, " f"scatter = {scatter:.4f} dex") # preregistered decision gate (fixed before the test; reported verbatim) if mean_rs < 2.0 and f_good > 0.7: verdict = "PASS" elif mean_rs > 2.0 or f_good < 0.5: verdict = "FAIL" else: verdict = "INCONCLUSIVE" print("\npreregistered gate: PASS iff mean chi^2/N < 2.0 AND f(chi^2/N<2) > 0.7") print(f" FAIL iff mean chi^2/N > 2.0 OR f(chi^2/N<2) < 0.5") print(f"VERDICT: {verdict} (mean = {mean_rs:.3f}, f_good = {f_good:.3f})") for line in chart_lines: print(line) if __name__ == "__main__": main()