"""Independent radial-coordinate checks for T617; SageMath.
Run: sage -python independent_checks.py, with radial_check.py and
numerical-evidence.json beside this script. No writes to the website.

The primary solver evolves xi,q as functions of sqrt(theta). This check
independently evolves theta,q in xi, initializes with exact rational series,
and accumulates J=int xi^2 theta^(n+1) dxi. The pressure constant is recovered
from W=xi_1^3/(4*pi*(5-n)*J), using the virial identity rather than surface
mass alone. This code is adapted from T614's attached independent checks;
it shares no numerical routines with solver.py. It does not replace the
campaign's separate critique. Endpoint and truncation errors are heuristic.
"""
import json,time
from pathlib import Path
from radial_check import radial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ

# Pols (2011), Stellar Structure and Evolution, Table 4.1, printed p.47.
# Low precision normalization checks, not the source of high-precision digits.
POLS={'0':('0.119366','0.000001'),'1':('0.392699','0.000001'),
      '3/2':('0.770140','0.000001'),'2':('1.638183','0.000001'),
      '3':('11.05068','0.00001'),'4':('247.559','0.001'),
      '9/2':('4921.84','0.01')}

def main():
    expected=json.loads(Path(__file__).with_name('numerical-evidence.json').read_text())['values']
    output={}
    # Known nonzero controls run before the unfamiliar indices.
    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))
        W=x**3/(4*F.pi()*(5-nn)*J)
        surface=x**4/(4*F.pi()*(nn+1)*mass**2)
        assert abs((W-surface)/W)<F(10)**-40
        if n in ('0','1'):
            anchor=F(3)/(8*F.pi()) if n=='0' else F.pi()/8
            err=abs((W-anchor)/anchor);assert err<F(10)**-40,(n,err)
            assert abs((W*(1+F(10)**-25)-anchor)/anchor)>F(10)**-40
            meta['analytic_control_relative_error']=str(err)
            meta['displaced_control_rejected']=True
        err=abs((W-F(expected[n]['high']))/W)
        assert err<F(10)**-40,(n,err)
        meta['pressure']=str(W);meta['relative_difference_from_primary']=str(err)
        meta['surface_pressure_relative_difference']=str(abs((W-surface)/W))
        if n in POLS:
            center,ulp=map(F,POLS[n]);error=abs(W-center)/ulp
            assert error<2,('Pols printed precision',n,error)
            meta['Pols_table_4_1_error_in_last_place']=str(error)
        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()