back to table · edit · history · where entries came from · files · download
2578 bytes, as of the version from 2026-10-30 05:29 (current). Recorded here, not run.
"""Independent radial-coordinate checks for T618; run with SageMath.
Run: sage -python independent_checks.py, with radial_check.py and
numerical-evidence.json beside it. The primary solver evolves xi,q in
sqrt(theta). This check independently evolves theta,q in xi, initializes
with exact rational coefficients, and integrates J=int xi^2 theta^(n+1) dxi.
It recovers mass from the virial identity J=(n+1)*mu^2/((5-n)*xi_1), then
uses the physical scale equations with a=rho_c=G=1 to recover N_n.
No routines from the primary solver are imported. This is build validation,
not the separate campaign critique. Endpoint/truncation errors are heuristic.
"""
import json,time
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from radial_check import radial
# Pols (2011), Table 4.1, printed p.47. Normalization checks only.
POLS={'1':'0.63662','3/2':'0.42422','2':'0.36475',
'3':'0.36394','4':'0.47720','9/2':'0.65798'}
def main():
expected=json.loads(Path(__file__).with_name('numerical-evidence.json').read_text())['values']
output={}
# Nonzero known answers before any unfamiliar model, including a negative control.
for n in ('0','1','1/2','3/2','2','5/2','3','7/2','4','9/2'):
started=time.monotonic();mass,meta=radial(n);F=mass.parent()
x=F(meta['radius']);J=F(meta['pressure_integral']);nn=F(QQ(n))
qJ=((5-nn)*x*J/(nn+1)).sqrt()
if n in ('0','1'):
xa=F(6).sqrt() if n=='0' else F.pi()
qa=2*F(6).sqrt() if n=='0' else F.pi()
err=max(abs(x/xa-1),abs(qJ/qa-1));assert err<F(10)**-40
assert abs(qJ*(1+F(10)**-25)/qa-1)>F(10)**-40
meta['analytic_background_control_error']=str(err)
meta['displaced_control_rejected']=True
if n!='0':
K=4*F.pi()/(nn+1);M=4*F.pi()*qJ;R=x
value=K/(M**((nn-1)/nn)*R**((3-nn)/nn))
if n=='1':assert abs(value/(2/F.pi())-1)<F(10)**-40
err=abs(value/F(expected[n]['coefficient_210'])-1)
assert err<F(10)**-40,(n,err)
meta['coefficient']=str(value);meta['relative_difference_from_primary']=str(err)
if n in POLS:
ulps=abs(value-F(POLS[n]))/F('0.00001')
assert ulps<1,('Pols Table 4.1',n,ulps)
meta['Pols_error_in_last_place']=str(ulps)
meta['seconds']=time.monotonic()-started;output[n]=meta
print('CHECK '+n+' '+json.dumps(meta),flush=True)
print('DATA '+json.dumps(output),flush=True)
if __name__=='__main__':main()