"""Independent numerical evolution of actual JavaScript Gyr snapshots."""
from pathlib import Path
import json
import numpy as np
from scipy.integrate import solve_ivp
native=json.loads(Path('toy-rotation-reference.json').read_text())
errors=[]
for row in native['cases']:
    initial=2*np.pi/row['initialPeriodDays'];saturation=2*np.pi/3
    if row['ageGyr']==0:period=row['initialPeriodDays']
    else:
        solution=solve_ivp(lambda t,y:[-row['coefficient']*y[0]*min(y[0]**2,saturation**2)],[0,row['ageGyr']],[initial],method='DOP853',rtol=1e-11,atol=1e-12)
        assert solution.success
        period=2*np.pi/solution.y[0,-1]
    errors.append(abs(period/row['periodDays']-1))
assert max(errors)<1e-8
report={'status':'PASS40 actual JavaScript evolutionary snapshots versus independent differential-equation integration','maximum_relative_period_error':max(errors),'units':{'observingTime':'days','evolutionaryTime':'Gyr','coefficient':'day squared/Gyr'},'qualification':'Prescribed fixed inertia and saturation; no mass-dependent stellar structure evolution or empirical age inference'}
Path('results/evolution-reference-results.json').write_text(json.dumps(report,indent=2)+'\n')
print(json.dumps(report,indent=2))