The file src/32-dark/rar_shape.py of the validation package of the paper 32-dark, rendered with an anchor per line; the file itself is beside this page.
1#!/usr/bin/env python32"""3rar_shape.py -- the named test of row B8: the approach to Newton against the binned radial acceleration relation.45The paper's interpolation is g_obs = g_b / (1 - e^{-sqrt(x)}), x = g_b/a0 (the first-passage registration with the6amplitude barrier, row B8); the alternative with the same two limits is the 'simple' rational form7g_obs = g_b * (1 + sqrt(1 + 4/x))/2. The two differ by up to 0.05 in g_obs/g_b near g_b ~ 5 a0 (interpolation.py).8This script bins the 2693 SPARC points of Lelli et al. (2017, data/RAR.mrt) in log g_b (0.2 dex), takes the mean9log g_obs per bin with the bin's scatter (the standard deviation, the error bars of the binned relation) as its10uncertainty, and confronts both forms with the binned relation: (i) at the fixed11floor a0 = cH0/2pi (H0 = 67.4), (ii) with a0 free for each form (a grid search on the binned chi^2). A labelled12[approx] data comparison; the exact statements of the paper are elsewhere (firstpassage_finite.py, born_exact.py).13"""14import numpy as np, os15HERE = os.path.dirname(os.path.abspath(__file__))16c = 2.99792458e8; H0 = 67.4e3/3.0857e22; a0_floor = c*H0/(2*np.pi)17fails = 018def check(name, ok):19 global fails20 print((" PASS " if ok else " FAIL ") + name)21 if not ok: fails += 122def nu_exp(x): return 1/(1-np.exp(-np.sqrt(x)))23def nu_simple(x): return 0.5*(1+np.sqrt(1+4/x))2425d = np.loadtxt(os.path.join(HERE, "data", "RAR.mrt"), skiprows=14)26lgb, egb, lgo, ego = d.T27print(f"{len(lgb)} SPARC points; a0 = cH0/2pi = {a0_floor:.4e} m/s^2 (H0 = 67.4); the RAR fit g_dagger = 1.20e-10")2829# --- the binned relation: mean log g_obs in 0.2-dex bins of log g_b, the bin's scatter as its uncertainty30edges = np.arange(-12.0, -8.4, 0.2)31bins = []32for lo, hi in zip(edges[:-1], edges[1:]):33 m = (lgb >= lo) & (lgb < hi)34 if m.sum() < 10: continue35 mean = lgo[m].mean(); sd = lgo[m].std() # the bin's mean and scatter (the binned relation's error bar)36 bins.append((0.5*(lo+hi), m.sum(), mean, sd))37bins = np.array(bins); cb, nb, mb, eb = bins.T38print(f"{len(bins)} bins of 0.2 dex in log g_b with >= 10 points; the uncertainty of a bin is its scatter (0.03-0.13 dex)")3940def chi2(form, a0, sel=None):41 pred = np.log10(10**cb * form(10**cb/a0))42 r = (mb - pred)/eb43 return np.sum(r[sel]**2) if sel is not None else np.sum(r**2), pred4445print("\n== 1. both forms at the fixed floor a0 = cH0/2pi")46print(f" {'log g_b':>8} {'N':>4} {'<log g_obs>':>12} {'sd':>6} {'exp':>8} {'simple':>8} {'res exp':>8} {'res smp':>8}")47c_e0, p_e0 = chi2(nu_exp, a0_floor); c_s0, p_s0 = chi2(nu_simple, a0_floor)48for i in range(len(cb)):49 print(f" {cb[i]:8.1f} {int(nb[i]):4d} {mb[i]:12.3f} {eb[i]:6.3f} {p_e0[i]:8.3f} {p_s0[i]:8.3f} {mb[i]-p_e0[i]:+8.3f} {mb[i]-p_s0[i]:+8.3f}")50print(f" chi^2 over {len(cb)} bins: exponential {c_e0:.1f}, simple {c_s0:.1f}")5152print("\n== 2. a0 free for each form (grid 0.8-1.6e-10)")53grid = np.linspace(0.8e-10, 1.6e-10, 801)54best = {}55for name, form in (("exponential", nu_exp), ("simple", nu_simple)):56 cs = np.array([chi2(form, a)[0] for a in grid]); k = int(np.argmin(cs))57 best[name] = (grid[k], cs[k]); print(f" {name:12s}: best a0 = {grid[k]:.3e}, chi^2 = {cs[k]:.1f}")58a_e, c_e = best["exponential"]; a_s, c_s = best["simple"]59check("the exponential form's best-fit a0 lies within 10 % of the RAR fit g_dagger = 1.20e-10", abs(a_e/1.20e-10 - 1) < 0.10)60check("with a0 free, the exponential form fits the binned relation better than the simple form (lower chi^2)", c_e < c_s)6162print("\n== 3. the discriminating band 2 < x < 10 (g_b ~ 5 a0), each form at its own best a0")63knee_e = (10**cb/a_e > 2) & (10**cb/a_e < 10); knee_s = (10**cb/a_s > 2) & (10**cb/a_s < 10)64ce_k, pe = chi2(nu_exp, a_e, knee_e); cs_k, ps = chi2(nu_simple, a_s, knee_s)65print(f" {'log g_b':>8} {'<log g_obs>':>12} {'sd':>6} {'exp':>8} {'simple':>8} {'res exp':>8} {'res smp':>8}")66for i in np.where(knee_e | knee_s)[0]:67 print(f" {cb[i]:8.1f} {mb[i]:12.3f} {eb[i]:6.3f} {pe[i]:8.3f} {ps[i]:8.3f} {mb[i]-pe[i]:+8.3f} {mb[i]-ps[i]:+8.3f}")68print(f" chi^2 in the band: exponential {ce_k:.1f} ({knee_e.sum()} bins), simple {cs_k:.1f} ({knee_s.sum()} bins)")69print(f" the forms differ by up to {np.max(np.abs(np.log10(nu_simple(10**cb/a_e)/nu_exp(10**cb/a_e)))):.3f} dex over the bins")70check("in the band 2 < x < 10 the binned relation departs from the exponential form by less than the bin scatter in every bin (the B8 falsifier not triggered)",71 np.all(np.abs((mb - pe)[knee_e]) < eb[knee_e]))72off = (mb - pe)[knee_e]73print(f" note: the binned relation sits {-off.mean():.3f} dex below the exponential form on average in the band, within the scatter")74print(f" power: the two forms differ by <= {np.max(np.abs(np.log10(nu_simple(10**cb/a_e)/nu_exp(10**cb/a_e)))):.3f} dex per bin against a bin scatter of {eb[knee_e].mean():.3f} dex: the binned relation prefers the exponential form by chi^2 but does not exclude the rational form")75check("the simple rational form is disfavoured in the band (its chi^2 exceeds the exponential form's)", cs_k > ce_k)7677print("\n== 4. the fixed floor a0 = cH0/2pi against the free fit")78print(f" floor/fit = {a0_floor/a_e:.3f} (the paper: 13 % below the fitted value); Delta chi^2 (fixed - free) = {c_e0 - c_e:.1f}")79check("the free-fit a0 of the exponential form exceeds the floor cH0/2pi by 5-25 % (the paper's 13 %)", 1.05 < a_e/a0_floor < 1.25)80print("\nrar_shape.py:", "PASS" if fails == 0 else f"FAIL ({fails})")81raise SystemExit(1 if fails else 0)82